00001
00005 inline long zgematrix::zgesv(zgematrix& mat)
00006 {
00007 #ifdef CPPL_VERBOSE
00008 std::cerr << "# [MARK] zgematrix::zgesv(zgematrix&)"
00009 << std::endl;
00010 #endif//CPPL_VERBOSE
00011
00012 #ifdef CPPL_DEBUG
00013 if(M!=N || M!=mat.M){
00014 std::cerr << "[ERROR] zgematrix::zgesv(zgematrix&) " << std::endl
00015 << "These two matrices cannot be solved." << std::endl
00016 << "Your input was (" << M << "x" << N << ") and ("
00017 << mat.M << "x" << mat.N << ")." << std::endl;
00018 exit(1);
00019 }
00020 #endif//CPPL_DEBUG
00021
00022 long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.M), INFO(1);
00023 zgesv_(N, NRHS, Array, LDA, IPIV, mat.Array, LDB, INFO);
00024 delete [] IPIV;
00025
00026 if(INFO!=0){
00027 std::cerr << "[WARNING] zgematrix::zgesv(zgematrix&) "
00028 << "Serious trouble happend. INFO = " << INFO << "."
00029 << std::endl;
00030 }
00031 return INFO;
00032 }
00033
00034
00038 inline long zgematrix::zgesv(zcovector& vec)
00039 {
00040 #ifdef CPPL_VERBOSE
00041 std::cerr << "# [MARK] zgematrix::zgesv(zcovector&)"
00042 << std::endl;
00043 #endif//CPPL_VERBOSE
00044
00045 #ifdef CPPL_DEBUG
00046 if(M!=N || M!=vec.L){
00047 std::cerr << "[ERROR] zgematrix::zgesv(zcovector&) " << std::endl
00048 << "These matrix and vector cannot be solved." << std::endl
00049 << "Your input was (" << M << "x" << N << ") and ("
00050 << vec.L << ")." << std::endl;
00051 exit(1);
00052 }
00053 #endif//CPPL_DEBUG
00054 long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), INFO(1);
00055 zgesv_(N, NRHS, Array, LDA, IPIV, vec.Array, LDB, INFO);
00056 delete [] IPIV;
00057
00058 if(INFO!=0){
00059 std::cerr << "[WARNING] zgematrix::zgesv(zcovector&) "
00060 << "Serious trouble happend. INFO = " << INFO << "."
00061 << std::endl;
00062 }
00063 return INFO;
00064 }
00065
00069
00070
00072 inline long zgematrix::zgels(zgematrix& mat)
00073 {
00074 #ifdef CPPL_VERBOSE
00075 std::cerr << "# [MARK] zgematrix::zgels(zgematrix&)"
00076 << std::endl;
00077 #endif//CPPL_VERBOSE
00078
00079 #ifdef CPPL_DEBUG
00080 if(M!=mat.M){
00081 std::cerr << "[ERROR] zgematrix::zgels(zgematrix&) " << std::endl
00082 << "These two matrices cannot be solved." << std::endl
00083 << "Your input was (" << M << "x" << N << ") and ("
00084 << mat.M << "x" << mat.N << ")." << std::endl;
00085 exit(1);
00086 }
00087 #endif//CPPL_DEBUG
00088
00089 if(M<N){
00090 zgematrix tmp(N,mat.N);
00091 for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
00092 tmp(i,j) =mat(i,j);
00093 }}
00094 mat.clear();
00095 swap(mat,tmp);
00096 }
00097
00098 char TRANS('N');
00099 long NRHS(mat.N), LDA(M), LDB(mat.M),
00100 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00101 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00102 zgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
00103 delete [] WORK;
00104
00105 if(M>N){
00106 zgematrix tmp(N,mat.N);
00107 for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
00108 tmp(i,j) =mat(i,j);
00109 }}
00110 mat.clear();
00111 swap(mat,tmp);
00112 }
00113
00114 if(INFO!=0){
00115 std::cerr << "[WARNING] zgematrix::zgels(zgematrix&) "
00116 << "Serious trouble happend. INFO = " << INFO << "."
00117 << std::endl;
00118 }
00119 return INFO;
00120 }
00121
00122
00124 inline long zgematrix::zgels(zcovector& vec)
00125 {
00126 #ifdef CPPL_VERBOSE
00127 std::cerr << "# [MARK] zgematrix::zgels(zcovector&)"
00128 << std::endl;
00129 #endif//CPPL_VERBOSE
00130
00131 #ifdef CPPL_DEBUG
00132 if(M!=vec.L){
00133 std::cerr << "[ERROR] zgematrix::zgels(zcovector&) " << std::endl
00134 << "These matrix and vector cannot be solved." << std::endl
00135 << "Your input was (" << M << "x" << N << ") and ("
00136 << vec.L << ")." << std::endl;
00137 exit(1);
00138 }
00139 #endif//CPPL_DEBUG
00140
00141 if(M<N){
00142 zcovector tmp(N);
00143 for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
00144 vec.clear();
00145 swap(vec,tmp);
00146 }
00147
00148 char TRANS('N');
00149 long NRHS(1), LDA(M), LDB(vec.L),
00150 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00151 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00152 zgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
00153 delete [] WORK;
00154
00155 if(M>N){
00156 zcovector tmp(N);
00157 for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
00158 vec.clear();
00159 swap(vec,tmp);
00160 }
00161
00162 if(INFO!=0){
00163 std::cerr << "[WARNING] zgematrix::zgels(zcovector&) "
00164 << "Serious trouble happend. INFO = " << INFO << "."
00165 << std::endl;
00166 }
00167 return INFO;
00168 }
00169
00170
00177 inline long zgematrix::zgels(zgematrix& mat, drovector& residual)
00178 {
00179 #ifdef CPPL_VERBOSE
00180 std::cerr << "# [MARK] zgematrix::zgels(zgematrix&, drovector&)"
00181 << std::endl;
00182 #endif//CPPL_VERBOSE
00183
00184 #ifdef CPPL_DEBUG
00185 if(M!=mat.M){
00186 std::cerr << "[ERROR] zgematrix::zgels(zgematrix&, drovector&) "
00187 << std::endl
00188 << "These two matrices cannot be solved." << std::endl
00189 << "Your input was (" << M << "x" << N << ") and ("
00190 << mat.M << "x" << mat.N << ")." << std::endl;
00191 exit(1);
00192 }
00193 #endif//CPPL_DEBUG
00194
00195 residual.resize(mat.N); residual.zero();
00196
00197 if(M<N){
00198 zgematrix tmp(N,mat.N);
00199 for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
00200 tmp(i,j) =mat(i,j);
00201 }}
00202 mat.clear();
00203 swap(mat,tmp);
00204 }
00205
00206 char TRANS('N');
00207 long NRHS(mat.N), LDA(M), LDB(mat.M),
00208 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00209 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00210 zgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
00211 delete [] WORK;
00212
00213 if(M>N){
00214 for(long i=0; i<residual.l; i++){ for(long j=0; j<M-N; j++){
00215 residual(i) += std::norm(mat(N+j,i));
00216 }}
00217
00218 zgematrix tmp(N,mat.N);
00219 for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
00220 tmp(i,j) =mat(i,j);
00221 }}
00222 mat.clear();
00223 swap(mat,tmp);
00224 }
00225
00226 if(INFO!=0){
00227 std::cerr << "[WARNING] zgematrix::zgels(zgematrix&, drovector&) "
00228 << "Serious trouble happend. INFO = " << INFO << "."
00229 << std::endl;
00230 }
00231 return INFO;
00232 }
00233
00234
00240 inline long zgematrix::zgels(zcovector& vec, double& residual)
00241 {
00242 #ifdef CPPL_VERBOSE
00243 std::cerr << "# [MARK] zgematrix::zgels(zcovector&, double&)"
00244 << std::endl;
00245 #endif//CPPL_VERBOSE
00246
00247 #ifdef CPPL_DEBUG
00248 if(M!=vec.L){
00249 std::cerr << "[ERROR] zgematrix::zgels(zcovector&, double&) " << std::endl
00250 << "These matrix and vector cannot be solved." << std::endl
00251 << "Your input was (" << M << "x" << N << ") and ("
00252 << vec.L << ")." << std::endl;
00253 exit(1);
00254 }
00255 #endif//CPPL_DEBUG
00256
00257 residual=0.0;
00258
00259 if(M<N){
00260 zcovector tmp(n);
00261 for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
00262 vec.clear();
00263 swap(vec,tmp);
00264 }
00265
00266 char TRANS('N');
00267 long NRHS(1), LDA(m), LDB(vec.L),
00268 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00269 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00270 zgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
00271 delete [] WORK;
00272
00273 if(M>N){
00274 for(long i=0; i<M-N; i++){ residual+=std::norm(vec(N+i)); }
00275
00276 zcovector tmp(N);
00277 for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
00278 vec.clear();
00279 swap(vec,tmp);
00280 }
00281
00282 if(INFO!=0){
00283 std::cerr << "[WARNING] zgematrix::zgels(zcovector&, double&) "
00284 << "Serious trouble happend. INFO = " << INFO << "."
00285 << std::endl;
00286 }
00287 return INFO;
00288 }
00289
00293
00294
00297 inline long zgematrix::zgelss(zcovector& B, dcovector& S, long& RANK,
00298 const double RCOND =-1. )
00299 {
00300 #ifdef CPPL_VERBOSE
00301 std::cerr << "# [MARK] zgematrix::zgelss(zcovector&, dcovector&, long&, const double)"
00302 << std::endl;
00303 #endif//CPPL_VERBOSE
00304
00305 #ifdef CPPL_DEBUG
00306 if(M!=B.L){
00307 std::cerr << "[ERROR] zgematrix::zgelss"
00308 << "(zcovector&, dcovector&, long&, const double) " << std::endl
00309 << "These matrix and vector cannot be solved." << std::endl
00310 << "Your input was (" << M << "x" << N << ") and ("
00311 << B.L << ")." << std::endl;
00312 exit(1);
00313 }
00314 #endif//CPPL_DEBUG
00315
00316 if(M<N){
00317 zcovector tmp(N);
00318 for(long i=0; i<B.L; i++){ tmp(i)=B(i); }
00319 B.clear();
00320 swap(B,tmp);
00321 }
00322
00323 S.resize(min(M,N));
00324
00325 long NRHS(1), LDA(M), LDB(B.L),
00326 LWORK(2*min(M,N) +max(max(M,N),NRHS)), INFO(1);
00327 double *RWORK(new double[5*min(M,N)]);
00328 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00329 zgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.array, RCOND, RANK,
00330 WORK, LWORK, RWORK, INFO);
00331 delete [] RWORK; delete [] WORK;
00332
00333 if(INFO!=0){
00334 std::cerr << "[WARNING] zgematrix::zgelss"
00335 << "(zcovector&, dcovector&, long&, double) "
00336 << "Serious trouble happend. INFO = " << INFO << "."
00337 << std::endl;
00338 }
00339 return INFO;
00340 }
00341
00342
00345 inline long zgematrix::zgelss(zgematrix& B, dcovector& S, long& RANK,
00346 const double RCOND =-1. )
00347 {
00348 #ifdef CPPL_VERBOSE
00349 std::cerr << "# [MARK] zgematrix::zgelss(zgematrix&, dcovector&, long&, const double)"
00350 << std::endl;
00351 #endif//CPPL_VERBOSE
00352
00353 #ifdef CPPL_DEBUG
00354 if(M!=B.M){
00355 std::cerr << "[ERROR] zgematrix::zgelss"
00356 << "(zgematrix&, dcovector&, long&, const double) " << std::endl
00357 << "These matrix and vector cannot be solved." << std::endl
00358 << "Your input was (" << M << "x" << N << ") and ("
00359 << B.M << "x" << B.N << ")." << std::endl;
00360 exit(1);
00361 }
00362 #endif//CPPL_DEBUG
00363
00364 if(M<N){
00365 zgematrix tmp(N,B.N);
00366 for(long i=0; i<B.M; i++){
00367 for(long j=0; j<B.N; j++){
00368 tmp(i,j)=B(i,j);
00369 }
00370 }
00371 B.clear();
00372 swap(B,tmp);
00373 }
00374
00375 S.resize(min(M,N));
00376
00377 long NRHS(B.N), LDA(M), LDB(B.M),
00378 LWORK(2*min(M,N) +max(max(M,N),NRHS)), INFO(1);
00379 double *RWORK(new double[5*min(M,N)]);
00380 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00381 zgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.array, RCOND, RANK,
00382 WORK, LWORK, RWORK, INFO);
00383 delete [] RWORK; delete [] WORK;
00384
00385 if(INFO!=0){
00386 std::cerr << "[WARNING] zgematrix::zgelss"
00387 << "(zgematrix&, dcovector&, long&, const double) "
00388 << "Serious trouble happend. INFO = " << INFO << "."
00389 << std::endl;
00390 }
00391 return INFO;
00392 }
00393
00397
00398
00403 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w)
00404 {
00405 #ifdef CPPL_VERBOSE
00406 std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&)"
00407 << std::endl;
00408 #endif//CPPL_VERBOSE
00409
00410 #ifdef CPPL_DEBUG
00411 if(M!=N){
00412 std::cerr << "[ERROR] zgematrix::zgeev"
00413 << "(std::vector<std::complex<double>>&) "
00414 << std::endl
00415 << "This matrix cannot have eigenvalues." << std::endl
00416 << "Your input was (" << M << "x" << N << ")." << std::endl;
00417 exit(1);
00418 }
00419 #endif//CPPL_DEBUG
00420
00421 w.resize(N);
00422 char JOBVL('N'), JOBVR('N');
00423 long LDA(N), LDVL(1), LDVR(1), LWORK(4*n), INFO(1);
00424 double *RWORK(new double[2*n]);
00425 std::complex<double> *VL(NULL), *VR(NULL),
00426 *WORK(new std::complex<double>[LWORK]);
00427 zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
00428 VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO);
00429 delete [] RWORK; delete [] WORK; delete [] VL; delete [] VR;
00430
00431 if(INFO!=0){
00432 std::cerr << "[WARNING] zgematrix::zgesv"
00433 << "(std::vector<std::complex<double>&, std::vector<zcovector>&) "
00434 << "Serious trouble happend. INFO = " << INFO << "."
00435 << std::endl;
00436 }
00437 return INFO;
00438 }
00439
00440
00446 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w,
00447 std::vector<zcovector>& vr)
00448 {
00449 #ifdef CPPL_VERBOSE
00450 std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&, std::vector<zcovector>&)"
00451 << std::endl;
00452 #endif//CPPL_VERBOSE
00453
00454 #ifdef CPPL_DEBUG
00455 if(M!=N){
00456 std::cerr << "[ERROR] zgematrix::zgeev"
00457 << "(std::vector<std::complex<double>>&, std::vector<zcovector>&)"
00458 << std::endl
00459 << "This matrix cannot have eigenvalues." << std::endl
00460 << "Your input was (" << M << "x" << N << ")." << std::endl;
00461 exit(1);
00462 }
00463 #endif//CPPL_DEBUG
00464
00465 w.resize(N); vr.resize(N);
00466 for(long i=0; i<N; i++){ vr[i].resize(N); }
00467 zgematrix VR(N,N);
00468 char JOBVL('N'), JOBVR('V');
00469 long LDA(N), LDVL(1), LDVR(N), LWORK(4*n), INFO(1);
00470 double *RWORK(new double[2*n]);
00471 std::complex<double> *VL(NULL), *WORK(new std::complex<double>[LWORK]);
00472 zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
00473 VL, LDVL, VR.Array, LDVR, WORK, LWORK, RWORK, INFO);
00474 delete [] RWORK; delete [] WORK; delete [] VL;
00475
00477 for(long j=0; j<N; j++){ for(long i=0; i<N; i++){
00478 vr[j](i) = VR(i,j);
00479 }}
00480
00481 if(INFO!=0){
00482 std::cerr << "[WARNING] zgematrix::zgesv"
00483 << "(std::vector<std::complex<double>&, std::vector<zcovector>&) "
00484 << "Serious trouble happend. INFO = " << INFO << "."
00485 << std::endl;
00486 }
00487 return INFO;
00488 }
00489
00490
00496 inline long zgematrix::zgeev(std::vector< std::complex<double> >& w,
00497 std::vector<zrovector>& vl)
00498 {
00499 #ifdef CPPL_VERBOSE
00500 std::cerr << "# [MARK] zgematrix::zgeev(std::vector< std::complex<double> >&, std::vector<zrovector>&)"
00501 << std::endl;
00502 #endif//CPPL_VERBOSE
00503
00504 #ifdef CPPL_DEBUG
00505 if(M!=N){
00506 std::cerr << "[ERROR] zgematrix::zgeev"
00507 << "(std::vector<std::complex<double>>&, std::vector<zrovector>&)"
00508 << std::endl
00509 << "This matrix cannot have eigenvalues." << std::endl
00510 << "Your input was (" << M << "x" << N << ")." << std::endl;
00511 exit(1);
00512 }
00513 #endif//CPPL_DEBUG
00514
00515 w.resize(N); vl.resize(N);
00516 for(long i=0; i<N; i++){ vl[i].resize(N); }
00517 zgematrix VL(N,N);
00518 char JOBVL('V'), JOBVR('N');
00519 long LDA(N), LDVL(N), LDVR(1), LWORK(4*n), INFO(1);
00520 double *RWORK(new double[2*n]);
00521 std::complex<double> *VR(NULL), *WORK(new std::complex<double>[LWORK]);
00522 zgeev_(JOBVL, JOBVR, N, Array, LDA, &w[0],
00523 VL.Array, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO);
00524 delete [] RWORK; delete [] WORK; delete [] VR;
00525
00527 for(long j=0; j<N; j++){ for(long i=0; i<N; i++){
00528 vl[j](i) = std::conj(VL(i,j));
00529 }}
00530
00531 if(INFO!=0){
00532 std::cerr << "[WARNING] zgematrix::zgesv"
00533 << "(std::vector<std::complex<double>&, std::vector<zrovector>&) "
00534 << "Serious trouble happend. INFO = " << INFO << "."
00535 << std::endl;
00536 }
00537 return INFO;
00538 }
00539
00540
00544
00545
00546
00547
00548
00552
00553
00561 inline long zgematrix::zgesvd(dcovector& S, zgematrix& U, zgematrix& VT)
00562 {
00563 #ifdef CPPL_VERBOSE
00564 std::cerr << "# [MARK] zgematrix::zgesvd(dcovector&, zgematrix&, zgematrix&)"
00565 << std::endl;
00566 #endif//CPPL_VERBOSE
00567
00568 char JOBU('A'), JOBVT('A');
00569 long LDA(M), LDU(M), LDVT(N),
00570 LWORK(max(3*min(m,n)+max(m,n),5*min(m,n))), INFO(1);
00571 double *RWORK(new double[5*min(m,n)]);
00572 std::complex<double> *WORK(new std::complex<double>[LWORK]);
00573 S.resize(min(M,N)); U.resize(LDU,M); VT.resize(LDVT,N);
00574
00575 zgesvd_(JOBU, JOBVT, M, N, Array, LDA, S.array, U.Array,
00576 LDU, VT.Array, LDVT, WORK, LWORK, RWORK, INFO);
00577 delete [] RWORK; delete [] WORK;
00578
00579 if(INFO!=0){
00580 std::cerr << "[WARNING] zgematrix::zgesvd"
00581 << "(dceovector&, zgematrix&, zcovector&) "
00582 << "Serious trouble happend. INFO = " << INFO << "."
00583 << std::endl;
00584 }
00585 return INFO;
00586 }