IT++ Logo

eigen.cpp

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

Generated on Sun Dec 9 17:30:58 2007 for IT++ by Doxygen 1.5.4