My Project
zgematrix-lapack.hpp
1 //=============================================================================
5 inline long zgematrix::zgesv(zgematrix& mat)
6 {
7 #ifdef CPPL_VERBOSE
8  std::cerr << "# [MARK] zgematrix::zgesv(zgematrix&)"
9  << std::endl;
10 #endif//CPPL_VERBOSE
11 
12 #ifdef CPPL_DEBUG
13  if(M!=N || M!=mat.M){
14  std::cerr << "[ERROR] zgematrix::zgesv(zgematrix&) " << std::endl
15  << "These two matrices cannot be solved." << std::endl
16  << "Your input was (" << M << "x" << N << ") and ("
17  << mat.M << "x" << mat.N << ")." << std::endl;
18  exit(1);
19  }
20 #endif//CPPL_DEBUG
21 
22  long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.M), INFO(1);
23  zgesv_(N, NRHS, Array, LDA, IPIV, mat.Array, LDB, INFO);
24  delete [] IPIV;
25 
26  if(INFO!=0){
27  std::cerr << "[WARNING] zgematrix::zgesv(zgematrix&) "
28  << "Serious trouble happend. INFO = " << INFO << "."
29  << std::endl;
30  }
31  return INFO;
32 }
33 
34 //=============================================================================
38 inline long zgematrix::zgesv(zcovector& vec)
39 {
40 #ifdef CPPL_VERBOSE
41  std::cerr << "# [MARK] zgematrix::zgesv(zcovector&)"
42  << std::endl;
43 #endif//CPPL_VERBOSE
44 
45 #ifdef CPPL_DEBUG
46  if(M!=N || M!=vec.L){
47  std::cerr << "[ERROR] zgematrix::zgesv(zcovector&) " << std::endl
48  << "These matrix and vector cannot be solved." << std::endl
49  << "Your input was (" << M << "x" << N << ") and ("
50  << vec.L << ")." << std::endl;
51  exit(1);
52  }
53 #endif//CPPL_DEBUG
54  long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), INFO(1);
55  zgesv_(N, NRHS, Array, LDA, IPIV, vec.Array, LDB, INFO);
56  delete [] IPIV;
57 
58  if(INFO!=0){
59  std::cerr << "[WARNING] zgematrix::zgesv(zcovector&) "
60  << "Serious trouble happend. INFO = " << INFO << "."
61  << std::endl;
62  }
63  return INFO;
64 }
65 
69 
70 //=============================================================================
72 inline long zgematrix::zgels(zgematrix& mat)
73 {
74 #ifdef CPPL_VERBOSE
75  std::cerr << "# [MARK] zgematrix::zgels(zgematrix&)"
76  << std::endl;
77 #endif//CPPL_VERBOSE
78 
79 #ifdef CPPL_DEBUG
80  if(M!=mat.M){
81  std::cerr << "[ERROR] zgematrix::zgels(zgematrix&) " << std::endl
82  << "These two matrices cannot be solved." << std::endl
83  << "Your input was (" << M << "x" << N << ") and ("
84  << mat.M << "x" << mat.N << ")." << std::endl;
85  exit(1);
86  }
87 #endif//CPPL_DEBUG
88 
89  if(M<N){ //underdetermined
90  zgematrix tmp(N,mat.N);
91  for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
92  tmp(i,j) =mat(i,j);
93  }}
94  mat.clear();
95  swap(mat,tmp);
96  }
97 
98  char TRANS('N');
99  long NRHS(mat.N), LDA(M), LDB(mat.M),
100  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
101  std::complex<double> *WORK(new std::complex<double>[LWORK]);
102  zgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
103  delete [] WORK;
104 
105  if(M>N){ //overdetermined
106  zgematrix tmp(N,mat.N);
107  for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
108  tmp(i,j) =mat(i,j);
109  }}
110  mat.clear();
111  swap(mat,tmp);
112  }
113 
114  if(INFO!=0){
115  std::cerr << "[WARNING] zgematrix::zgels(zgematrix&) "
116  << "Serious trouble happend. INFO = " << INFO << "."
117  << std::endl;
118  }
119  return INFO;
120 }
121 
122 //=============================================================================
124 inline long zgematrix::zgels(zcovector& vec)
125 {
126 #ifdef CPPL_VERBOSE
127  std::cerr << "# [MARK] zgematrix::zgels(zcovector&)"
128  << std::endl;
129 #endif//CPPL_VERBOSE
130 
131 #ifdef CPPL_DEBUG
132  if(M!=vec.L){
133  std::cerr << "[ERROR] zgematrix::zgels(zcovector&) " << std::endl
134  << "These matrix and vector cannot be solved." << std::endl
135  << "Your input was (" << M << "x" << N << ") and ("
136  << vec.L << ")." << std::endl;
137  exit(1);
138  }
139 #endif//CPPL_DEBUG
140 
141  if(M<N){ //underdetermined
142  zcovector tmp(N);
143  for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
144  vec.clear();
145  swap(vec,tmp);
146  }
147 
148  char TRANS('N');
149  long NRHS(1), LDA(M), LDB(vec.L),
150  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
151  std::complex<double> *WORK(new std::complex<double>[LWORK]);
152  zgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
153  delete [] WORK;
154 
155  if(M>N){ //overdetermined
156  zcovector tmp(N);
157  for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
158  vec.clear();
159  swap(vec,tmp);
160  }
161 
162  if(INFO!=0){
163  std::cerr << "[WARNING] zgematrix::zgels(zcovector&) "
164  << "Serious trouble happend. INFO = " << INFO << "."
165  << std::endl;
166  }
167  return INFO;
168 }
169 
170 //=============================================================================
177 inline long zgematrix::zgels(zgematrix& mat, drovector& residual)
178 {
179 #ifdef CPPL_VERBOSE
180  std::cerr << "# [MARK] zgematrix::zgels(zgematrix&, drovector&)"
181  << std::endl;
182 #endif//CPPL_VERBOSE
183 
184 #ifdef CPPL_DEBUG
185  if(M!=mat.M){
186  std::cerr << "[ERROR] zgematrix::zgels(zgematrix&, drovector&) "
187  << std::endl
188  << "These two matrices cannot be solved." << std::endl
189  << "Your input was (" << M << "x" << N << ") and ("
190  << mat.M << "x" << mat.N << ")." << std::endl;
191  exit(1);
192  }
193 #endif//CPPL_DEBUG
194 
195  residual.resize(mat.N); residual.zero();
196 
197  if(M<N){ //underdetermined
198  zgematrix tmp(N,mat.N);
199  for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
200  tmp(i,j) =mat(i,j);
201  }}
202  mat.clear();
203  swap(mat,tmp);
204  }
205 
206  char TRANS('N');
207  long NRHS(mat.N), LDA(M), LDB(mat.M),
208  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
209  std::complex<double> *WORK(new std::complex<double>[LWORK]);
210  zgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
211  delete [] WORK;
212 
213  if(M>N){ //overdetermined
214  for(long i=0; i<residual.l; i++){ for(long j=0; j<M-N; j++){
215  residual(i) += std::norm(mat(N+j,i));
216  }}
217 
218  zgematrix tmp(N,mat.N);
219  for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
220  tmp(i,j) =mat(i,j);
221  }}
222  mat.clear();
223  swap(mat,tmp);
224  }
225 
226  if(INFO!=0){
227  std::cerr << "[WARNING] zgematrix::zgels(zgematrix&, drovector&) "
228  << "Serious trouble happend. INFO = " << INFO << "."
229  << std::endl;
230  }
231  return INFO;
232 }
233 
234 //=============================================================================
240 inline long zgematrix::zgels(zcovector& vec, double& residual)
241 {
242 #ifdef CPPL_VERBOSE
243  std::cerr << "# [MARK] zgematrix::zgels(zcovector&, double&)"
244  << std::endl;
245 #endif//CPPL_VERBOSE
246 
247 #ifdef CPPL_DEBUG
248  if(M!=vec.L){
249  std::cerr << "[ERROR] zgematrix::zgels(zcovector&, double&) " << std::endl
250  << "These matrix and vector cannot be solved." << std::endl
251  << "Your input was (" << M << "x" << N << ") and ("
252  << vec.L << ")." << std::endl;
253  exit(1);
254  }
255 #endif//CPPL_DEBUG
256 
257  residual=0.0;
258 
259  if(M<N){ //underdetermined
260  zcovector tmp(n);
261  for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
262  vec.clear();
263  swap(vec,tmp);
264  }
265 
266  char TRANS('N');
267  long NRHS(1), LDA(m), LDB(vec.L),
268  LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
269  std::complex<double> *WORK(new std::complex<double>[LWORK]);
270  zgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
271  delete [] WORK;
272 
273  if(M>N){ //overdetermined
274  for(long i=0; i<M-N; i++){ residual+=std::norm(vec(N+i)); }
275 
276  zcovector tmp(N);
277  for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
278  vec.clear();
279  swap(vec,tmp);
280  }
281 
282  if(INFO!=0){
283  std::cerr << "[WARNING] zgematrix::zgels(zcovector&, double&) "
284  << "Serious trouble happend. INFO = " << INFO << "."
285  << std::endl;
286  }
287  return INFO;
288 }
289 
293 
294 //=============================================================================
297 inline long zgematrix::zgelss(zcovector& B, dcovector& S, long& RANK,
298  const double RCOND =-1. )
299 {
300 #ifdef CPPL_VERBOSE
301  std::cerr << "# [MARK] zgematrix::zgelss(zcovector&, dcovector&, long&, const double)"
302  << std::endl;
303 #endif//CPPL_VERBOSE
304 
305 #ifdef CPPL_DEBUG
306  if(M!=B.L){
307  std::cerr << "[ERROR] zgematrix::zgelss"
308  << "(zcovector&, dcovector&, long&, const double) " << std::endl
309  << "These matrix and vector cannot be solved." << std::endl
310  << "Your input was (" << M << "x" << N << ") and ("
311  << B.L << ")." << std::endl;
312  exit(1);
313  }
314 #endif//CPPL_DEBUG
315 
316  if(M<N){ //underdetermined
317  zcovector tmp(N);
318  for(long i=0; i<B.L; i++){ tmp(i)=B(i); }
319  B.clear();
320  swap(B,tmp);
321  }
322 
323  S.resize(min(M,N));
324 
325  long NRHS(1), LDA(M), LDB(B.L),
326  LWORK(2*min(M,N) +max(max(M,N),NRHS)), INFO(1);
327  double *RWORK(new double[5*min(M,N)]);
328  std::complex<double> *WORK(new std::complex<double>[LWORK]);
329  zgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.array, RCOND, RANK,
330  WORK, LWORK, RWORK, INFO);
331  delete [] RWORK; delete [] WORK;
332 
333  if(INFO!=0){
334  std::cerr << "[WARNING] zgematrix::zgelss"
335  << "(zcovector&, dcovector&, long&, double) "
336  << "Serious trouble happend. INFO = " << INFO << "."
337  << std::endl;
338  }
339  return INFO;
340 }
341 
342 //=============================================================================
345 inline long zgematrix::zgelss(zgematrix& B, dcovector& S, long& RANK,
346  const double RCOND =-1. )
347 {
348 #ifdef CPPL_VERBOSE
349  std::cerr << "# [MARK] zgematrix::zgelss(zgematrix&, dcovector&, long&, const double)"
350  << std::endl;
351 #endif//CPPL_VERBOSE
352 
353 #ifdef CPPL_DEBUG
354  if(M!=B.M){
355  std::cerr << "[ERROR] zgematrix::zgelss"
356  << "(zgematrix&, dcovector&, long&, const double) " << std::endl
357  << "These matrix and vector cannot be solved." << std::endl
358  << "Your input was (" << M << "x" << N << ") and ("
359  << B.M << "x" << B.N << ")." << std::endl;
360  exit(1);
361  }
362 #endif//CPPL_DEBUG
363 
364  if(M<N){ //underdetermined
365  zgematrix tmp(N,B.N);
366  for(long i=0; i<B.M; i++){
367  for(long j=0; j<B.N; j++){
368  tmp(i,j)=B(i,j);
369  }
370  }
371  B.clear();
372  swap(B,tmp);
373  }
374 
375  S.resize(min(M,N));
376 
377  long NRHS(B.N), LDA(M), LDB(B.M),
378  LWORK(2*min(M,N) +max(max(M,N),NRHS)), INFO(1);
379  double *RWORK(new double[5*min(M,N)]);
380  std::complex<double> *WORK(new std::complex<double>[LWORK]);
381  zgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.array, RCOND, RANK,
382  WORK, LWORK, RWORK, INFO);
383  delete [] RWORK; delete [] WORK;
384 
385  if(INFO!=0){
386  std::cerr << "[WARNING] zgematrix::zgelss"
387  << "(zgematrix&, dcovector&, long&, const double) "
388  << "Serious trouble happend. INFO = " << INFO << "."
389  << std::endl;
390  }
391  return INFO;
392 }
393 
397 
398 //=============================================================================
403 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w)
404 {
405 #ifdef CPPL_VERBOSE
406  std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&)"
407  << std::endl;
408 #endif//CPPL_VERBOSE
409 
410 #ifdef CPPL_DEBUG
411  if(M!=N){
412  std::cerr << "[ERROR] zgematrix::zgeev"
413  << "(std::vector<std::complex<double>>&) "
414  << std::endl
415  << "This matrix cannot have eigenvalues." << std::endl
416  << "Your input was (" << M << "x" << N << ")." << std::endl;
417  exit(1);
418  }
419 #endif//CPPL_DEBUG
420 
421  w.resize(N);
422  char JOBVL('N'), JOBVR('N');
423  long LDA(N), LDVL(1), LDVR(1), LWORK(4*n), INFO(1);
424  double *RWORK(new double[2*n]);
425  std::complex<double> *VL(NULL), *VR(NULL),
426  *WORK(new std::complex<double>[LWORK]);
427  zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
428  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO);
429  delete [] RWORK; delete [] WORK; delete [] VL; delete [] VR;
430 
431  if(INFO!=0){
432  std::cerr << "[WARNING] zgematrix::zgesv"
433  << "(std::vector<std::complex<double>&, std::vector<zcovector>&) "
434  << "Serious trouble happend. INFO = " << INFO << "."
435  << std::endl;
436  }
437  return INFO;
438 }
439 
440 //=============================================================================
446 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w,
447  std::vector<zcovector>& vr)
448 {
449 #ifdef CPPL_VERBOSE
450  std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&, std::vector<zcovector>&)"
451  << std::endl;
452 #endif//CPPL_VERBOSE
453 
454 #ifdef CPPL_DEBUG
455  if(M!=N){
456  std::cerr << "[ERROR] zgematrix::zgeev"
457  << "(std::vector<std::complex<double>>&, std::vector<zcovector>&)"
458  << std::endl
459  << "This matrix cannot have eigenvalues." << std::endl
460  << "Your input was (" << M << "x" << N << ")." << std::endl;
461  exit(1);
462  }
463 #endif//CPPL_DEBUG
464 
465  w.resize(N); vr.resize(N);
466  for(long i=0; i<N; i++){ vr[i].resize(N); }
467  zgematrix VR(N,N);
468  char JOBVL('N'), JOBVR('V');
469  long LDA(N), LDVL(1), LDVR(N), LWORK(4*n), INFO(1);
470  double *RWORK(new double[2*n]);
471  std::complex<double> *VL(NULL), *WORK(new std::complex<double>[LWORK]);
472  zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
473  VL, LDVL, VR.Array, LDVR, WORK, LWORK, RWORK, INFO);
474  delete [] RWORK; delete [] WORK; delete [] VL;
475 
477  for(long j=0; j<N; j++){ for(long i=0; i<N; i++){
478  vr[j](i) = VR(i,j);
479  }}
480 
481  if(INFO!=0){
482  std::cerr << "[WARNING] zgematrix::zgesv"
483  << "(std::vector<std::complex<double>&, std::vector<zcovector>&) "
484  << "Serious trouble happend. INFO = " << INFO << "."
485  << std::endl;
486  }
487  return INFO;
488 }
489 
490 //=============================================================================
496 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w,
497  std::vector<zrovector>& vl)
498 {
499 #ifdef CPPL_VERBOSE
500  std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&, std::vector<zrovector>&)"
501  << std::endl;
502 #endif//CPPL_VERBOSE
503 
504 #ifdef CPPL_DEBUG
505  if(M!=N){
506  std::cerr << "[ERROR] zgematrix::zgeev"
507  << "(std::vector<std::complex<double>>&, std::vector<zrovector>&)"
508  << std::endl
509  << "This matrix cannot have eigenvalues." << std::endl
510  << "Your input was (" << M << "x" << N << ")." << std::endl;
511  exit(1);
512  }
513 #endif//CPPL_DEBUG
514 
515  w.resize(N); vl.resize(N);
516  for(long i=0; i<N; i++){ vl[i].resize(N); }
517  zgematrix VL(N,N);
518  char JOBVL('V'), JOBVR('N');
519  long LDA(N), LDVL(N), LDVR(1), LWORK(4*n), INFO(1);
520  double *RWORK(new double[2*n]);
521  std::complex<double> *VR(NULL), *WORK(new std::complex<double>[LWORK]);
522  zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
523  VL.Array, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO);
524  delete [] RWORK; delete [] WORK; delete [] VR;
525 
527  for(long j=0; j<N; j++){ for(long i=0; i<N; i++){
528  vl[j](i) = std::conj(VL(i,j));
529  }}
530 
531  if(INFO!=0){
532  std::cerr << "[WARNING] zgematrix::zgesv"
533  << "(std::vector<std::complex<double>&, std::vector<zrovector>&) "
534  << "Serious trouble happend. INFO = " << INFO << "."
535  << std::endl;
536  }
537  return INFO;
538 }
539 
540 
544 
545 //=============================================================================
546 //inline long zgematrix::zgegv()
547 
548 
552 
553 //=============================================================================
561 inline long zgematrix::zgesvd(dcovector& S, zgematrix& U, zgematrix& VT)
562 {
563 #ifdef CPPL_VERBOSE
564  std::cerr << "# [MARK] zgematrix::zgesvd(dcovector&, zgematrix&, zgematrix&)"
565  << std::endl;
566 #endif//CPPL_VERBOSE
567 
568  char JOBU('A'), JOBVT('A');
569  long LDA(M), LDU(M), LDVT(N),
570  LWORK(max(3*min(m,n)+max(m,n),5*min(m,n))), INFO(1);
571  double *RWORK(new double[5*min(m,n)]);
572  std::complex<double> *WORK(new std::complex<double>[LWORK]);
573  S.resize(min(M,N)); U.resize(LDU,M); VT.resize(LDVT,N);
574 
575  zgesvd_(JOBU, JOBVT, M, N, Array, LDA, S.array, U.Array,
576  LDU, VT.Array, LDVT, WORK, LWORK, RWORK, INFO);
577  delete [] RWORK; delete [] WORK;
578 
579  if(INFO!=0){
580  std::cerr << "[WARNING] zgematrix::zgesvd"
581  << "(dceovector&, zgematrix&, zcovector&) "
582  << "Serious trouble happend. INFO = " << INFO << "."
583  << std::endl;
584  }
585  return INFO;
586 }
long zgesvd(dcovector &, zgematrix &, zgematrix &)
Definition: zgematrix-lapack.hpp:561
void resize(const long &, const long &)
Definition: zgematrix-misc.hpp:126
void clear()
Definition: zcovector-misc.hpp:3
void resize(const long &)
Definition: dcovector-misc.hpp:93
double *const & array
1D array to store vector data (readable)
Definition: dcovector.hpp:13
long zgelss(zcovector &, dcovector &, long &, const double)
Definition: zgematrix-lapack.hpp:297
friend _zgematrix i(const zgematrix &)
Definition: zgematrix-calc.hpp:21
void zero()
Definition: drovector-misc.hpp:23
Complex Double-precision General Dence Matrix Class.
Definition: zgematrix.hpp:3
long const & l
vector size (readable)
Definition: drovector.hpp:12
Real Double-precision Row Vector Class.
Definition: drovector.hpp:3
long zgels(zgematrix &)
Definition: zgematrix-lapack.hpp:72
void clear()
Definition: zgematrix-misc.hpp:3
long zgesv(zgematrix &)
Definition: zgematrix-lapack.hpp:5
long const & m
matrix row size (readable)
Definition: zgematrix.hpp:14
friend void swap(zgematrix &, zgematrix &)
Definition: zgematrix-misc.hpp:154
Real Double-precision Column Vector Class.
Definition: dcovector.hpp:3
long const & n
matrix column size (readable)
Definition: zgematrix.hpp:15
Complex Double-precision Column Vector Class.
Definition: zcovector.hpp:3
void resize(const long &)
Definition: drovector-misc.hpp:93
std::complex< double > * Array
1D Array to store vector data
Definition: _zrovector.hpp:8
long zgeev(std::vector< std::complex< double > > &)
Definition: zgematrix-lapack.hpp:403