CPPLapack
 All Classes Files Functions Variables Friends Pages
zgematrix-lapack.hpp
Go to the documentation of this file.
1 //=============================================================================
2 /*! solve A*X=Y using zgesv\n
3  The argument is zgematrix Y. Y is overwritten and become the solution X.
4  A is also overwritten and become P*l*U. */
5 inline CPPL_INT zgematrix::zgesv(zgematrix& mat)
6 {CPPL_VERBOSE_REPORT;
7 #ifdef CPPL_DEBUG
8  if(m!=n || m!=mat.m){
9  ERROR_REPORT;
10  std::cerr << "These two matrices cannot be solved." << std::endl
11  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
12  exit(1);
13  }
14 #endif//CPPL_DEBUG
15 
16  CPPL_INT NRHS(mat.n), LDA(n), *IPIV(new CPPL_INT[n]), LDB(mat.m), INFO(1);
17  zgesv_(&n, &NRHS, array, &LDA, IPIV, mat.array, &LDB, &INFO);
18  delete [] IPIV;
19 
20  if(INFO!=0){
21  WARNING_REPORT;
22  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
23  }
24  return INFO;
25 }
26 
27 //=============================================================================
28 /*! solve A*x=y using zgesv\n
29  The argument is zcovector y. y is overwritten and become the solution x.
30  A is also overwritten and become P*l*U. */
31 inline CPPL_INT zgematrix::zgesv(zcovector& vec)
32 {CPPL_VERBOSE_REPORT;
33 #ifdef CPPL_DEBUG
34  if(m!=n || m!=vec.l){
35  ERROR_REPORT;
36  std::cerr << "These matrix and vector cannot be solved." << std::endl
37  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
38  exit(1);
39  }
40 #endif//CPPL_DEBUG
41 
42  CPPL_INT NRHS(1), LDA(n), *IPIV(new CPPL_INT[n]), LDB(vec.l), INFO(1);
43  zgesv_(&n, &NRHS, array, &LDA, IPIV, vec.array, &LDB, &INFO);
44  delete [] IPIV;
45 
46  if(INFO!=0){
47  WARNING_REPORT;
48  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
49  }
50  return INFO;
51 }
52 
53 ///////////////////////////////////////////////////////////////////////////////
54 ///////////////////////////////////////////////////////////////////////////////
55 ///////////////////////////////////////////////////////////////////////////////
56 
57 //=============================================================================
58 /*! solve overdetermined or underdetermined A*X=Y using zgels\n*/
59 inline CPPL_INT zgematrix::zgels(zgematrix& mat)
60 {CPPL_VERBOSE_REPORT;
61 #ifdef CPPL_DEBUG
62  if(m!=mat.m){
63  ERROR_REPORT;
64  std::cerr << "These two matrices cannot be solved." << std::endl
65  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
66  exit(1);
67  }
68 #endif//CPPL_DEBUG
69 
70  if(m<n){ //underdetermined
71  zgematrix tmp(n,mat.n);
72  for(CPPL_INT i=0; i<mat.m; i++){
73  for(CPPL_INT j=0; j<mat.n; j++){
74  tmp(i,j) =mat(i,j);
75  }
76  }
77  mat.clear();
78  swap(mat,tmp);
79  }
80 
81  char TRANS('n');
82  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
83  comple *WORK(new comple[LWORK]);
84  zgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
85  delete [] WORK;
86 
87  if(m>n){ //overdetermined
88  zgematrix tmp(n,mat.n);
89  for(CPPL_INT i=0; i<tmp.m; i++){
90  for(CPPL_INT j=0; j<tmp.n; j++){
91  tmp(i,j) =mat(i,j);
92  }
93  }
94  mat.clear();
95  swap(mat,tmp);
96  }
97 
98  if(INFO!=0){
99  WARNING_REPORT;
100  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
101  }
102  return INFO;
103 }
104 
105 //=============================================================================
106 /*! solve overdetermined or underdetermined A*x=y using zgels\n*/
107 inline CPPL_INT zgematrix::zgels(zcovector& vec)
108 {CPPL_VERBOSE_REPORT;
109 #ifdef CPPL_DEBUG
110  if(m!=vec.l){
111  ERROR_REPORT;
112  std::cerr << "These matrix and vector cannot be solved." << std::endl
113  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
114  exit(1);
115  }
116 #endif//CPPL_DEBUG
117 
118  if(m<n){ //underdetermined
119  zcovector tmp(n);
120  for(CPPL_INT i=0; i<vec.l; i++){
121  tmp(i)=vec(i);
122  }
123  vec.clear();
124  swap(vec,tmp);
125  }
126 
127  char TRANS('n');
128  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
129  comple *WORK(new comple[LWORK]);
130  zgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
131  delete [] WORK;
132 
133  if(m>n){ //overdetermined
134  zcovector tmp(n);
135  for(CPPL_INT i=0; i<tmp.l; i++){
136  tmp(i)=vec(i);
137  }
138  vec.clear();
139  swap(vec,tmp);
140  }
141 
142  if(INFO!=0){
143  WARNING_REPORT;
144  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
145  }
146  return INFO;
147 }
148 
149 //=============================================================================
150 /*! solve overdetermined or underdetermined A*X=Y using zgels
151  with the sum of residual squares output\n
152  The residual is set as the columnwise sum of residual squares
153  for overdetermined problems
154  while it is always zero for underdetermined problems.
155 */
156 inline CPPL_INT zgematrix::zgels(zgematrix& mat, drovector& residual)
157 {CPPL_VERBOSE_REPORT;
158 #ifdef CPPL_DEBUG
159  if(m!=mat.m){
160  ERROR_REPORT;
161  std::cerr << "These two matrices cannot be solved." << std::endl
162  << "Your input was (" << m << "x" << n << ") and (" << mat.m << "x" << mat.n << ")." << std::endl;
163  exit(1);
164  }
165 #endif//CPPL_DEBUG
166 
167  residual.resize(mat.n);
168  residual.zero();
169 
170  if(m<n){ //underdetermined
171  zgematrix tmp(n,mat.n);
172  for(CPPL_INT i=0; i<mat.m; i++){
173  for(CPPL_INT j=0; j<mat.n; j++){
174  tmp(i,j) =mat(i,j);
175  }
176  }
177  mat.clear();
178  swap(mat,tmp);
179  }
180 
181  char TRANS('n');
182  CPPL_INT NRHS(mat.n), LDA(m), LDB(mat.m), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
183  comple *WORK(new comple[LWORK]);
184  zgels_(&TRANS, &m, &n, &NRHS, array, &LDA, mat.array, &LDB, WORK, &LWORK, &INFO);
185  delete [] WORK;
186 
187  if(m>n){ //overdetermined
188  for(CPPL_INT i=0; i<residual.l; i++){
189  for(CPPL_INT j=0; j<m-n; j++){
190  residual(i) += std::norm(mat(n+j,i));
191  }
192  }
193 
194  zgematrix tmp(n,mat.n);
195  for(CPPL_INT i=0; i<tmp.m; i++){
196  for(CPPL_INT j=0; j<tmp.n; j++){
197  tmp(i,j) =mat(i,j);
198  }
199  }
200  mat.clear();
201  swap(mat,tmp);
202  }
203 
204  if(INFO!=0){
205  WARNING_REPORT;
206  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
207  }
208  return INFO;
209 }
210 
211 //=============================================================================
212 /*! solve overdetermined or underdetermined A*x=y using zgels
213  with the sum of residual squares output\n
214  The residual is set as the sum of residual squares for overdetermined problems
215  while it is always zero for underdetermined problems.
216 */
217 inline CPPL_INT zgematrix::zgels(zcovector& vec, double& residual)
218 {CPPL_VERBOSE_REPORT;
219 #ifdef CPPL_DEBUG
220  if(m!=vec.l){
221  ERROR_REPORT;
222  std::cerr << "These matrix and vector cannot be solved." << std::endl
223  << "Your input was (" << m << "x" << n << ") and (" << vec.l << ")." << std::endl;
224  exit(1);
225  }
226 #endif//CPPL_DEBUG
227 
228  residual=0.0;
229 
230  if(m<n){ //underdetermined
231  zcovector tmp(n);
232  for(CPPL_INT i=0; i<vec.l; i++){
233  tmp(i)=vec(i);
234  }
235  vec.clear();
236  swap(vec,tmp);
237  }
238 
239  char TRANS('n');
240  CPPL_INT NRHS(1), LDA(m), LDB(vec.l), LWORK(std::min(m,n)+std::max(std::min(m,n),NRHS)), INFO(1);
241  comple *WORK(new comple[LWORK]);
242  zgels_(&TRANS, &m, &n, &NRHS, array, &LDA, vec.array, &LDB, WORK, &LWORK, &INFO);
243  delete [] WORK;
244 
245  if(m>n){ //overdetermined
246  for(CPPL_INT i=0; i<m-n; i++){
247  residual += std::norm(vec(n+i));
248  }
249 
250  zcovector tmp(n);
251  for(CPPL_INT i=0; i<tmp.l; i++){
252  tmp(i) =vec(i);
253  }
254  vec.clear();
255  swap(vec,tmp);
256  }
257 
258  if(INFO!=0){
259  WARNING_REPORT;
260  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
261  }
262  return INFO;
263 }
264 
265 ///////////////////////////////////////////////////////////////////////////////
266 ///////////////////////////////////////////////////////////////////////////////
267 ///////////////////////////////////////////////////////////////////////////////
268 
269 //=============================================================================
270 /*! calculate the least-squares-least-norm solution
271  for overdetermined or underdetermined A*x=y using zgelss\n */
272 inline CPPL_INT zgematrix::zgelss(zcovector& B, dcovector& S, CPPL_INT& RANK,
273  const double RCOND =-1. )
274 {CPPL_VERBOSE_REPORT;
275 #ifdef CPPL_DEBUG
276  if(m!=B.l){
277  ERROR_REPORT;
278  std::cerr << "These matrix and vector cannot be solved." << std::endl
279  << "Your input was (" << m << "x" << n << ") and (" << B.l << ")." << std::endl;
280  exit(1);
281  }
282 #endif//CPPL_DEBUG
283 
284  if(m<n){ //underdetermined
285  zcovector tmp(n);
286  for(CPPL_INT i=0; i<B.l; i++){
287  tmp(i)=B(i);
288  }
289  B.clear();
290  swap(B,tmp);
291  }
292 
293  S.resize(std::min(m,n));
294 
295  CPPL_INT NRHS(1), LDA(m), LDB(B.l), LWORK(2*std::min(m,n) +std::max(std::max(m,n),NRHS)), INFO(1);
296  double *RWORK(new double[5*std::min(m,n)]);
297  comple *WORK(new comple[LWORK]);
298  zgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, RWORK, &INFO);
299  delete [] RWORK;
300  delete [] WORK;
301 
302  if(INFO!=0){
303  WARNING_REPORT;
304  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
305  }
306  return INFO;
307 }
308 
309 //=============================================================================
310 /*! calculate the least-squares-least-norm solution
311  for overdetermined or underdetermined A*x=y using zgelss\n */
312 inline CPPL_INT zgematrix::zgelss(zgematrix& B, dcovector& S, CPPL_INT& RANK,
313  const double RCOND =-1. )
314 {CPPL_VERBOSE_REPORT;
315 #ifdef CPPL_DEBUG
316  if(m!=B.m){
317  ERROR_REPORT;
318  std::cerr << "These matrix and vector cannot be solved." << std::endl
319  << "Your input was (" << m << "x" << n << ") and (" << B.m << "x" << B.n << ")." << std::endl;
320  exit(1);
321  }
322 #endif//CPPL_DEBUG
323 
324  if(m<n){ //underdetermined
325  zgematrix tmp(n,B.n);
326  for(CPPL_INT i=0; i<B.m; i++){
327  for(CPPL_INT j=0; j<B.n; j++){
328  tmp(i,j)=B(i,j);
329  }
330  }
331  B.clear();
332  swap(B,tmp);
333  }
334 
335  S.resize(std::min(m,n));
336 
337  CPPL_INT NRHS(B.n), LDA(m), LDB(B.m), LWORK(2*std::min(m,n) +std::max(std::max(m,n),NRHS)), INFO(1);
338  double *RWORK(new double[5*std::min(m,n)]);
339  comple *WORK(new comple[LWORK]);
340  zgelss_(&m, &n, &NRHS, array, &LDA, B.array, &LDB, S.array, &RCOND, &RANK, WORK, &LWORK, RWORK, &INFO);
341  delete [] RWORK;
342  delete [] WORK;
343 
344  if(INFO!=0){
345  WARNING_REPORT;
346  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
347  }
348  return INFO;
349 }
350 
351 ///////////////////////////////////////////////////////////////////////////////
352 ///////////////////////////////////////////////////////////////////////////////
353 ///////////////////////////////////////////////////////////////////////////////
354 
355 //=============================================================================
356 /*! calculate eigenvalues\n
357  The argument need not to be initialized.
358  w is overwitten and become eigenvalues.
359  This matrix is also overwritten. */
360 inline CPPL_INT zgematrix::zgeev(std::vector< comple >& w)
361 {CPPL_VERBOSE_REPORT;
362 #ifdef CPPL_DEBUG
363  if(m!=n){
364  ERROR_REPORT;
365  std::cerr << "This matrix cannot have eigenvalues." << std::endl
366  << "Your input was (" << m << "x" << n << ")." << std::endl;
367  exit(1);
368  }
369 #endif//CPPL_DEBUG
370 
371  w.resize(n);
372 
373  char JOBVL('n'), JOBVR('n');
374  CPPL_INT LDA(n), LDVL(1), LDVR(1), LWORK(4*n), INFO(1);
375  double *RWORK(new double[2*n]);
376  comple *VL(NULL), *VR(NULL), *WORK(new comple[LWORK]);
377  zgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &w[0], VL, &LDVL, VR, &LDVR, WORK, &LWORK, RWORK, &INFO);
378  delete [] RWORK;
379  delete [] WORK;
380  delete [] VL;
381  delete [] VR;
382 
383  if(INFO!=0){
384  WARNING_REPORT;
385  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
386  }
387  return INFO;
388 }
389 
390 //=============================================================================
391 /*! calculate eigenvalues and right eigenvectors\n
392  All of the arguments need not to be initialized.
393  w, vr are overwitten and become eigenvalues and right eigenvectors,
394  respectively.
395  This matrix is also overwritten. */
396 inline CPPL_INT zgematrix::zgeev(std::vector< comple >& w,
397  std::vector<zcovector>& vr)
398 {CPPL_VERBOSE_REPORT;
399 #ifdef CPPL_DEBUG
400  if(m!=n){
401  ERROR_REPORT;
402  std::cerr << "This matrix cannot have eigenvalues." << std::endl
403  << "Your input was (" << m << "x" << n << ")." << std::endl;
404  exit(1);
405  }
406 #endif//CPPL_DEBUG
407 
408  w.resize(n);
409  vr.resize(n);
410  for(CPPL_INT i=0; i<n; i++){
411  vr[i].resize(n);
412  }
413 
414  zgematrix VR(n,n);
415  char JOBVL('n'), JOBVR('V');
416  CPPL_INT LDA(n), LDVL(1), LDVR(n), LWORK(4*n), INFO(1);
417  double *RWORK(new double[2*n]);
418  comple *VL(NULL), *WORK(new comple[LWORK]);
419  zgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &w[0], VL, &LDVL, VR.array, &LDVR, WORK, &LWORK, RWORK, &INFO);
420  delete [] RWORK;
421  delete [] WORK;
422  delete [] VL;
423 
424  //// forming ////
425  for(CPPL_INT j=0; j<n; j++){
426  for(CPPL_INT i=0; i<n; i++){
427  vr[j](i) = VR(i,j);
428  }
429  }
430 
431  if(INFO!=0){
432  WARNING_REPORT;
433  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
434  }
435  return INFO;
436 }
437 
438 //=============================================================================
439 /*! calculate eigenvalues and left eigenvectors\n
440  All of the arguments need not to be initialized.
441  w, vr are overwitten and become eigenvalues and left eigenvectors,
442  respectively.
443  This matrix is also overwritten. */
444 inline CPPL_INT zgematrix::zgeev(std::vector< comple >& w,
445  std::vector<zrovector>& vl)
446 {CPPL_VERBOSE_REPORT;
447 #ifdef CPPL_DEBUG
448  if(m!=n){
449  ERROR_REPORT;
450  std::cerr << "This matrix cannot have eigenvalues." << std::endl
451  << "Your input was (" << m << "x" << n << ")." << std::endl;
452  exit(1);
453  }
454 #endif//CPPL_DEBUG
455 
456  w.resize(n);
457  vl.resize(n);
458  for(CPPL_INT i=0; i<n; i++){
459  vl[i].resize(n);
460  }
461  zgematrix VL(n,n);
462 
463  char JOBVL('V'), JOBVR('n');
464  CPPL_INT LDA(n), LDVL(n), LDVR(1), LWORK(4*n), INFO(1);
465  double *RWORK(new double[2*n]);
466  comple *VR(NULL), *WORK(new comple[LWORK]);
467  zgeev_(&JOBVL, &JOBVR, &n, array, &LDA, &w[0], VL.array, &LDVL, VR, &LDVR, WORK, &LWORK, RWORK, &INFO);
468  delete [] RWORK;
469  delete [] WORK;
470  delete [] VR;
471 
472  //// forming ////
473  for(CPPL_INT j=0; j<n; j++){
474  for(CPPL_INT i=0; i<n; i++){
475  vl[j](i) = std::conj(VL(i,j));
476  }
477  }
478 
479  if(INFO!=0){
480  WARNING_REPORT;
481  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
482  }
483  return INFO;
484 }
485 
486 
487 ///////////////////////////////////////////////////////////////////////////////
488 ///////////////////////////////////////////////////////////////////////////////
489 ///////////////////////////////////////////////////////////////////////////////
490 
491 //=============================================================================
492 //inline CPPL_INT zgematrix::zgegv()
493 
494 
495 ///////////////////////////////////////////////////////////////////////////////
496 ///////////////////////////////////////////////////////////////////////////////
497 ///////////////////////////////////////////////////////////////////////////////
498 
499 //=============================================================================
500 /*! compute the singular value decomposition (SVD)\n
501  The arguments are zcocector S, zgematrix U and VT.
502  All of them need not to be initialized.
503  S, U and VT are overwitten and become singular values, left singular vectors,
504  and right singular vectors respectively.
505  This matrix also overwritten.
506 */
507 inline CPPL_INT zgematrix::zgesvd(dcovector& S, zgematrix& U, zgematrix& VT)
508 {CPPL_VERBOSE_REPORT;
509  char JOBU('A'), JOBVT('A');
510  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);
511  double *RWORK(new double[5*std::min(m,n)]);
512  comple *WORK(new comple[LWORK]);
513  S.resize(std::min(m,n));
514  U.resize(LDU,m);
515  VT.resize(LDVT,n);
516 
517  zgesvd_(&JOBU, &JOBVT, &m, &n, array, &LDA, S.array, U.array, &LDU, VT.array, &LDVT, WORK, &LWORK, RWORK, &INFO);
518  delete [] RWORK;
519  delete [] WORK;
520 
521  if(INFO!=0){
522  WARNING_REPORT;
523  std::cerr << "Serious trouble happend. INFO = " << INFO << "." << std::endl;
524  }
525  return INFO;
526 }
void clear()
comple * array
1D array to store vector data
Definition: zcovector.hpp:10
void resize(const CPPL_INT &, const CPPL_INT &)
CPPL_INT l
vector size
Definition: zcovector.hpp:9
dcovector & resize(const CPPL_INT &, const CPPL_INT=0)
CPPL_INT l
vector size
Definition: drovector.hpp:9
friend _zgematrix i(const zgematrix &)
double * array
1D array to store vector data
Definition: dcovector.hpp:11
_zcovector conj(const _zcovector &vec)
CPPL_INT zgels(zgematrix &)
CPPL_INT n
matrix column size
Definition: zgematrix.hpp:10
Complex Double-precision General Dence Matrix Class.
Definition: zgematrix.hpp:3
Real Double-precision Row Vector Class.
Definition: drovector.hpp:3
comple * array
1D array to store matrix data
Definition: zgematrix.hpp:11
CPPL_INT zgeev(std::vector< comple > &)
drovector & zero()
void clear()
drovector & resize(const CPPL_INT &, const CPPL_INT=0)
CPPL_INT zgesvd(dcovector &, zgematrix &, zgematrix &)
CPPL_INT m
matrix row size
Definition: zgematrix.hpp:9
friend void swap(zgematrix &, zgematrix &)
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
Complex Double-precision Column Vector Class.
Definition: zcovector.hpp:3
CPPL_INT zgelss(zcovector &, dcovector &, CPPL_INT &, const double)
CPPL_INT zgesv(zgematrix &)