00001
00006 inline long dgematrix::dgesv(dgematrix& mat)
00007 {
00008 #ifdef CPPL_VERBOSE
00009 std::cerr << "# [MARK] dgematrix::dgesv(dgematrix&)"
00010 << std::endl;
00011 #endif//CPPL_VERBOSE
00012
00013 #ifdef CPPL_DEBUG
00014 if(M!=N || M!=mat.M){
00015 std::cerr << "[ERROR] dgematrix::dgesv(dgematrix&) " << std::endl
00016 << "These two matrices cannot be solved." << std::endl
00017 << "Your input was (" << M << "x" << N << ") and ("
00018 << mat.M << "x" << mat.N << ")." << std::endl;
00019 exit(1);
00020 }
00021 #endif//CPPL_DEBUG
00022
00023 long NRHS(mat.N), LDA(N), *IPIV(new long[N]), LDB(mat.M), INFO(1);
00024 dgesv_(N, NRHS, Array, LDA, IPIV, mat.Array, LDB, INFO);
00025 delete [] IPIV;
00026
00027 if(INFO!=0){
00028 std::cerr << "[WARNING] dgematrix::dgesv(dgematrix&) "
00029 << "Serious trouble happend. INFO = " << INFO << "."
00030 << std::endl;
00031 }
00032 return INFO;
00033 }
00034
00035
00040 inline long dgematrix::dgesv(dcovector& vec)
00041 {
00042 #ifdef CPPL_VERBOSE
00043 std::cerr << "# [MARK] dgematrix::dgesv(dcovector&)"
00044 << std::endl;
00045 #endif//CPPL_VERBOSE
00046
00047 #ifdef CPPL_DEBUG
00048 if(M!=N || M!=vec.L){
00049 std::cerr << "[ERROR] dgematrix::dgesv(dcovector&) " << std::endl
00050 << "These matrix and vector cannot be solved." << std::endl
00051 << "Your input was (" << M << "x" << N << ") and ("
00052 << vec.L << ")." << std::endl;
00053 exit(1);
00054 }
00055 #endif//CPPL_DEBUG
00056
00057 long NRHS(1), LDA(N), *IPIV(new long[N]), LDB(vec.L), INFO(1);
00058 dgesv_(N, NRHS, Array, LDA, IPIV, vec.Array, LDB, INFO);
00059 delete [] IPIV;
00060
00061 if(INFO!=0){
00062 std::cerr << "[WARNING] dgematrix::dgesv(dcovector&) "
00063 << "Serious trouble happend. INFO = " << INFO << "."
00064 << std::endl;
00065 }
00066 return INFO;
00067 }
00068
00072
00073
00076 inline long dgematrix::dgels(dgematrix& mat)
00077 {
00078 #ifdef CPPL_VERBOSE
00079 std::cerr << "# [MARK] dgematrix::dgels(dgematrix&)"
00080 << std::endl;
00081 #endif//CPPL_VERBOSE
00082
00083 #ifdef CPPL_DEBUG
00084 if(M!=mat.M){
00085 std::cerr << "[ERROR] dgematrix::dgels(dgematrix&) " << std::endl
00086 << "These two matrices cannot be solved." << std::endl
00087 << "Your input was (" << M << "x" << N << ") and ("
00088 << mat.M << "x" << mat.N << ")." << std::endl;
00089 exit(1);
00090 }
00091 #endif//CPPL_DEBUG
00092
00093 if(M<N){
00094 dgematrix tmp(N,mat.N);
00095 for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
00096 tmp(i,j) =mat(i,j);
00097 }}
00098 mat.clear();
00099 swap(mat,tmp);
00100 }
00101
00102 char TRANS('N');
00103 long NRHS(mat.N), LDA(M), LDB(mat.M),
00104 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00105 double *WORK(new double[LWORK]);
00106 dgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
00107 delete [] WORK;
00108
00109 if(M>N){
00110 dgematrix tmp(N,mat.N);
00111 for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
00112 tmp(i,j) =mat(i,j);
00113 }}
00114 mat.clear();
00115 swap(mat,tmp);
00116 }
00117
00118 if(INFO!=0){
00119 std::cerr << "[WARNING] dgematrix::dgels(dgematrix&) "
00120 << "Serious trouble happend. INFO = " << INFO << "."
00121 << std::endl;
00122 }
00123 return INFO;
00124 }
00125
00126
00129 inline long dgematrix::dgels(dcovector& vec)
00130 {
00131 #ifdef CPPL_VERBOSE
00132 std::cerr << "# [MARK] dgematrix::dgels(dcovector&)"
00133 << std::endl;
00134 #endif//CPPL_VERBOSE
00135
00136 #ifdef CPPL_DEBUG
00137 if(M!=vec.L){
00138 std::cerr << "[ERROR] dgematrix::dgels(dcovector&) " << std::endl
00139 << "These matrix and vector cannot be solved." << std::endl
00140 << "Your input was (" << M << "x" << N << ") and ("
00141 << vec.L << ")." << std::endl;
00142 exit(1);
00143 }
00144 #endif//CPPL_DEBUG
00145
00146 if(M<N){
00147 dcovector tmp(N);
00148 for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
00149 vec.clear();
00150 swap(vec,tmp);
00151 }
00152
00153 char TRANS('N');
00154 long NRHS(1), LDA(M), LDB(vec.L),
00155 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00156 double *WORK(new double[LWORK]);
00157 dgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
00158 delete [] WORK;
00159
00160 if(M>N){
00161 dcovector tmp(N);
00162 for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
00163 vec.clear();
00164 swap(vec,tmp);
00165 }
00166
00167 if(INFO!=0){
00168 std::cerr << "[WARNING] dgematrix::dgels(dcovector&) "
00169 << "Serious trouble happend. INFO = " << INFO << "."
00170 << std::endl;
00171 }
00172 return INFO;
00173 }
00174
00175
00182 inline long dgematrix::dgels(dgematrix& mat, drovector& residual)
00183 {
00184 #ifdef CPPL_VERBOSE
00185 std::cerr << "# [MARK] dgematrix::dgels(dgematrix&, drovector&)"
00186 << std::endl;
00187 #endif//CPPL_VERBOSE
00188
00189 #ifdef CPPL_DEBUG
00190 if(M!=mat.M){
00191 std::cerr << "[ERROR] dgematrix::dgels(dgematrix&, drovector&) "
00192 << std::endl
00193 << "These two matrices cannot be solved." << std::endl
00194 << "Your input was (" << M << "x" << N << ") and ("
00195 << mat.M << "x" << mat.N << ")." << std::endl;
00196 exit(1);
00197 }
00198 #endif//CPPL_DEBUG
00199
00200 residual.resize(mat.N); residual.zero();
00201
00202 if(M<N){
00203 dgematrix tmp(N,mat.N);
00204 for(long i=0; i<mat.M; i++){ for(long j=0; j<mat.N; j++){
00205 tmp(i,j) =mat(i,j);
00206 }}
00207 mat.clear();
00208 swap(mat,tmp);
00209 }
00210
00211 char TRANS('N');
00212 long NRHS(mat.N), LDA(M), LDB(mat.M),
00213 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00214 double *WORK(new double[LWORK]);
00215 dgels_(TRANS, M, N, NRHS, Array, LDA, mat.Array, LDB, WORK, LWORK, INFO);
00216 delete [] WORK;
00217
00218 if(M>N){
00219 for(long i=0; i<residual.L; i++){ for(long j=0; j<M-N; j++){
00220 residual(i) += std::pow(mat(N+j,i), 2.0);
00221 }}
00222
00223 dgematrix tmp(N,mat.N);
00224 for(long i=0; i<tmp.M; i++){ for(long j=0; j<tmp.N; j++){
00225 tmp(i,j) =mat(i,j);
00226 }}
00227 mat.clear();
00228 swap(mat,tmp);
00229 }
00230
00231 if(INFO!=0){
00232 std::cerr << "[WARNING] dgematrix::dgels(dgematrix&, drovector&) "
00233 << "Serious trouble happend. INFO = " << INFO << "."
00234 << std::endl;
00235 }
00236 return INFO;
00237 }
00238
00239
00246 inline long dgematrix::dgels(dcovector& vec, double& residual)
00247 {
00248 #ifdef CPPL_VERBOSE
00249 std::cerr << "# [MARK] dgematrix::dgels(dcovector&, double&)"
00250 << std::endl;
00251 #endif//CPPL_VERBOSE
00252
00253 #ifdef CPPL_DEBUG
00254 if(M!=vec.L){
00255 std::cerr << "[ERROR] dgematrix::dgels(dcovector&, double&) " << std::endl
00256 << "These matrix and vector cannot be solved." << std::endl
00257 << "Your input was (" << M << "x" << N << ") and ("
00258 << vec.L << ")." << std::endl;
00259 exit(1);
00260 }
00261 #endif//CPPL_DEBUG
00262
00263 residual=0.0;
00264
00265 if(M<N){
00266 dcovector tmp(N);
00267 for(long i=0; i<vec.L; i++){ tmp(i)=vec(i); }
00268 vec.clear();
00269 swap(vec,tmp);
00270 }
00271
00272 char TRANS('N');
00273 long NRHS(1), LDA(M), LDB(vec.L),
00274 LWORK(min(M,N)+max(min(M,N),NRHS)), INFO(1);
00275 double *WORK(new double[LWORK]);
00276 dgels_(TRANS, M, N, NRHS, Array, LDA, vec.Array, LDB, WORK, LWORK, INFO);
00277 delete [] WORK;
00278
00279 if(M>N){
00280 for(long i=0; i<M-N; i++){ residual+=std::pow(vec(N+i),2.0); }
00281
00282 dcovector tmp(N);
00283 for(long i=0; i<tmp.L; i++){ tmp(i)=vec(i); }
00284 vec.clear();
00285 swap(vec,tmp);
00286 }
00287
00288 if(INFO!=0){
00289 std::cerr << "[WARNING] dgematrix::dgels(dcovector&, double&) "
00290 << "Serious trouble happend. INFO = " << INFO << "."
00291 << std::endl;
00292 }
00293 return INFO;
00294 }
00295
00299
00300
00304 inline long dgematrix::dgelss(dcovector& B, dcovector& S, long& RANK,
00305 const double RCOND =-1. )
00306 {
00307 #ifdef CPPL_VERBOSE
00308 std::cerr << "# [MARK] dgematrix::dgelss(dcovector&, dcovector&, long&, const double)"
00309 << std::endl;
00310 #endif//CPPL_VERBOSE
00311
00312 #ifdef CPPL_DEBUG
00313 if(M!=B.L){
00314 std::cerr << "[ERROR] dgematrix::dgelss"
00315 << "(dcovector&, dcovector&, long, double) " << std::endl
00316 << "These matrix and vector cannot be solved." << std::endl
00317 << "Your input was (" << M << "x" << N << ") and ("
00318 << B.L << ")." << std::endl;
00319 exit(1);
00320 }
00321 #endif//CPPL_DEBUG
00322
00323 if(M<N){
00324 dcovector tmp(N);
00325 for(long i=0; i<B.L; i++){ tmp(i)=B(i); }
00326 B.clear();
00327 swap(B,tmp);
00328 }
00329
00330 S.resize(min(M,N));
00331
00332 long NRHS(1), LDA(M), LDB(B.L),
00333 LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
00334 double *WORK(new double[LWORK]);
00335 dgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.Array,
00336 RCOND, RANK, WORK, LWORK, INFO);
00337 delete [] WORK;
00338
00339 if(INFO!=0){
00340 std::cerr << "[WARNING] dgematrix::dgelss"
00341 << "(dcovector&, docvector&, long, const double) "
00342 << "Serious trouble happend. INFO = " << INFO << "."
00343 << std::endl;
00344 }
00345 return INFO;
00346 }
00347
00348
00352 inline long dgematrix::dgelss(dgematrix& B, dcovector& S, long& RANK,
00353 const double RCOND =-1. )
00354 {
00355 #ifdef CPPL_VERBOSE
00356 std::cerr << "# [MARK] dgematrix::dgelss(dgematrix&, dcovector&, long& const double)"
00357 << std::endl;
00358 #endif//CPPL_VERBOSE
00359
00360 #ifdef CPPL_DEBUG
00361 if(M!=B.M){
00362 std::cerr << "[ERROR] dgematrix::dgelss"
00363 << "(dgematrix&, dcovector&, long&, const double) " << std::endl
00364 << "These matrix and vector cannot be solved." << std::endl
00365 << "Your input was (" << M << "x" << N << ") and ("
00366 << B.M << "x" << B.N << ")." << std::endl;
00367 exit(1);
00368 }
00369 #endif//CPPL_DEBUG
00370
00371 if(M<N){
00372 dgematrix tmp(N,B.N);
00373 for(long i=0; i<B.M; i++){
00374 for(long j=0; j<B.N; j++){
00375 tmp(i,j)=B(i,j);
00376 }
00377 }
00378 B.clear();
00379 swap(B,tmp);
00380 }
00381
00382 S.resize(min(M,N));
00383
00384 long NRHS(B.N), LDA(M), LDB(B.M),
00385 LWORK(3*min(M,N)+max(max(2*min(M,N),max(M,N)), NRHS)), INFO(1);
00386 double *WORK(new double[LWORK]);
00387 dgelss_(M, N, NRHS, Array, LDA, B.Array, LDB, S.Array,
00388 RCOND, RANK, WORK, LWORK, INFO);
00389 delete [] WORK;
00390
00391 if(INFO!=0){
00392 std::cerr << "[WARNING] dgematrix::dgelss"
00393 << "(dgematrix&, docvector&, long, const double) "
00394 << "Serious trouble happend. INFO = " << INFO << "."
00395 << std::endl;
00396 }
00397 return INFO;
00398 }
00399
00403
00410 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi)
00411 {
00412 #ifdef CPPL_VERBOSE
00413 std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&)"
00414 << std::endl;
00415 #endif//CPPL_VERBOSE
00416
00417 #ifdef CPPL_DEBUG
00418 if(M!=N){
00419 std::cerr << "[ERROR] dgematrix::dgeev"
00420 << "(vector<double>&, vector<double>&) " << std::endl
00421 << "This matrix is not a square matrix." << std::endl
00422 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00423 exit(1);
00424 }
00425 #endif//CPPL_DEBUG
00426
00427 wr.resize(N); wi.resize(N);
00428 char JOBVL('N'), JOBVR('N');
00429 long LDA(N), LDVL(1), LDVR(1), LWORK(3*N), INFO(1);
00430 double *VL(NULL), *VR(NULL), *WORK(new double[LWORK]);
00431 dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
00432 VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
00433 delete [] WORK; delete [] VL; delete [] VL;
00434
00435 if(INFO!=0){
00436 std::cerr << "[WARNING] dgematrix::dgeev"
00437 << "(vector<double>&, vector<double>&)"
00438 << std::endl
00439 << "Serious trouble happend. INFO = " << INFO << "."
00440 << std::endl;
00441 }
00442 return INFO;
00443 }
00444
00445
00453 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
00454 std::vector<dcovector>& vrr,
00455 std::vector<dcovector>& vri)
00456 {
00457 #ifdef CPPL_VERBOSE
00458 std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
00459 << std::endl;
00460 #endif//CPPL_VERBOSE
00461
00462 #ifdef CPPL_DEBUG
00463 if(M!=N){
00464 std::cerr << "[ERROR] dgematrix::dgeev"
00465 << "(vector<double>&, vector<double>&, "
00466 << "vector<dcovector>&, vector<dcovector>&) "
00467 << std::endl
00468 << "This matrix is not a square matrix." << std::endl
00469 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00470 exit(1);
00471 }
00472 #endif//CPPL_DEBUG
00473
00474 wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
00475 for(long i=0; i<N; i++){ vrr[i].resize(N); vri[i].resize(N); }
00476 dgematrix VR(N,N);
00477 char JOBVL('N'), JOBVR('V');
00478 long LDA(N), LDVL(1), LDVR(N), LWORK(4*N), INFO(1);
00479 double *VL(NULL), *WORK(new double[LWORK]);
00480 dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
00481 VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
00482 delete [] WORK; delete [] VL;
00483
00485 for(long j=0; j<N; j++){
00486 if(fabs(wi[j])<1e-10){
00487 for(long i=0; i<N; i++){
00488 vrr[j](i) = VR(i,j); vri[j](i) = 0.0;
00489 }
00490 }
00491 else{
00492 for(long i=0; i<N; i++){
00493 vrr[j](i) = VR(i,j); vri[j](i) = VR(i,j+1);
00494 vrr[j+1](i) = VR(i,j); vri[j+1](i) =-VR(i,j+1);
00495 }
00496 j++;
00497 }
00498 }
00499
00500 if(INFO!=0){
00501 std::cerr << "[WARNING] dgematrix::dgeev"
00502 << "(vector<double>&, vector<double>&, "
00503 << "vector<dcovector>&, vector<dcovector>&) "
00504 << std::endl
00505 << "Serious trouble happend. INFO = " << INFO << "."
00506 << std::endl;
00507 }
00508 return INFO;
00509 }
00510
00511
00519 inline long dgematrix::dgeev(std::vector<double>& wr, std::vector<double>& wi,
00520 std::vector<drovector>& vlr,
00521 std::vector<drovector>& vli)
00522 {
00523 #ifdef CPPL_VERBOSE
00524 std::cerr << "# [MARK] dgematrix::dgeev(std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
00525 << std::endl;
00526 #endif//CPPL_VERBOSE
00527
00528 #ifdef CPPL_DEBUG
00529 if(M!=N){
00530 std::cerr << "[ERROR] dgematrix::dgeev"
00531 << "(vector<double>&, vector<double>&, "
00532 << "vector<drovector>&, vector<drovector>&) "
00533 << std::endl
00534 << "This matrix is not a square matrix." << std::endl
00535 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00536 exit(1);
00537 }
00538 #endif//CPPL_DEBUG
00539
00540 wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
00541 for(long i=0; i<N; i++){ vlr[i].resize(N); vli[i].resize(N); }
00542 dgematrix VL(N,N);
00543 char JOBVL('V'), JOBVR('N');
00544 long LDA(N), LDVL(N), LDVR(1), LWORK(4*N), INFO(1);
00545 double *VR(NULL), *WORK(new double[LWORK]);
00546 dgeev_(JOBVL, JOBVR, N, Array, LDA, &wr[0], &wi[0],
00547 VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
00548 delete [] WORK; delete [] VR;
00549
00551 for(long j=0; j<N; j++){
00552 if(fabs(wi[j])<1e-10){
00553 for(long i=0; i<N; i++){
00554 vlr[j](i) = VL(i,j); vli[j](i) = 0.0;
00555 }
00556 }
00557 else{
00558 for(long i=0; i<N; i++){
00559 vlr[j](i) = VL(i,j); vli[j](i) =-VL(i,j+1);
00560 vlr[j+1](i) = VL(i,j); vli[j+1](i) = VL(i,j+1);
00561 }
00562 j++;
00563 }
00564 }
00565
00566
00567 if(INFO!=0){
00568 std::cerr << "[WARNING] dgematrix::dgeev"
00569 << "(vector<double>&, vector<double>&, "
00570 << "vector<drovector>&, vector<drovector>&) "
00571 << std::endl
00572 << "Serious trouble happend. INFO = " << INFO << "."
00573 << std::endl;
00574 }
00575 return INFO;
00576 }
00577
00578
00579
00583
00584
00591 inline long dgematrix::dggev(dgematrix& matB,
00592 std::vector<double>& wr, std::vector<double>& wi)
00593 {
00594 #ifdef CPPL_VERBOSE
00595 std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&)"
00596 << std::endl;
00597 #endif//CPPL_VERBOSE
00598
00599 #ifdef CPPL_DEBUG
00600 if(M!=N){
00601 std::cerr << "[ERROR] dgematrix::dggev"
00602 << "(dgematrix&, vector<double>&, vector<double>&) " << std::endl
00603 << "This matrix is not a square matrix." << std::endl
00604 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00605 exit(1);
00606 }
00607 if(matB.M!=N || matB.N!=N){
00608 std::cerr << "[ERROR] dgematrix::dggev"
00609 << "(dgematrix&, vector<double>&, vector<double>&) " << std::endl
00610 << "The matrix B is not a square matrix "
00611 << "having the same size as \"this\" matrix." << std::endl
00612 << "The B matrix is (" << matB.M << "x" << matB.N << ")."
00613 << std::endl;
00614 exit(1);
00615 }
00616 #endif//CPPL_DEBUG
00617
00618 wr.resize(N); wi.resize(N);
00619 char JOBVL('N'), JOBVR('N');
00620 long LDA(N), LDB(N), LDVL(1), LDVR(1), LWORK(8*N), INFO(1);
00621 double *BETA(new double[N]), *VL(NULL), *VR(NULL),
00622 *WORK(new double[LWORK]);
00623 dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
00624 BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO);
00625 delete [] WORK; delete [] VL; delete [] VR;
00626
00628 for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
00629 delete [] BETA;
00630
00631 if(INFO!=0){
00632 std::cerr << "[WARNING] dgematrix::dgeev"
00633 << "(dgematrix&, vector<double>&, vector<double>&)"
00634 << std::endl
00635 << "Serious trouble happend. INFO = " << INFO << "."
00636 << std::endl;
00637 }
00638 return INFO;
00639 }
00640
00641
00649 inline long dgematrix::dggev(dgematrix& matB,
00650 std::vector<double>& wr, std::vector<double>& wi,
00651 std::vector<dcovector>& vrr,
00652 std::vector<dcovector>& vri)
00653 {
00654 #ifdef CPPL_VERBOSE
00655 std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<dcovector>&, std::vector<dcovector>&)"
00656 << std::endl;
00657 #endif//CPPL_VERBOSE
00658
00659 #ifdef CPPL_DEBUG
00660 if(M!=N){
00661 std::cerr << "[ERROR] dgematrix::dggev"
00662 << "(dgematrix&, vector<double>&, vector<double>&, "
00663 << "vector<dcovector>&, vector<dcovector>&)" << std::endl
00664 << "This matrix is not a square matrix." << std::endl
00665 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00666 exit(1);
00667 }
00668 if(matB.M!=N || matB.N!=N){
00669 std::cerr << "[ERROR] dgematrix::dggev"
00670 << "(dgematrix&, vector<double>&, vector<double>&, "
00671 << "vector<dcovector>&, vector<dcovector>&)" << std::endl
00672 << "The matrix B is not a square matrix "
00673 << "having the same size as \"this\" matrix." << std::endl
00674 << "The B matrix is (" << matB.M << "x" << matB.N << ")."
00675 << std::endl;
00676 exit(1);
00677 }
00678 #endif//CPPL_DEBUG
00679
00680 wr.resize(N); wi.resize(N); vrr.resize(N); vri.resize(N);
00681 for(long i=0; i<N; i++){ vrr[i].resize(N); vri[i].resize(N); }
00682 dgematrix VR(N,N);
00683 char JOBVL('N'), JOBVR('V');
00684 long LDA(N), LDB(N), LDVL(1), LDVR(N), LWORK(8*N), INFO(1);
00685 double *BETA(new double[N]), *VL(NULL), *WORK(new double[LWORK]);
00686 dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
00687 BETA, VL, LDVL, VR.Array, LDVR, WORK, LWORK, INFO);
00688 delete [] WORK; delete [] VL;
00689
00691 for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
00692 delete [] BETA;
00693
00695 for(long j=0; j<N; j++){
00696 if(fabs(wi[j])<1e-10){
00697 for(long i=0; i<N; i++){
00698 vrr[j](i) = VR(i,j); vri[j](i) = 0.0;
00699 }
00700 }
00701 else{
00702 for(long i=0; i<N; i++){
00703 vrr[j](i) = VR(i,j); vri[j](i) = VR(i,j+1);
00704 vrr[j+1](i) = VR(i,j); vri[j+1](i) =-VR(i,j+1);
00705 }
00706 j++;
00707 }
00708 }
00709
00710 if(INFO!=0){
00711 std::cerr << "[WARNING] dgematrix::dggev"
00712 << "(dgematrix&, vector<double>&, vector<double>&, "
00713 << "vector<dcovector>&, vector<dcovector>&)" << std::endl
00714 << "Serious trouble happend. INFO = " << INFO << "."
00715 << std::endl;
00716 }
00717 return INFO;
00718 }
00719
00720
00728 inline long dgematrix::dggev(dgematrix& matB,
00729 std::vector<double>& wr, std::vector<double>& wi,
00730 std::vector<drovector>& vlr,
00731 std::vector<drovector>& vli)
00732 {
00733 #ifdef CPPL_VERBOSE
00734 std::cerr << "# [MARK] dgematrix::dggev(dgematrix&, std::vector<double>&, std::vector<double>&, std::vector<drovector>&, std::vector<drovector>&)"
00735 << std::endl;
00736 #endif//CPPL_VERBOSE
00737
00738 #ifdef CPPL_DEBUG
00739 if(M!=N){
00740 std::cerr << "[ERROR] dgematrix::dggev"
00741 << "(dgematrix&, vector<double>&, vector<double>&, "
00742 << "vector<drovector>&, vector<drovector>&)" << std::endl
00743 << "This matrix is not a square matrix." << std::endl
00744 << "This matrix is (" << M << "x" << N << ")." << std::endl;
00745 exit(1);
00746 }
00747 if(matB.M!=N || matB.N!=N){
00748 std::cerr << "[ERROR] dgematrix::dggev"
00749 << "(dgematrix&, vector<double>&, vector<double>&, "
00750 << "vector<drovector>&, vector<drovector>&)" << std::endl
00751 << "The matrix B is not a square matrix "
00752 << "having the same size as \"this\" matrix." << std::endl
00753 << "The B matrix is (" << matB.M << "x" << matB.N << ")."
00754 << std::endl;
00755 exit(1);
00756 }
00757 #endif//CPPL_DEBUG
00758
00759 wr.resize(N); wi.resize(N); vlr.resize(N); vli.resize(N);
00760 for(long i=0; i<N; i++){ vlr[i].resize(N); vli[i].resize(N); }
00761 dgematrix VL(N,N);
00762 char JOBVL('V'), JOBVR('N');
00763 long LDA(N), LDB(N), LDVL(N), LDVR(1), LWORK(8*N), INFO(1);
00764 double *BETA(new double[N]), *VR(NULL), *WORK(new double[LWORK]);
00765 dggev_(JOBVL, JOBVR, N, Array, LDA, matB.Array, LDB, &wr[0], &wi[0],
00766 BETA, VL.Array, LDVL, VR, LDVR, WORK, LWORK, INFO);
00767 delete [] WORK; delete [] VR;
00768
00770 for(long i=0; i<N; i++){ wr[i]/=BETA[i]; wi[i]/=BETA[i]; }
00771 delete [] BETA;
00772
00774 for(long j=0; j<N; j++){
00775 if(fabs(wi[j])<1e-10){
00776 for(long i=0; i<N; i++){
00777 vlr[j](i) = VL(i,j); vli[j](i) = 0.0;
00778 }
00779 }
00780 else{
00781 for(long i=0; i<N; i++){
00782 vlr[j](i) = VL(i,j); vli[j](i) =-VL(i,j+1);
00783 vlr[j+1](i) = VL(i,j); vli[j+1](i) = VL(i,j+1);
00784 }
00785 j++;
00786 }
00787 }
00788
00789 if(INFO!=0){
00790 std::cerr << "[WARNING] dgematrix::dggev"
00791 << "(dgematrix&, vector<double>&, vector<double>&, "
00792 << "vector<drovector>&, vector<drovector>&)" << std::endl
00793 << "Serious trouble happend. INFO = " << INFO << "."
00794 << std::endl;
00795 }
00796 return INFO;
00797 }
00798
00802
00803
00812 inline long dgematrix::dgesvd(dcovector& S, dgematrix& U, dgematrix& VT)
00813 {
00814 #ifdef CPPL_VERBOSE
00815 std::cerr << "# [MARK] dgematrix::dgesvd(dcovector&, dgematrix&, dgematrix&)"
00816 << std::endl;
00817 #endif//CPPL_VERBOSE
00818
00819 char JOBU('A'), JOBVT('A');
00820 long LDA(M), LDU(M), LDVT(N),
00821 LWORK(max(3*min(M,N)+max(M,N),5*min(M,N))), INFO(1);
00822 double *WORK(new double[LWORK]);
00823 S.resize(min(M,N)); U.resize(LDU,M); VT.resize(LDVT,N);
00824
00825 dgesvd_(JOBU, JOBVT, M, N, Array, LDA, S.Array, U.Array,
00826 LDU, VT.Array, LDVT, WORK, LWORK, INFO);
00827 delete [] WORK;
00828
00829 if(INFO!=0){
00830 std::cerr << "[WARNING] dgematrix::dgesvd"
00831 << "(dceovector&, dgematrix&, dcovector&) "
00832 << "Serious trouble happend. INFO = " << INFO << "."
00833 << std::endl;
00834 }
00835 return INFO;
00836 }