VERB_code_2.3
dgematrix-lapack.hpp
1 //=============================================================================
6 inline long dgematrix::dgesv(dgematrix& mat)
7 {
8 #ifdef CPPL_VERBOSE
9  std::cerr << "# [MARK] dgematrix::dgesv(dgematrix&)"
10  << std::endl;
11 #endif//CPPL_VERBOSE
12 
13 #ifdef CPPL_DEBUG
14  if(M!=N || M!=mat.M){
15  std::cerr << "[ERROR] dgematrix::dgesv(dgematrix&) " << std::endl
16  << "These two matrices cannot be solved." << std::endl
17  << "Your input was (" << M << "x" << N << ") and ("
18  << mat.M << "x" << mat.N << ")." << std::endl;
19  exit(1);
20  }
21 #endif//CPPL_DEBUG
22 
23  long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.M), INFO(1);
24  dgesv_(N, NRHS, Array, LDA, IPIV, mat.Array, LDB, INFO);
25  delete [] IPIV;
26 
27  if(INFO!=0){
28  std::cerr << "[WARNING] dgematrix::dgesv(dgematrix&) "
29  << "Serious trouble happend. INFO = " << INFO << "."
30  << std::endl;
31  }
32  return INFO;
33 }
34 
35 //=============================================================================
40 inline long dgematrix::dgesv(dcovector& vec)
41 {
42 #ifdef CPPL_VERBOSE
43  std::cerr << "# [MARK] dgematrix::dgesv(dcovector&)"
44  << std::endl;
45 #endif//CPPL_VERBOSE
46 
47 #ifdef CPPL_DEBUG
48  if(M!=N || M!=vec.L){
49  std::cerr << "[ERROR] dgematrix::dgesv(dcovector&) " << std::endl
50  << "These matrix and vector cannot be solved." << std::endl
51  << "Your input was (" << M << "x" << N << ") and ("
52  << vec.L << ")." << std::endl;
53  exit(1);
54  }
55 #endif//CPPL_DEBUG
56 
57  long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), INFO(1);
58  dgesv_(N, NRHS, Array, LDA, IPIV, vec.Array, LDB, INFO);
59  delete [] IPIV;
60 
61  if(INFO!=0){
62  std::cerr << "[WARNING] dgematrix::dgesv(dcovector&) "
63  << "Serious trouble happend. INFO = " << INFO << "."
64  << std::endl;
65  }
66  return INFO;
67 }
68 
72 
73 //=============================================================================
76 inline long dgematrix::dgels(dgematrix& mat)
77 {
78 #ifdef CPPL_VERBOSE
79  std::cerr << "# [MARK] dgematrix::dgels(dgematrix&)"
80  << std::endl;
81 #endif//CPPL_VERBOSE
82 
83 #ifdef CPPL_DEBUG
84  if(M!=mat.M){
85  std::cerr << "[ERROR] dgematrix::dgels(dgematrix&) " << std::endl
86  << "These two matrices cannot be solved." << std::endl
87  << "Your input was (" << M << "x" << N << ") and ("
88  << mat.M << "x" << mat.N << ")." << std::endl;
89  exit(1);
90  }
91 #endif//CPPL_DEBUG
92 
93  if(M<N){ //underdetermined
94  dgematrix tmp(N,mat.N);
95  for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
96  tmp(i,j) =mat(i,j);
97  }}
98  mat.clear();
99  swap(mat,tmp);
100  }
101 
102  char TRANS('N');
103  long NRHS(mat.N), LDA(M), LDB(mat.M),
104  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
105  double *WORK(new double[LWORK]);
106  dgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
107  delete [] WORK;
108 
109  if(M>N){ //overdetermined
110  dgematrix tmp(N,mat.N);
111  for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
112  tmp(i,j) =mat(i,j);
113  }}
114  mat.clear();
115  swap(mat,tmp);
116  }
117 
118  if(INFO!=0){
119  std::cerr << "[WARNING] dgematrix::dgels(dgematrix&) "
120  << "Serious trouble happend. INFO = " << INFO << "."
121  << std::endl;
122  }
123  return INFO;
124 }
125 
126 //=============================================================================
129 inline long dgematrix::dgels(dcovector& vec)
130 {
131 #ifdef CPPL_VERBOSE
132  std::cerr << "# [MARK] dgematrix::dgels(dcovector&)"
133  << std::endl;
134 #endif//CPPL_VERBOSE
135 
136 #ifdef CPPL_DEBUG
137  if(M!=vec.L){
138  std::cerr << "[ERROR] dgematrix::dgels(dcovector&) " << std::endl
139  << "These matrix and vector cannot be solved." << std::endl
140  << "Your input was (" << M << "x" << N << ") and ("
141  << vec.L << ")." << std::endl;
142  exit(1);
143  }
144 #endif//CPPL_DEBUG
145 
146  if(M<N){ //underdetermined
147  dcovector tmp(N);
148  for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
149  vec.clear();
150  swap(vec,tmp);
151  }
152 
153  char TRANS('N');
154  long NRHS(1), LDA(M), LDB(vec.L),
155  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
156  double *WORK(new double[LWORK]);
157  dgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
158  delete [] WORK;
159 
160  if(M>N){ //overdetermined
161  dcovector tmp(N);
162  for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
163  vec.clear();
164  swap(vec,tmp);
165  }
166 
167  if(INFO!=0){
168  std::cerr << "[WARNING] dgematrix::dgels(dcovector&) "
169  << "Serious trouble happend. INFO = " << INFO << "."
170  << std::endl;
171  }
172  return INFO;
173 }
174 
175 //=============================================================================
182 inline long dgematrix::dgels(dgematrix& mat, drovector& residual)
183 {
184 #ifdef CPPL_VERBOSE
185  std::cerr << "# [MARK] dgematrix::dgels(dgematrix&, drovector&)"
186  << std::endl;
187 #endif//CPPL_VERBOSE
188 
189 #ifdef CPPL_DEBUG
190  if(M!=mat.M){
191  std::cerr << "[ERROR] dgematrix::dgels(dgematrix&, drovector&) "
192  << std::endl
193  << "These two matrices cannot be solved." << std::endl
194  << "Your input was (" << M << "x" << N << ") and ("
195  << mat.M << "x" << mat.N << ")." << std::endl;
196  exit(1);
197  }
198 #endif//CPPL_DEBUG
199 
200  residual.resize(mat.N); residual.zero();
201 
202  if(M<N){ //underdetermined
203  dgematrix tmp(N,mat.N);
204  for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
205  tmp(i,j) =mat(i,j);
206  }}
207  mat.clear();
208  swap(mat,tmp);
209  }
210 
211  char TRANS('N');
212  long NRHS(mat.N), LDA(M), LDB(mat.M),
213  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
214  double *WORK(new double[LWORK]);
215  dgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
216  delete [] WORK;
217 
218  if(M>N){ //overdetermined
219  for(long i=0; i<residual.L; i++){ for(long j=0; j<M-N; j++){
220  residual(i) += std::pow(mat(N+j,i), 2.0);
221  }}
222 
223  dgematrix tmp(N,mat.N);
224  for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
225  tmp(i,j) =mat(i,j);
226  }}
227  mat.clear();
228  swap(mat,tmp);
229  }
230 
231  if(INFO!=0){
232  std::cerr << "[WARNING] dgematrix::dgels(dgematrix&, drovector&) "
233  << "Serious trouble happend. INFO = " << INFO << "."
234  << std::endl;
235  }
236  return INFO;
237 }
238 
239 //=============================================================================
246 inline long dgematrix::dgels(dcovector& vec, double& residual)
247 {
248 #ifdef CPPL_VERBOSE
249  std::cerr << "# [MARK] dgematrix::dgels(dcovector&, double&)"
250  << std::endl;
251 #endif//CPPL_VERBOSE
252 
253 #ifdef CPPL_DEBUG
254  if(M!=vec.L){
255  std::cerr << "[ERROR] dgematrix::dgels(dcovector&, double&) " << std::endl
256  << "These matrix and vector cannot be solved." << std::endl
257  << "Your input was (" << M << "x" << N << ") and ("
258  << vec.L << ")." << std::endl;
259  exit(1);
260  }
261 #endif//CPPL_DEBUG
262 
263  residual=0.0;
264 
265  if(M<N){ //underdetermined
266  dcovector tmp(N);
267  for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
268  vec.clear();
269  swap(vec,tmp);
270  }
271 
272  char TRANS('N');
273  long NRHS(1), LDA(M), LDB(vec.L),
274  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
275  double *WORK(new double[LWORK]);
276  dgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
277  delete [] WORK;
278 
279  if(M>N){ //overdetermined
280  for(long i=0; i<M-N; i++){ residual+=std::pow(vec(N+i),2.0); }
281 
282  dcovector tmp(N);
283  for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
284  vec.clear();
285  swap(vec,tmp);
286  }
287 
288  if(INFO!=0){
289  std::cerr << "[WARNING] dgematrix::dgels(dcovector&, double&) "
290  << "Serious trouble happend. INFO = " << INFO << "."
291  << std::endl;
292  }
293  return INFO;
294 }
295 
299 
300 //=============================================================================
304 inline long dgematrix::dgelss(dcovector& B, dcovector& S, long& RANK,
305  const double RCOND =-1. )
306 {
307 #ifdef CPPL_VERBOSE
308  std::cerr << "# [MARK] dgematrix::dgelss(dcovector&, dcovector&, long&, const double)"
309  << std::endl;
310 #endif//CPPL_VERBOSE
311 
312 #ifdef CPPL_DEBUG
313  if(M!=B.L){
314  std::cerr << "[ERROR] dgematrix::dgelss"
315  << "(dcovector&, dcovector&, long, double) " << std::endl
316  << "These matrix and vector cannot be solved." << std::endl
317  << "Your input was (" << M << "x" << N << ") and ("
318  << B.L << ")." << std::endl;
319  exit(1);
320  }
321 #endif//CPPL_DEBUG
322 
323  if(M<N){ //underdetermined
324  dcovector tmp(N);
325  for(long i=0; i<B.L; i++){ tmp(i)=B(i); }
326  B.clear();
327  swap(B,tmp);
328  }
329 
330  S.resize(min(M,N));
331 
332  long NRHS(1), LDA(M), LDB(B.L),
333  LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
334  double *WORK(new double[LWORK]);
335  dgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.Array,
336  RCOND, RANK, WORK, LWORK, INFO);
337  delete [] WORK;
338 
339  if(INFO!=0){
340  std::cerr << "[WARNING] dgematrix::dgelss"
341  << "(dcovector&, docvector&, long, const double) "
342  << "Serious trouble happend. INFO = " << INFO << "."
343  << std::endl;
344  }
345  return INFO;
346 }
347 
348 //=============================================================================
352 inline long dgematrix::dgelss(dgematrix& B, dcovector& S, long& RANK,
353  const double RCOND =-1. )
354 {
355 #ifdef CPPL_VERBOSE
356  std::cerr << "# [MARK] dgematrix::dgelss(dgematrix&, dcovector&, long& const double)"
357  << std::endl;
358 #endif//CPPL_VERBOSE
359 
360 #ifdef CPPL_DEBUG
361  if(M!=B.M){
362  std::cerr << "[ERROR] dgematrix::dgelss"
363  << "(dgematrix&, dcovector&, long&, const double) " << std::endl
364  << "These matrix and vector cannot be solved." << std::endl
365  << "Your input was (" << M << "x" << N << ") and ("
366  << B.M << "x" << B.N << ")." << std::endl;
367  exit(1);
368  }
369 #endif//CPPL_DEBUG
370 
371  if(M<N){ //underdetermined
372  dgematrix tmp(N,B.N);
373  for(long i=0; i<B.M; i++){
374  for(long j=0; j<B.N; j++){
375  tmp(i,j)=B(i,j);
376  }
377  }
378  B.clear();
379  swap(B,tmp);
380  }
381 
382  S.resize(min(M,N));
383 
384  long NRHS(B.N), LDA(M), LDB(B.M),
385  LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
386  double *WORK(new double[LWORK]);
387  dgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.Array,
388  RCOND, RANK, WORK, LWORK, INFO);
389  delete [] WORK;
390 
391  if(INFO!=0){
392  std::cerr << "[WARNING] dgematrix::dgelss"
393  << "(dgematrix&, docvector&, long, const double) "
394  << "Serious trouble happend. INFO = " << INFO << "."
395  << std::endl;
396  }
397  return INFO;
398 }
399 
403 //=============================================================================
410 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi)
411 {
412 #ifdef CPPL_VERBOSE
413  std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&)"
414  << std::endl;
415 #endif//CPPL_VERBOSE
416 
417 #ifdef CPPL_DEBUG
418  if(M!=N){
419  std::cerr << "[ERROR] dgematrix::dgeev"
420  << "(vector<double>&, vector<double>&) " << std::endl
421  << "This matrix is not a square matrix." << std::endl
422  << "This matrix is (" << M << "x" << N << ")." << std::endl;
423  exit(1);
424  }
425 #endif//CPPL_DEBUG
426 
427  wr.resize(N); wi.resize(N);
428  char JOBVL('N'), JOBVR('N');
429  long LDA(N), LDVL(1), LDVR(1), LWORK(3*N), INFO(1);
430  double *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
431  dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
432  VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
433  delete [] WORK; delete [] VL; delete [] VL;
434 
435  if(INFO!=0){
436  std::cerr << "[WARNING] dgematrix::dgeev"
437  << "(vector<double>&, vector<double>&)"
438  << std::endl
439  << "Serious trouble happend. INFO = " << INFO << "."
440  << std::endl;
441  }
442  return INFO;
443 }
444 
445 //=============================================================================
453 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
454  std::vector<dcovector>& vrr,
455  std::vector<dcovector>& vri)
456 {
457 #ifdef CPPL_VERBOSE
458  std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
459  << std::endl;
460 #endif//CPPL_VERBOSE
461 
462 #ifdef CPPL_DEBUG
463  if(M!=N){
464  std::cerr << "[ERROR] dgematrix::dgeev"
465  << "(vector<double>&, vector<double>&, "
466  << "vector<dcovector>&, vector<dcovector>&) "
467  << std::endl
468  << "This matrix is not a square matrix." << std::endl
469  << "This matrix is (" << M << "x" << N << ")." << std::endl;
470  exit(1);
471  }
472 #endif//CPPL_DEBUG
473 
474  wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
475  for(long i=0; i<N; i++){ vrr[i].resize(N); vri[i].resize(N); }
476  dgematrix VR(N,N);
477  char JOBVL('N'), JOBVR('V');
478  long LDA(N), LDVL(1), LDVR(N), LWORK(4*N), INFO(1);
479  double *VL(NULL), *WORK(new double[LWORK]);
480  dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
481  VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
482  delete [] WORK; delete [] VL;
483 
485  for(long j=0; j<N; j++){
486  if(fabs(wi[j])<1e-10){
487  for(long i=0; i<N; i++){
488  vrr[j](i) = VR(i,j); vri[j](i) = 0.0;
489  }
490  }
491  else{
492  for(long i=0; i<N; i++){
493  vrr[j](i) = VR(i,j); vri[j](i) = VR(i,j+1);
494  vrr[j+1](i) = VR(i,j); vri[j+1](i) =-VR(i,j+1);
495  }
496  j++;
497  }
498  }
499 
500  if(INFO!=0){
501  std::cerr << "[WARNING] dgematrix::dgeev"
502  << "(vector<double>&, vector<double>&, "
503  << "vector<dcovector>&, vector<dcovector>&) "
504  << std::endl
505  << "Serious trouble happend. INFO = " << INFO << "."
506  << std::endl;
507  }
508  return INFO;
509 }
510 
511 //=============================================================================
519 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
520  std::vector<drovector>& vlr,
521  std::vector<drovector>& vli)
522 {
523 #ifdef CPPL_VERBOSE
524  std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
525  << std::endl;
526 #endif//CPPL_VERBOSE
527 
528 #ifdef CPPL_DEBUG
529  if(M!=N){
530  std::cerr << "[ERROR] dgematrix::dgeev"
531  << "(vector<double>&, vector<double>&, "
532  << "vector<drovector>&, vector<drovector>&) "
533  << std::endl
534  << "This matrix is not a square matrix." << std::endl
535  << "This matrix is (" << M << "x" << N << ")." << std::endl;
536  exit(1);
537  }
538 #endif//CPPL_DEBUG
539 
540  wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
541  for(long i=0; i<N; i++){ vlr[i].resize(N); vli[i].resize(N); }
542  dgematrix VL(N,N);
543  char JOBVL('V'), JOBVR('N');
544  long LDA(N), LDVL(N), LDVR(1), LWORK(4*N), INFO(1);
545  double *VR(NULL), *WORK(new double[LWORK]);
546  dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
547  VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
548  delete [] WORK; delete [] VR;
549 
551  for(long j=0; j<N; j++){
552  if(fabs(wi[j])<1e-10){
553  for(long i=0; i<N; i++){
554  vlr[j](i) = VL(i,j); vli[j](i) = 0.0;
555  }
556  }
557  else{
558  for(long i=0; i<N; i++){
559  vlr[j](i) = VL(i,j); vli[j](i) =-VL(i,j+1);
560  vlr[j+1](i) = VL(i,j); vli[j+1](i) = VL(i,j+1);
561  }
562  j++;
563  }
564  }
565 
566 
567  if(INFO!=0){
568  std::cerr << "[WARNING] dgematrix::dgeev"
569  << "(vector<double>&, vector<double>&, "
570  << "vector<drovector>&, vector<drovector>&) "
571  << std::endl
572  << "Serious trouble happend. INFO = " << INFO << "."
573  << std::endl;
574  }
575  return INFO;
576 }
577 
578 
579 
583 
584 //=============================================================================
591 inline long dgematrix::dggev(dgematrix& matB,
592  std::vector<double>& wr, std::vector<double>& wi)
593 {
594 #ifdef CPPL_VERBOSE
595  std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&)"
596  << std::endl;
597 #endif//CPPL_VERBOSE
598 
599 #ifdef CPPL_DEBUG
600  if(M!=N){
601  std::cerr << "[ERROR] dgematrix::dggev"
602  << "(dgematrix&, vector<double>&, vector<double>&) " << std::endl
603  << "This matrix is not a square matrix." << std::endl
604  << "This matrix is (" << M << "x" << N << ")." << std::endl;
605  exit(1);
606  }
607  if(matB.M!=N || matB.N!=N){
608  std::cerr << "[ERROR] dgematrix::dggev"
609  << "(dgematrix&, vector<double>&, vector<double>&) " << std::endl
610  << "The matrix B is not a square matrix "
611  << "having the same size as \"this\" matrix." << std::endl
612  << "The B matrix is (" << matB.M << "x" << matB.N << ")."
613  << std::endl;
614  exit(1);
615  }
616 #endif//CPPL_DEBUG
617 
618  wr.resize(N); wi.resize(N);
619  char JOBVL('N'), JOBVR('N');
620  long LDA(N), LDB(N), LDVL(1), LDVR(1), LWORK(8*N), INFO(1);
621  double *BETA(new double[N]), *VL(NULL), *VR(NULL),
622  *WORK(new double[LWORK]);
623  dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
624  BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
625  delete [] WORK; delete [] VL; delete [] VR;
626 
628  for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
629  delete [] BETA;
630 
631  if(INFO!=0){
632  std::cerr << "[WARNING] dgematrix::dgeev"
633  << "(dgematrix&, vector<double>&, vector<double>&)"
634  << std::endl
635  << "Serious trouble happend. INFO = " << INFO << "."
636  << std::endl;
637  }
638  return INFO;
639 }
640 
641 //=============================================================================
649 inline long dgematrix::dggev(dgematrix& matB,
650  std::vector<double>& wr, std::vector<double>& wi,
651  std::vector<dcovector>& vrr,
652  std::vector<dcovector>& vri)
653 {
654 #ifdef CPPL_VERBOSE
655  std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
656  << std::endl;
657 #endif//CPPL_VERBOSE
658 
659 #ifdef CPPL_DEBUG
660  if(M!=N){
661  std::cerr << "[ERROR] dgematrix::dggev"
662  << "(dgematrix&, vector<double>&, vector<double>&, "
663  << "vector<dcovector>&, vector<dcovector>&)" << std::endl
664  << "This matrix is not a square matrix." << std::endl
665  << "This matrix is (" << M << "x" << N << ")." << std::endl;
666  exit(1);
667  }
668  if(matB.M!=N || matB.N!=N){
669  std::cerr << "[ERROR] dgematrix::dggev"
670  << "(dgematrix&, vector<double>&, vector<double>&, "
671  << "vector<dcovector>&, vector<dcovector>&)" << std::endl
672  << "The matrix B is not a square matrix "
673  << "having the same size as \"this\" matrix." << std::endl
674  << "The B matrix is (" << matB.M << "x" << matB.N << ")."
675  << std::endl;
676  exit(1);
677  }
678 #endif//CPPL_DEBUG
679 
680  wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
681  for(long i=0; i<N; i++){ vrr[i].resize(N); vri[i].resize(N); }
682  dgematrix VR(N,N);
683  char JOBVL('N'), JOBVR('V');
684  long LDA(N), LDB(N), LDVL(1), LDVR(N), LWORK(8*N), INFO(1);
685  double *BETA(new double[N]), *VL(NULL), *WORK(new double[LWORK]);
686  dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
687  BETA, VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
688  delete [] WORK; delete [] VL;
689 
691  for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
692  delete [] BETA;
693 
695  for(long j=0; j<N; j++){
696  if(fabs(wi[j])<1e-10){
697  for(long i=0; i<N; i++){
698  vrr[j](i) = VR(i,j); vri[j](i) = 0.0;
699  }
700  }
701  else{
702  for(long i=0; i<N; i++){
703  vrr[j](i) = VR(i,j); vri[j](i) = VR(i,j+1);
704  vrr[j+1](i) = VR(i,j); vri[j+1](i) =-VR(i,j+1);
705  }
706  j++;
707  }
708  }
709 
710  if(INFO!=0){
711  std::cerr << "[WARNING] dgematrix::dggev"
712  << "(dgematrix&, vector<double>&, vector<double>&, "
713  << "vector<dcovector>&, vector<dcovector>&)" << std::endl
714  << "Serious trouble happend. INFO = " << INFO << "."
715  << std::endl;
716  }
717  return INFO;
718 }
719 
720 //=============================================================================
728 inline long dgematrix::dggev(dgematrix& matB,
729  std::vector<double>& wr, std::vector<double>& wi,
730  std::vector<drovector>& vlr,
731  std::vector<drovector>& vli)
732 {
733 #ifdef CPPL_VERBOSE
734  std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
735  << std::endl;
736 #endif//CPPL_VERBOSE
737 
738 #ifdef CPPL_DEBUG
739  if(M!=N){
740  std::cerr << "[ERROR] dgematrix::dggev"
741  << "(dgematrix&, vector<double>&, vector<double>&, "
742  << "vector<drovector>&, vector<drovector>&)" << std::endl
743  << "This matrix is not a square matrix." << std::endl
744  << "This matrix is (" << M << "x" << N << ")." << std::endl;
745  exit(1);
746  }
747  if(matB.M!=N || matB.N!=N){
748  std::cerr << "[ERROR] dgematrix::dggev"
749  << "(dgematrix&, vector<double>&, vector<double>&, "
750  << "vector<drovector>&, vector<drovector>&)" << std::endl
751  << "The matrix B is not a square matrix "
752  << "having the same size as \"this\" matrix." << std::endl
753  << "The B matrix is (" << matB.M << "x" << matB.N << ")."
754  << std::endl;
755  exit(1);
756  }
757 #endif//CPPL_DEBUG
758 
759  wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
760  for(long i=0; i<N; i++){ vlr[i].resize(N); vli[i].resize(N); }
761  dgematrix VL(N,N);
762  char JOBVL('V'), JOBVR('N');
763  long LDA(N), LDB(N), LDVL(N), LDVR(1), LWORK(8*N), INFO(1);
764  double *BETA(new double[N]), *VR(NULL), *WORK(new double[LWORK]);
765  dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
766  BETA, VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
767  delete [] WORK; delete [] VR;
768 
770  for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
771  delete [] BETA;
772 
774  for(long j=0; j<N; j++){
775  if(fabs(wi[j])<1e-10){
776  for(long i=0; i<N; i++){
777  vlr[j](i) = VL(i,j); vli[j](i) = 0.0;
778  }
779  }
780  else{
781  for(long i=0; i<N; i++){
782  vlr[j](i) = VL(i,j); vli[j](i) =-VL(i,j+1);
783  vlr[j+1](i) = VL(i,j); vli[j+1](i) = VL(i,j+1);
784  }
785  j++;
786  }
787  }
788 
789  if(INFO!=0){
790  std::cerr << "[WARNING] dgematrix::dggev"
791  << "(dgematrix&, vector<double>&, vector<double>&, "
792  << "vector<drovector>&, vector<drovector>&)" << std::endl
793  << "Serious trouble happend. INFO = " << INFO << "."
794  << std::endl;
795  }
796  return INFO;
797 }
798 
802 
803 //=============================================================================
813 {
814 #ifdef CPPL_VERBOSE
815  std::cerr << "# [MARK] dgematrix::dgesvd(dcovector&, dgematrix&, dgematrix&)"
816  << std::endl;
817 #endif//CPPL_VERBOSE
818 
819  char JOBU('A'), JOBVT('A');
820  long LDA(M), LDU(M), LDVT(N),
821  LWORK(max(3*min(M,N)+max(M,N),5*min(M,N))), INFO(1);
822  double *WORK(new double[LWORK]);
823  S.resize(min(M,N)); U.resize(LDU,M); VT.resize(LDVT,N);
824 
825  dgesvd_(JOBU, JOBVT, M, N, Array, LDA, S.Array, U.Array,
826  LDU, VT.Array, LDVT, WORK, LWORK, INFO);
827  delete [] WORK;
828 
829  if(INFO!=0){
830  std::cerr << "[WARNING] dgematrix::dgesvd"
831  << "(dceovector&, dgematrix&, dcovector&) "
832  << "Serious trouble happend. INFO = " << INFO << "."
833  << std::endl;
834  }
835  return INFO;
836 }
void clear()
Definition: dgematrix-misc.hpp:3
long dgelss(dcovector &, dcovector &, long &, const double)
Definition: dgematrix-lapack.hpp:304
long dgels(dgematrix &)
Definition: dgematrix-lapack.hpp:76
void resize(const long &)
Definition: dcovector-misc.hpp:93
double max(double v1, double v2)
Return maximum.
Definition: variousFunctions.cpp:355
void resize(const long &, const long &)
Definition: dgematrix-misc.hpp:126
friend void swap(dgematrix &, dgematrix &)
Definition: dgematrix-misc.hpp:154
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
void zero()
Definition: drovector-misc.hpp:23
void clear()
Definition: dcovector-misc.hpp:3
Real Double-precision Row Vector Class.
Definition: drovector.hpp:3
long dgesv(dgematrix &)
Definition: dgematrix-lapack.hpp:6
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
long dgeev(std::vector< double > &, std::vector< double > &)
Definition: dgematrix-lapack.hpp:410
friend _dgematrix i(const dgematrix &)
Definition: dgematrix-calc.hpp:21
double * Array
1D Array to store vector data
Definition: _drovector.hpp:8
void resize(const long &)
Definition: drovector-misc.hpp:93
long dgesvd(dcovector &, dgematrix &, dgematrix &)
Definition: dgematrix-lapack.hpp:812
long dggev(dgematrix &, std::vector< double > &, std::vector< double > &)
Definition: dgematrix-lapack.hpp:591