CPPLapack
 All Classes Files Functions Variables Friends Pages
dgematrix-lapack.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*X=Y using dgesv\n
3  The argument is dgematrix Y. Y is overwritten and become the solution X.
4  A is also overwritten and become P*l*U.
5 */
6 inline CPPL_INT dgematrix::dgesv(dgematrix& mat)
7 {CPPL_VERBOSE_REPORT;
8 #ifdef CPPL_DEBUG
9  if(m!=n || m!=mat.m){
10  ERROR_REPORT;
11  std::cerr << "These two matrices cannot be solved." << std::endl
12  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
13  exit(1);
14  }
15 #endif//CPPL_DEBUG
16 
17  CPPL_INT NRHS(mat.n), LDA(n), *IPIV(new CPPL_INT[n]), LDB(mat.m), INFO(1);
18  dgesv_(&n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, &INFO);
19  delete [] IPIV;
20 
21  if(INFO!=0){
22  WARNING_REPORT;
23  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
24  }
25  return INFO;
26 }
27 
28 //=============================================================================
29 /*! solve A*x=y using dgesv\n
30  The argument is dcovector y. y is overwritten and become the solution x.
31  A is also overwritten and become P*l*U.
32 */
33 inline CPPL_INT dgematrix::dgesv(dcovector& vec)
34 {CPPL_VERBOSE_REPORT;
35 #ifdef CPPL_DEBUG
36  if(m!=n || m!=vec.l){
37  ERROR_REPORT;
38  std::cerr << "These matrix and vector cannot be solved." << std::endl
39  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
40  exit(1);
41  }
42 #endif//CPPL_DEBUG
43 
44  CPPL_INT NRHS(1), LDA(n), *IPIV(new CPPL_INT[n]), LDB(vec.l), INFO(1);
45  dgesv_(&n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, &INFO);
46  delete [] IPIV;
47 
48  if(INFO!=0){
49  WARNING_REPORT;
50  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
51  }
52  return INFO;
53 }
54 
55 ///////////////////////////////////////////////////////////////////////////////
56 ///////////////////////////////////////////////////////////////////////////////
57 ///////////////////////////////////////////////////////////////////////////////
58 
59 //=============================================================================
60 /*! solve overdetermined or underdetermined A*X=Y using dgels\n
61  */
62 inline CPPL_INT dgematrix::dgels(dgematrix& mat)
63 {CPPL_VERBOSE_REPORT;
64 #ifdef CPPL_DEBUG
65  if(m!=mat.m){
66  ERROR_REPORT;
67  std::cerr << "These two matrices cannot be solved." << std::endl
68  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
69  exit(1);
70  }
71 #endif//CPPL_DEBUG
72 
73  if(m<n){ //underdetermined
74  dgematrix tmp(n,mat.n);
75  for(CPPL_INT i=0; i<mat.m; i++){
76  for(CPPL_INT j=0; j<mat.n; j++){
77  tmp(i,j) =mat(i,j);
78  }
79  }
80  mat.clear();
81  swap(mat,tmp);
82  }
83 
84  char TRANS('n');
85  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
86  double *WORK(new double[LWORK]);
87  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
88  delete [] WORK;
89 
90  if(m>n){ //overdetermined
91  dgematrix tmp(n,mat.n);
92  for(CPPL_INT i=0; i<tmp.m; i++){
93  for(CPPL_INT j=0; j<tmp.n; j++){
94  tmp(i,j) =mat(i,j);
95  }
96  }
97  mat.clear();
98  swap(mat,tmp);
99  }
100 
101  if(INFO!=0){
102  WARNING_REPORT;
103  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
104  }
105  return INFO;
106 }
107 
108 //=============================================================================
109 /*! solve overdetermined or underdetermined A*x=y using dgels\n
110  */
111 inline CPPL_INT dgematrix::dgels(dcovector& vec)
112 {CPPL_VERBOSE_REPORT;
113 #ifdef CPPL_DEBUG
114  if(m!=vec.l){
115  ERROR_REPORT;
116  std::cerr << "These matrix and vector cannot be solved." << std::endl
117  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
118  exit(1);
119  }
120 #endif//CPPL_DEBUG
121 
122  if(m<n){ //underdetermined
123  dcovector tmp(n);
124  for(CPPL_INT i=0; i<vec.l; i++){
125  tmp(i)=vec(i);
126  }
127  vec.clear();
128  swap(vec,tmp);
129  }
130 
131  char TRANS('n');
132  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
133  double *WORK(new double[LWORK]);
134  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
135  delete [] WORK;
136 
137  if(m>n){ //overdetermined
138  dcovector tmp(n);
139  for(CPPL_INT i=0; i<tmp.l; i++){
140  tmp(i)=vec(i);
141  }
142  vec.clear();
143  swap(vec,tmp);
144  }
145 
146  if(INFO!=0){
147  WARNING_REPORT;
148  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
149  }
150  return INFO;
151 }
152 
153 //=============================================================================
154 /*! solve overdetermined or underdetermined A*X=Y using dgels
155  with the sum of residual squares output\n
156  The residual is set as the columnwise sum of residual squares
157  for overdetermined problems
158  while it is always zero for underdetermined problems.
159 */
160 inline CPPL_INT dgematrix::dgels(dgematrix& mat, drovector& residual)
161 {CPPL_VERBOSE_REPORT;
162 #ifdef CPPL_DEBUG
163  if(m!=mat.m){
164  ERROR_REPORT;
165  std::cerr << "These two matrices cannot be solved." << std::endl
166  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
167  exit(1);
168  }
169 #endif//CPPL_DEBUG
170 
171  residual.resize(mat.n);
172  residual.zero();
173 
174  if(m<n){ //underdetermined
175  dgematrix tmp(n,mat.n);
176  for(CPPL_INT i=0; i<mat.m; i++){
177  for(CPPL_INT j=0; j<mat.n; j++){
178  tmp(i,j) =mat(i,j);
179  }
180  }
181  mat.clear();
182  swap(mat,tmp);
183  }
184 
185  char TRANS('n');
186  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
187  double *WORK(new double[LWORK]);
188  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
189  delete [] WORK;
190 
191  if(m>n){ //overdetermined
192  for(CPPL_INT i=0; i<residual.l; i++){
193  for(CPPL_INT j=0; j<m-n; j++){
194  residual(i) += std::pow(mat(n+j,i), 2.0);
195  }
196  }
197 
198  dgematrix tmp(n,mat.n);
199  for(CPPL_INT i=0; i<tmp.m; i++){
200  for(CPPL_INT j=0; j<tmp.n; j++){
201  tmp(i,j) =mat(i,j);
202  }
203  }
204  mat.clear();
205  swap(mat,tmp);
206  }
207 
208  if(INFO!=0){
209  WARNING_REPORT;
210  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
211  }
212  return INFO;
213 }
214 
215 //=============================================================================
216 /*! solve overdetermined or underdetermined A*x=y using dgels
217  with the sum of residual squares output\n
218  The residual is set as the sum of residual squares
219  for overdetermined problems
220  while it is always zero for underdetermined problems.
221 */
222 inline CPPL_INT dgematrix::dgels(dcovector& vec, double& residual)
223 {CPPL_VERBOSE_REPORT;
224 #ifdef CPPL_DEBUG
225  if(m!=vec.l){
226  ERROR_REPORT;
227  std::cerr << "These matrix and vector cannot be solved." << std::endl
228  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
229  exit(1);
230  }
231 #endif//CPPL_DEBUG
232 
233  residual=0.0;
234 
235  if(m<n){ //underdetermined
236  dcovector tmp(n);
237  for(CPPL_INT i=0; i<vec.l; i++){
238  tmp(i)=vec(i);
239  }
240  vec.clear();
241  swap(vec,tmp);
242  }
243 
244  char TRANS('n');
245  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
246  double *WORK(new double[LWORK]);
247  dgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
248  delete [] WORK;
249 
250  if(m>n){ //overdetermined
251  for(CPPL_INT i=0; i<m-n; i++){
252  residual+=std::pow(vec(n+i),2.0);
253  }
254 
255  dcovector tmp(n);
256  for(CPPL_INT i=0; i<tmp.l; i++){
257  tmp(i)=vec(i);
258  }
259  vec.clear();
260  swap(vec,tmp);
261  }
262 
263  if(INFO!=0){
264  WARNING_REPORT;
265  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
266  }
267  return INFO;
268 }
269 
270 ///////////////////////////////////////////////////////////////////////////////
271 ///////////////////////////////////////////////////////////////////////////////
272 ///////////////////////////////////////////////////////////////////////////////
273 
274 //=============================================================================
275 /*! calculate the least-squares-least-norm solution for overdetermined or
276  underdetermined A*x=y using dgelss\n
277 */
278 inline CPPL_INT dgematrix::dgelss
279 (
280  dcovector& B,
281  dcovector& S,
282  CPPL_INT& RANK,
283  const double RCOND =-1.
284  )
285 {CPPL_VERBOSE_REPORT;
286 #ifdef CPPL_DEBUG
287  if(m!=B.l){
288  ERROR_REPORT;
289  std::cerr << "These matrix and vector cannot be solved." << std::endl
290  << "Your input was (" << m << "x" << n << ") and (" << B.l << ")." << std::endl;
291  exit(1);
292  }
293 #endif//CPPL_DEBUG
294 
295  if(m<n){ //underdetermined
296  dcovector tmp(n);
297  for(CPPL_INT i=0; i<B.l; i++){
298  tmp(i)=B(i);
299  }
300  B.clear();
301  swap(B,tmp);
302  }
303 
304  S.resize(std::min(m,n));
305 
306  CPPL_INT NRHS(1), LDA(m), LDB(B.l), LWORK(3*std::min(m,n)+std::max(std::max(2*std::min(m,n),std::max(m,n)), NRHS)), INFO(1);
307  double *WORK(new double[LWORK]);
308  dgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, &INFO);
309  delete [] WORK;
310 
311  if(INFO!=0){
312  WARNING_REPORT;
313  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
314  }
315  return INFO;
316 }
317 
318 //=============================================================================
319 /*! calculate the least-squares-least-norm solution for overdetermined or
320  underdetermined A*x=y using dgelss\n
321 */
322 inline CPPL_INT dgematrix::dgelss
323 (
324  dgematrix& B,
325  dcovector& S,
326  CPPL_INT& RANK,
327  const double RCOND =-1.
328  )
329 {CPPL_VERBOSE_REPORT;
330 #ifdef CPPL_DEBUG
331  if(m!=B.m){
332  ERROR_REPORT;
333  std::cerr << "These matrix and vector cannot be solved." << std::endl
334  << "Your input was (" << m << "x" << n << ") and (" << B.m << "x" << B.n << ")." << std::endl;
335  exit(1);
336  }
337 #endif//CPPL_DEBUG
338 
339  if(m<n){ //underdetermined
340  dgematrix tmp(n,B.n);
341  for(CPPL_INT i=0; i<B.m; i++){
342  for(CPPL_INT j=0; j<B.n; j++){
343  tmp(i,j)=B(i,j);
344  }
345  }
346  B.clear();
347  swap(B,tmp);
348  }
349 
350  S.resize(std::min(m,n));
351 
352  CPPL_INT NRHS(B.n), LDA(m), LDB(B.m), LWORK(3*std::min(m,n)+std::max(std::max(2*std::min(m,n),std::max(m,n)), NRHS)), INFO(1);
353  double *WORK(new double[LWORK]);
354  dgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, &INFO);
355  delete [] WORK;
356 
357  if(INFO!=0){
358  WARNING_REPORT;
359  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
360  }
361  return INFO;
362 }
363 
364 ///////////////////////////////////////////////////////////////////////////////
365 ///////////////////////////////////////////////////////////////////////////////
366 ///////////////////////////////////////////////////////////////////////////////
367 
368 //=============================================================================
369 /*! calculate the least-squares-least-norm solution for overdetermined or underdetermined A*x=y using dgelsd\n
370 */
371 inline CPPL_INT dgematrix::dgelsd
372 (
373  dcovector& B,
374  dcovector& S,
375  CPPL_INT& RANK,
376  const double RCOND =-1.
377  )
378 {CPPL_VERBOSE_REPORT;
379 #ifdef CPPL_DEBUG
380  if(m!=B.l){
381  ERROR_REPORT;
382  std::cerr << "These matrix and vector cannot be solved." << std::endl
383  << "Your input was (" << m << "x" << n << ") and (" << B.l << ")." << std::endl;
384  exit(1);
385  }
386 #endif//CPPL_DEBUG
387 
388  //////// resize ////////
389  if(m<n){ //underdetermined
390  dcovector tmp(n);
391  for(CPPL_INT i=0; i<B.l; i++){
392  tmp(i)=B(i);
393  }
394  B.clear();
395  swap(B,tmp);
396  }
397 
398  S.resize(std::min(m,n));
399 
400  //////// prerun ////////
401  CPPL_INT M =m;
402  CPPL_INT N =n;
403  CPPL_INT NRHS =1;
404  CPPL_INT LDA =m;
405  CPPL_INT LDB =B.l;
406  std::vector<double> WORK(1);
407  CPPL_INT LWORK =-1;
408  CPPL_INT MINMN =std::min(M,N);
409  CPPL_INT SMLSIZ =25;
410  CPPL_INT NLVL =CPPL_INT( std::log(double(MINMN)/double(SMLSIZ+1))/std::log(2.) );
411  CPPL_INT LIWORK =3*MINMN*NLVL+11*MINMN;
412  std::vector<CPPL_INT> IWORK(LIWORK);
413  CPPL_INT INFO =1;
414  dgelsd_(&M, &N, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, &WORK[0], &LWORK, &IWORK[0], &INFO);
415  LWORK =CPPL_INT(WORK[0]);
416  WORK.resize(LWORK);
417 
418  //////// run ////////
419  INFO =1;
420  RANK =1;
421  dgelsd_(&M, &N, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, &WORK[0], &LWORK, &IWORK[0], &INFO);
422 
423  //////// info ////////
424  if(INFO!=0){
425  WARNING_REPORT;
426  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
427  }
428  return INFO;
429 }
430 
431 ///////////////////////////////////////////////////////////////////////////////
432 ///////////////////////////////////////////////////////////////////////////////
433 ///////////////////////////////////////////////////////////////////////////////
434 
435 //=============================================================================
436 /*! calculate eigenvalues\n
437  All of the arguments need not to be initialized.
438  wr and wi are overwitten and become
439  real and imaginary part of eigenvalues, respectively.
440  This matrix is also overwritten.
441 */
442 inline CPPL_INT dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi)
443 {CPPL_VERBOSE_REPORT;
444 #ifdef CPPL_DEBUG
445  if(m!=n){
446  ERROR_REPORT;
447  std::cerr << "This matrix is not a square matrix." << std::endl
448  << "This matrix is (" << m << "x" << n << ")." << std::endl;
449  exit(1);
450  }
451 #endif//CPPL_DEBUG
452 
453  wr.resize(n);
454  wi.resize(n);
455  char JOBVL('n'), JOBVR('n');
456  CPPL_INT LDA(n), LDVL(1), LDVR(1), LWORK(3*n), INFO(1);
457  double *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
458  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
459  delete [] WORK;
460  delete [] VL;
461  delete [] VR;
462 
463  if(INFO!=0){
464  WARNING_REPORT;
465  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
466  }
467  return INFO;
468 }
469 
470 //=============================================================================
471 /*! calculate eigenvalues\n
472  All of the arguments need not to be initialized.
473  w are overwitten and become eigenvalues, respectively.
474  This matrix is also overwritten.
475 */
476 inline CPPL_INT dgematrix::dgeev(zcovector& w)
477 {CPPL_VERBOSE_REPORT;
478  //// call dgeev ////
479  std::vector<double> wr, wi;
480  CPPL_INT INFO =dgeev(wr,wi);
481 
482  //// assign ////
483  w.resize(n);
484  for(CPPL_INT i=0; i<n; i++){
485  w(i) =comple(wr[i],wi[i]);
486  }
487 
488  return INFO;
489 }
490 
491 //=============================================================================
492 /*! calculate right eigenvalues and right eigenvectors\n
493  All of the arguments need not to be initialized.
494  wr, wi, vrr, vri are overwitten and become
495  real and imaginary part of right eigenvalues and right eigenvectors,
496  respectively.
497  This matrix is also overwritten.
498 */
499 inline CPPL_INT dgematrix::dgeev
500 (
501  std::vector<double>& wr,
502  std::vector<double>& wi,
503  std::vector<dcovector>& vrr,
504  std::vector<dcovector>& vri
505  )
506 {CPPL_VERBOSE_REPORT;
507 #ifdef CPPL_DEBUG
508  if(m!=n){
509  ERROR_REPORT;
510  std::cerr << "This matrix is not a square matrix." << std::endl
511  << "This matrix is (" << m << "x" << n << ")." << std::endl;
512  exit(1);
513  }
514 #endif//CPPL_DEBUG
515 
516  wr.resize(n);
517  wi.resize(n);
518  vrr.resize(n);
519  vri.resize(n);
520  for(CPPL_INT i=0; i<n; i++){
521  vrr[i].resize(n);
522  vri[i].resize(n);
523  }
524 
525  dgematrix VR(n,n);
526  char JOBVL('n'), JOBVR('V');
527  CPPL_INT LDA(n), LDVL(1), LDVR(n), LWORK(4*n), INFO(1);
528  double *VL(NULL), *WORK(new double[LWORK]);
529  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL, &LDVL, VR.array, &LDVR, WORK, &LWORK, &INFO);
530  delete [] WORK;
531  delete [] VL;
532 
533  //// forming ////
534  for(CPPL_INT j=0; j<n; j++){
535  if(fabs(wi[j])<DBL_MIN){
536  for(CPPL_INT i=0; i<n; i++){
537  vrr[j](i) = VR(i,j);
538  vri[j](i) = 0.0;
539  }
540  }
541  else{
542  for(CPPL_INT i=0; i<n; i++){
543  vrr[j](i) = VR(i,j);
544  vri[j](i) = VR(i,j+1);
545  vrr[j+1](i) = VR(i,j);
546  vri[j+1](i) =-VR(i,j+1);
547  }
548  j++;
549  }
550  }
551 
552  if(INFO!=0){
553  WARNING_REPORT;
554  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
555  }
556  return INFO;
557 }
558 
559 //=============================================================================
560 /*! calculate left eigenvalues and left eigenvectors\n
561  All of the arguments need not to be initialized.
562  wr, wi, vrr, vri are overwitten and become
563  real and imaginary part of left eigenvalues and left eigenvectors,
564  respectively.
565  This matrix is also overwritten.
566 */
567 inline CPPL_INT dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
568  std::vector<drovector>& vlr,
569  std::vector<drovector>& vli)
570 {CPPL_VERBOSE_REPORT;
571 #ifdef CPPL_DEBUG
572  if(m!=n){
573  ERROR_REPORT;
574  std::cerr << "This matrix is not a square matrix." << std::endl
575  << "This matrix is (" << m << "x" << n << ")." << std::endl;
576  exit(1);
577  }
578 #endif//CPPL_DEBUG
579 
580  wr.resize(n);
581  wi.resize(n);
582  vlr.resize(n);
583  vli.resize(n);
584  for(CPPL_INT i=0; i<n; i++){
585  vlr[i].resize(n);
586  vli[i].resize(n);
587  }
588 
589  dgematrix VL(n,n);
590  char JOBVL('V'), JOBVR('n');
591  CPPL_INT LDA(n), LDVL(n), LDVR(1), LWORK(4*n), INFO(1);
592  double *VR(NULL), *WORK(new double[LWORK]);
593  dgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &wr[0], &wi[0], VL.array, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
594  delete [] WORK;
595  delete [] VR;
596 
597  //// forming ////
598  for(CPPL_INT j=0; j<n; j++){
599  if(fabs(wi[j])<DBL_MIN){
600  for(CPPL_INT i=0; i<n; i++){
601  vlr[j](i) = VL(i,j);
602  vli[j](i) = 0.0;
603  }
604  }
605  else{
606  for(CPPL_INT i=0; i<n; i++){
607  vlr[j](i) = VL(i,j);
608  vli[j](i) =-VL(i,j+1);
609  vlr[j+1](i) = VL(i,j);
610  vli[j+1](i) = VL(i,j+1);
611  }
612  j++;
613  }
614  }
615 
616 
617  if(INFO!=0){
618  WARNING_REPORT;
619  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
620  }
621  return INFO;
622 }
623 
624 
625 
626 ///////////////////////////////////////////////////////////////////////////////
627 ///////////////////////////////////////////////////////////////////////////////
628 ///////////////////////////////////////////////////////////////////////////////
629 
630 //=============================================================================
631 /*! calculate generalized eigenvalues\n
632  All of the arguments don't need to be initialized.
633  wr and wi are overwitten and become
634  real and imaginary part of generalized eigenvalues, respectively.
635  This matrix and matB are also overwritten.
636 */
637 inline CPPL_INT dgematrix::dggev(dgematrix& matB,
638  std::vector<double>& wr, std::vector<double>& wi)
639 {CPPL_VERBOSE_REPORT;
640 #ifdef CPPL_DEBUG
641  if(m!=n){
642  ERROR_REPORT;
643  std::cerr << "This matrix is not a square matrix." << std::endl
644  << "This matrix is (" << m << "x" << n << ")." << std::endl;
645  exit(1);
646  }
647  if(matB.m!=n || matB.n!=n){
648  ERROR_REPORT;
649  std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
650  << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
651  exit(1);
652  }
653 #endif//CPPL_DEBUG
654 
655  wr.resize(n);
656  wi.resize(n);
657  char JOBVL('n'), JOBVR('n');
658  CPPL_INT LDA(n), LDB(n), LDVL(1), LDVR(1), LWORK(8*n), INFO(1);
659  double *BETA(new double[n]), *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
660  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
661  delete [] WORK;
662  delete [] VL;
663  delete [] VR;
664 
665  //// reforming ////
666  for(CPPL_INT i=0; i<n; i++){
667  wr[i]/=BETA[i];
668  wi[i]/=BETA[i];
669  }
670  delete [] BETA;
671 
672  if(INFO!=0){
673  WARNING_REPORT;
674  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
675  }
676  return INFO;
677 }
678 
679 //=============================================================================
680 /*! calculate generalized eigenvalues and generalized right eigenvectors\n
681  All of the arguments don't need to be initialized.
682  wr, wi, vrr and vri are overwitten and become
683  real and imaginary part of generalized eigenvalue
684  and generalized right eigenvector, respectively.
685  This matrix and matB are also overwritten.
686 */
687 inline CPPL_INT dgematrix::dggev(dgematrix& matB,
688  std::vector<double>& wr, std::vector<double>& wi,
689  std::vector<dcovector>& vrr,
690  std::vector<dcovector>& vri)
691 {CPPL_VERBOSE_REPORT;
692 #ifdef CPPL_DEBUG
693  if(m!=n){
694  ERROR_REPORT;
695  std::cerr << "This matrix is not a square matrix." << std::endl
696  << "This matrix is (" << m << "x" << n << ")." << std::endl;
697  exit(1);
698  }
699  if(matB.m!=n || matB.n!=n){
700  ERROR_REPORT;
701  std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
702  << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
703  exit(1);
704  }
705 #endif//CPPL_DEBUG
706 
707  wr.resize(n);
708  wi.resize(n);
709  vrr.resize(n);
710  vri.resize(n);
711  for(CPPL_INT i=0; i<n; i++){
712  vrr[i].resize(n);
713  vri[i].resize(n);
714  }
715 
716  dgematrix VR(n,n);
717  char JOBVL('n'), JOBVR('V');
718  CPPL_INT LDA(n), LDB(n), LDVL(1), LDVR(n), LWORK(8*n), INFO(1);
719  double *BETA(new double[n]), *VL(NULL), *WORK(new double[LWORK]);
720  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL, &LDVL, VR.array, &LDVR, WORK, &LWORK, &INFO);
721  delete [] WORK;
722  delete [] VL;
723 
724  //// reforming ////
725  for(CPPL_INT i=0; i<n; i++){
726  wr[i]/=BETA[i];
727  wi[i]/=BETA[i];
728  }
729  delete [] BETA;
730 
731  //// forming ////
732  for(CPPL_INT j=0; j<n; j++){
733  if(fabs(wi[j])<DBL_MIN){
734  for(CPPL_INT i=0; i<n; i++){
735  vrr[j](i) = VR(i,j);
736  vri[j](i) = 0.0;
737  }
738  }
739  else{
740  for(CPPL_INT i=0; i<n; i++){
741  vrr[j](i) = VR(i,j);
742  vri[j](i) = VR(i,j+1);
743  vrr[j+1](i) = VR(i,j);
744  vri[j+1](i) =-VR(i,j+1);
745  }
746  j++;
747  }
748  }
749 
750  if(INFO!=0){
751  WARNING_REPORT;
752  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
753  }
754  return INFO;
755 }
756 
757 //=============================================================================
758 /*! calculate generalized eigenvalues and generalized left eigenvectors\n
759  All of the arguments don't need to be initialized.
760  wr, wi, vlr and vli are overwitten and become
761  real and imaginary part of generalized eigenvalue
762  and generalized left eigenvector, respectively.
763  This matrix and matB are also overwritten.
764 */
765 inline CPPL_INT dgematrix::dggev(dgematrix& matB,
766  std::vector<double>& wr, std::vector<double>& wi,
767  std::vector<drovector>& vlr,
768  std::vector<drovector>& vli)
769 {CPPL_VERBOSE_REPORT;
770 #ifdef CPPL_DEBUG
771  if(m!=n){
772  ERROR_REPORT;
773  std::cerr << "This matrix is not a square matrix." << std::endl
774  << "This matrix is (" << m << "x" << n << ")." << std::endl;
775  exit(1);
776  }
777  if(matB.m!=n || matB.n!=n){
778  ERROR_REPORT;
779  std::cerr << "The matrix B is not a square matrix having the same size as \"this\" matrix." << std::endl
780  << "The B matrix is (" << matB.m << "x" << matB.n << ")." << std::endl;
781  exit(1);
782  }
783 #endif//CPPL_DEBUG
784 
785  wr.resize(n);
786  wi.resize(n);
787  vlr.resize(n);
788  vli.resize(n);
789  for(CPPL_INT i=0; i<n; i++){
790  vlr[i].resize(n);
791  vli[i].resize(n);
792  }
793 
794  dgematrix VL(n,n);
795  char JOBVL('V'), JOBVR('n');
796  CPPL_INT LDA(n), LDB(n), LDVL(n), LDVR(1), LWORK(8*n), INFO(1);
797  double *BETA(new double[n]), *VR(NULL), *WORK(new double[LWORK]);
798  dggev_(&JOBVL, &JOBVR, &n, array, &LDA, matB.array, &LDB, &wr[0], &wi[0], BETA, VL.array, &LDVL, VR, &LDVR, WORK, &LWORK, &INFO);
799  delete [] WORK;
800  delete [] VR;
801 
802  //// reforming ////
803  for(CPPL_INT i=0; i<n; i++){
804  wr[i]/=BETA[i];
805  wi[i]/=BETA[i];
806  }
807  delete [] BETA;
808 
809  //// forming ////
810  for(CPPL_INT j=0; j<n; j++){
811  if(fabs(wi[j])<DBL_MIN){
812  for(CPPL_INT i=0; i<n; i++){
813  vlr[j](i) = VL(i,j);
814  vli[j](i) = 0.0;
815  }
816  }
817  else{
818  for(CPPL_INT i=0; i<n; i++){
819  vlr[j](i) = VL(i,j);
820  vli[j](i) =-VL(i,j+1);
821  vlr[j+1](i) = VL(i,j);
822  vli[j+1](i) = VL(i,j+1);
823  }
824  j++;
825  }
826  }
827 
828  if(INFO!=0){
829  WARNING_REPORT;
830  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
831  }
832  return INFO;
833 }
834 
835 ///////////////////////////////////////////////////////////////////////////////
836 ///////////////////////////////////////////////////////////////////////////////
837 ///////////////////////////////////////////////////////////////////////////////
838 
839 //=============================================================================
840 /*! compute the singular value decomposition (SVD)\n
841  The argument is dgbmatrix S.
842  S doesn't need to be initialized.
843  S is overwitten and become singular values.
844  This matrix also overwritten.
845 */
846 inline CPPL_INT dgematrix::dgesvd
847 (
848  dgbmatrix& S
849  )
850 {CPPL_VERBOSE_REPORT;
851  char JOBU ='N';
852  char JOBVT ='N';
853  // M
854  // N
855  // A
856  CPPL_INT LDA =m;
857  double* U =NULL;
858  S.resize(m,n,0,0);
859  CPPL_INT LDU =1;
860  double* VT =NULL;
861  CPPL_INT LDVT =1;
862  CPPL_INT LWORK =std::max(3*std::min(m,n)+std::max(m,n),5*std::min(m,n));
863  double *WORK =new double[LWORK];
864  CPPL_INT INFO =1;
865 
866  dgesvd_(&JOBU, &JOBVT, &m, &n, array, &LDA, S.array, U, &LDU, VT, &LDVT, WORK, &LWORK, &INFO);
867  delete [] WORK;
868 
869  if(INFO!=0){
870  WARNING_REPORT;
871  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
872  }
873  return INFO;
874 }
875 
876 //=============================================================================
877 /*! compute the singular value decomposition (SVD)\n
878  The arguments are dcocector S, dgematrix U and VT.
879  All of them need not to be initialized.
880  S, U and VT are overwitten and become singular values,
881  left singular vectors,
882  and right singular vectors respectively.
883  This matrix also overwritten.
884 */
885 inline CPPL_INT dgematrix::dgesvd(dcovector& S, dgematrix& U, dgematrix& VT)
886 {CPPL_VERBOSE_REPORT;
887  char JOBU('A'), JOBVT('A');
888  CPPL_INT LDA(m), LDU(m), LDVT(n), LWORK(std::max(3*std::min(m,n)+std::max(m,n),5*std::min(m,n))), INFO(1);
889  double *WORK(new double[LWORK]);
890  S.resize(std::min(m,n));
891  U.resize(LDU,m);
892  VT.resize(LDVT,n);
893 
894  dgesvd_(&JOBU, &JOBVT, &m, &n, array, &LDA, S.array, U.array, &LDU, VT.array, &LDVT, WORK, &LWORK, &INFO);
895  delete [] WORK;
896 
897  if(INFO!=0){
898  WARNING_REPORT;
899  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
900  }
901  return INFO;
902 }
903 
904 ///////////////////////////////////////////////////////////////////////////////
905 ///////////////////////////////////////////////////////////////////////////////
906 ///////////////////////////////////////////////////////////////////////////////
907 
908 //=============================================================================
909 /*! solve the linear equality-constrained least squares (LSE) problem\n
910  Input matrix and vectors, B, c, and d, are overwitten.
911  This matrix is also overwritten.
912  The solution vector x is to be automatically resized.
913 */
914 inline CPPL_INT dgematrix::dgglse
915 (
916  dgematrix& B,
917  dcovector& c,
918  dcovector& d,
919  dcovector& x
920  )
921 {CPPL_VERBOSE_REPORT;
922 #ifdef CPPL_DEBUG
923  if(m!=c.l){
924  ERROR_REPORT;
925  std::cerr << "A.m and c.l should be the same." << std::endl
926  << "Your input was A.m=" << m << " and c.l=" << c.l << std::endl;
927  exit(1);
928  }
929  if(B.m!=d.l){
930  ERROR_REPORT;
931  std::cerr << "B.m and d.l should be the same." << std::endl
932  << "Your input was B.m=" << B.m << " and d.l=" << d.l << std::endl;
933  exit(1);
934  }
935  if( !(B.m<=n) || !(n<=m+B.m) ){
936  ERROR_REPORT;
937  std::cerr << "B.m<=A.n<=A.m+B.m should be satisfied." << std::endl
938  << "Your input was B.m=" << B.m << ", A.n=" << n << ", and A.m+B.m=" << m+B.m << std::endl;
939  exit(1);
940  }
941 #endif//CPPL_DEBUG
942 
943  CPPL_INT lwork(-1), info(1);
944  dcovector work(1);
945  x.resize(n);
946 
947  //////// workspace query ////////
948  CPPL_INT lda =std::max(CPPL_INT(1),m);
949  CPPL_INT ldb =std::max(CPPL_INT(1),B.m);
950  dgglse_(&m, &n, &B.m, array, &lda, B.array, &ldb, c.array, d.array, x.array, work.array, &lwork, &info);
951  lwork =CPPL_INT(work(0));
952  work.resize(lwork);
953  info =1;
954 
955  //////// solve ////////
956  dgglse_(&m, &n, &B.m, array, &lda, B.array, &ldb, c.array, d.array, x.array, work.array, &lwork, &info);
957  work.clear();
958 
959  if(info!=0){
960  WARNING_REPORT;
961  std::cerr << "Serious trouble happend. info = " << info << "." << std::endl;
962  }
963  return info;
964 }
void clear()
dgematrix & resize(const CPPL_INT &, const CPPL_INT &)
CPPL_INT m
matrix row size
Definition: dgematrix.hpp:9
double * array
1D array to store matrix data
Definition: dgematrix.hpp:11
CPPL_INT dgeev(std::vector< double > &, std::vector< double > &)
CPPL_INT l
vector size
Definition: dcovector.hpp:9
_dgematrix i(const _dgbmatrix &mat)
dcovector & resize(const CPPL_INT &, const CPPL_INT=0)
friend void swap(dgematrix &, dgematrix &)
CPPL_INT dggev(dgematrix &, std::vector< double > &, std::vector< double > &)
CPPL_INT n
matrix column size
Definition: dgematrix.hpp:10
CPPL_INT l
vector size
Definition: drovector.hpp:9
Real Double-precision General Dence Matrix Class.
Definition: dgematrix.hpp:3
double * array
1D array to store vector data
Definition: dcovector.hpp:11
CPPL_INT dgesvd(dgbmatrix &)
double * array
1D array to store matrix data
Definition: dgbmatrix.hpp:13
CPPL_INT dgesv(dgematrix &)
void clear()
Real Double-precision Row Vector Class.
Definition: drovector.hpp:3
CPPL_INT dgelss(dcovector &, dcovector &, CPPL_INT &, const double)
void swap(dcovector &u, dcovector &v)
drovector & zero()
CPPL_INT dgglse(dgematrix &, dcovector &, dcovector &, dcovector &)
dgbmatrix & resize(const CPPL_INT &, const CPPL_INT &, const CPPL_INT &, const CPPL_INT &)
drovector & resize(const CPPL_INT &, const CPPL_INT=0)
Real Double-precision General Band Matrix Class.
Definition: dgbmatrix.hpp:3
CPPL_INT dgels(dgematrix &)
CPPL_INT dgelsd(dcovector &, dcovector &, CPPL_INT &, const double)
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
friend _dgematrix i(const dgematrix &)
void resize(const CPPL_INT &)
Complex Double-precision Column Vector Class.
Definition: zcovector.hpp:3