IT++ Logo Newcom Logo

eigen.cpp

Go to the documentation of this file.
00001 
00033 #ifndef _MSC_VER
00034 #  include <itpp/config.h>
00035 #else
00036 #  include <itpp/config_msvc.h>
00037 #endif
00038 
00039 #if defined(HAVE_LAPACK)
00040 #  include <itpp/base/lapack.h>
00041 #endif
00042 
00043 #include <itpp/base/eigen.h>
00044 #include <itpp/base/converters.h>
00045 
00046 
00047 namespace itpp { 
00048 
00049 #if defined(HAVE_LAPACK)
00050 
00051   bool eig_sym(const mat &A, vec &d, mat &V)
00052   {
00053     it_assert1(A.rows() == A.cols(), "eig_sym: Matrix is not symmetric");
00054 
00055     // Test for symmetric?
00056 
00057     char jobz='V', uplo='U';
00058     int n, lda, lwork, info;
00059     n = lda = A.rows();
00060     lwork = std::max(1,3*n-1); // This may be choosen better!
00061 
00062     d.set_size(n, false);
00063     vec work(lwork);
00064 
00065     V = A; // The routine overwrites input matrix with eigenvectors
00066 
00067     dsyev_(&jobz, &uplo, &n, V._data(), &lda, d._data(), work._data(), &lwork, &info);
00068 
00069     return (info==0);
00070   }
00071 
00072   bool eig_sym(const mat &A, vec &d)
00073   {
00074     it_assert1(A.rows() == A.cols(), "eig_sym: Matrix is not symmetric");
00075 
00076     // Test for symmetric?
00077 
00078     char jobz='N', uplo='U';
00079     int n, lda, lwork, info;
00080     n = lda = A.rows();
00081     lwork = std::max(1,3*n-1); // This may be choosen better!
00082 
00083     d.set_size(n, false);
00084     vec work(lwork);
00085 
00086     mat B(A); // The routine overwrites input matrix
00087 
00088     dsyev_(&jobz, &uplo, &n, B._data(), &lda, d._data(), work._data(), &lwork, &info);
00089 
00090     return (info==0);
00091   }
00092 
00093   bool eig_sym(const cmat &A, vec &d, cmat &V)
00094   {
00095     it_assert1(A.rows() == A.cols(), "eig_sym: Matrix is not hermitian");
00096 
00097     // Test for symmetric?
00098 
00099     char jobz='V', uplo='U';
00100     int n, lda, lwork, info;
00101     n = lda = A.rows();
00102     lwork = std::max(1,2*n-1); // This may be choosen better!
00103 
00104     d.set_size(n, false);
00105     cvec work(lwork);
00106     vec rwork(std::max(1,3*n-2)); // This may be choosen better!
00107 
00108     V = A; // The routine overwrites input matrix with eigenvectors
00109 
00110     zheev_(&jobz, &uplo, &n, V._data(), &lda, d._data(), work._data(), &lwork, rwork._data(), &info);
00111 
00112     return (info==0);
00113   }
00114 
00115   bool eig_sym(const cmat &A, vec &d)
00116   {
00117     it_assert1(A.rows() == A.cols(), "eig_sym: Matrix is not hermitian");
00118 
00119     // Test for symmetric?
00120 
00121     char jobz='N', uplo='U';
00122     int n, lda, lwork, info;
00123     n = lda = A.rows();
00124     lwork = std::max(1,2*n-1); // This may be choosen better!
00125 
00126     d.set_size(n, false);
00127     cvec work(lwork);
00128     vec rwork(std::max(1,3*n-2)); // This may be choosen better!
00129 
00130     cmat B(A); // The routine overwrites input matrix
00131 
00132     zheev_(&jobz, &uplo, &n, B._data(), &lda, d._data(), work._data(), &lwork, rwork._data(), &info);
00133 
00134     return (info==0);
00135   }
00136 
00137 
00138   // Non-symmetric matrix
00139   bool eig(const mat &A, cvec &d, cmat &V)
00140   {
00141     it_assert1(A.rows() == A.cols(), "eig: Matrix is not square");
00142 
00143     char jobvl='N', jobvr='V';
00144     int n, lda, ldvl, ldvr, lwork, info;
00145     n = lda = A.rows();
00146     ldvl = 1; ldvr = n;
00147     lwork = std::max(1,4*n); // This may be choosen better!
00148 
00149     vec work(lwork);
00150     vec rwork(std::max(1,2*n)); // This may be choosen better
00151     vec wr(n), wi(n);
00152     mat vl, vr(n,n);
00153 
00154     mat B(A); // The routine overwrites input matrix
00155 
00156     dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info);
00157 
00158     d = to_cvec(wr, wi);
00159 
00160     // Fix V
00161     V.set_size(n, n, false);
00162     for (int j=0; j<n; j++) {
00163       // if d(j) and d(j+1) are complex conjugate pairs, treat special
00164                 if( (j<n-1) && d(j) == std::conj(d(j+1))) {
00165         V.set_col(j, to_cvec(vr.get_col(j), vr.get_col(j+1)) );
00166         V.set_col(j+1, to_cvec(vr.get_col(j), -vr.get_col(j+1)) );
00167         j++;
00168       } else {
00169         V.set_col(j, to_cvec(vr.get_col(j)) );
00170       }
00171     }
00172 
00173     return (info==0);
00174   }
00175 
00176   // Non-symmetric matrix
00177   bool eig(const mat &A, cvec &d)
00178   {
00179     it_assert1(A.rows() == A.cols(), "eig: Matrix is not square");
00180 
00181     char jobvl='N', jobvr='N';
00182     int n, lda, ldvl, ldvr, lwork, info;
00183     n = lda = A.rows();
00184     ldvl = 1; ldvr = 1;
00185     lwork = std::max(1,4*n); // This may be choosen better!
00186 
00187     vec work(lwork);
00188     vec rwork(std::max(1,2*n)); // This may be choosen better
00189     vec wr(n), wi(n);
00190     mat vl, vr;
00191 
00192     mat B(A); // The routine overwrites input matrix
00193 
00194     dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info);
00195 
00196     d = to_cvec(wr, wi);
00197 
00198     return (info==0);
00199   }
00200 
00201   bool eig(const cmat &A, cvec &d, cmat &V)
00202   {
00203     it_assert1(A.rows() == A.cols(), "eig: Matrix is not square");
00204 
00205     char jobvl='N', jobvr='V';
00206     int n, lda, ldvl, ldvr, lwork, info;
00207     n = lda = A.rows();
00208     ldvl = 1; ldvr = n;
00209     lwork = std::max(1,2*n); // This may be choosen better!
00210 
00211     d.set_size(n, false);
00212     V.set_size(n, n, false);
00213     cvec work(lwork);
00214     vec rwork(std::max(1,2*n)); // This may be choosen better!
00215     cmat vl;
00216 
00217     cmat B(A); // The routine overwrites input matrix
00218 
00219     zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, V._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
00220 
00221 
00222     return (info==0);
00223   }
00224 
00225   bool eig(const cmat &A, cvec &d)
00226   {
00227     it_assert1(A.rows() == A.cols(), "eig: Matrix is not square");
00228 
00229     char jobvl='N', jobvr='N';
00230     int n, lda, ldvl, ldvr, lwork, info;
00231     n = lda = A.rows();
00232     ldvl = 1; ldvr = 1;
00233     lwork = std::max(1,2*n); // This may be choosen better!
00234 
00235     d.set_size(n, false);
00236     cvec work(lwork);
00237     vec rwork(std::max(1,2*n)); // This may be choosen better!
00238     cmat vl, vr;
00239 
00240     cmat B(A); // The routine overwrites input matrix
00241 
00242     zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
00243 
00244 
00245     return (info==0);
00246   }
00247 
00248 #else
00249 
00250   bool eig_sym(const mat &A, vec &d, mat &V)
00251   {
00252     it_error("LAPACK library is needed to use eig_sym() function");
00253     return false;   
00254   }
00255 
00256   bool eig_sym(const mat &A, vec &d)
00257   {
00258     it_error("LAPACK library is needed to use eig_sym() function");
00259     return false;   
00260   }
00261 
00262   bool eig_sym(const cmat &A, vec &d, cmat &V)
00263   { 
00264     it_error("LAPACK library is needed to use eig_sym() function");
00265     return false;   
00266   }
00267 
00268   bool eig_sym(const cmat &A, vec &d)
00269   {
00270     it_error("LAPACK library is needed to use eig_sym() function");
00271     return false;   
00272   }
00273 
00274 
00275   bool eig(const mat &A, cvec &d, cmat &V)
00276   {
00277     it_error("LAPACK library is needed to use eig() function");
00278     return false;   
00279   }
00280 
00281   bool eig(const mat &A, cvec &d)
00282   {
00283     it_error("LAPACK library is needed to use eig() function");
00284     return false;   
00285   }
00286 
00287   bool eig(const cmat &A, cvec &d, cmat &V)
00288   {
00289     it_error("LAPACK library is needed to use eig() function");
00290     return false;   
00291   }
00292 
00293   bool eig(const cmat &A, cvec &d)
00294   {
00295     it_error("LAPACK library is needed to use eig() function");
00296     return false;   
00297   }
00298 
00299 #endif // HAVE_LAPACK
00300 
00301   vec eig_sym(const mat &A)
00302   {
00303     vec d;
00304     eig_sym(A, d);
00305     return d;
00306   }
00307 
00308   vec eig_sym(const cmat &A)
00309   {
00310     vec d;
00311     eig_sym(A, d);
00312     return d;
00313   }
00314 
00315 
00316   cvec eig(const mat &A)
00317   {
00318     cvec d;
00319     eig(A, d);
00320     return d;
00321   }
00322 
00323   cvec eig(const cmat &A)
00324   {
00325     cvec d;
00326     eig(A, d);
00327     return d;
00328   }
00329 
00330 } //namespace itpp
SourceForge Logo

Generated on Sat Aug 25 23:40:51 2007 for IT++ by Doxygen 1.5.2