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

Generated on Tue Feb 2 09:33:28 2010 for IT++ by Doxygen 1.6.2