IT++ Logo
eigen.cpp
Go to the documentation of this file.
00001 
00029 #ifndef _MSC_VER
00030 #  include <itpp/config.h>
00031 #else
00032 #  include <itpp/config_msvc.h>
00033 #endif
00034 
00035 #if defined(HAVE_LAPACK)
00036 #  include <itpp/base/algebra/lapack.h>
00037 #endif
00038 
00039 #include <itpp/base/algebra/eigen.h>
00040 #include <itpp/base/converters.h>
00041 
00042 
00043 namespace itpp
00044 {
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;
00144   ldvr = n;
00145   lwork = std::max(1, 4 * n); // This may be choosen better!
00146 
00147   vec work(lwork);
00148   vec rwork(std::max(1, 2*n)); // This may be choosen better
00149   vec wr(n), wi(n);
00150   mat vl, vr(n, n);
00151 
00152   mat B(A); // The routine overwrites input matrix
00153 
00154   dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info);
00155 
00156   d = to_cvec(wr, wi);
00157 
00158   // Fix V
00159   V.set_size(n, n, false);
00160   for (int j = 0; j < n; j++) {
00161     // if d(j) and d(j+1) are complex conjugate pairs, treat special
00162     if ((j < n - 1) && d(j) == std::conj(d(j + 1))) {
00163       V.set_col(j, to_cvec(vr.get_col(j), vr.get_col(j + 1)));
00164       V.set_col(j + 1, to_cvec(vr.get_col(j), -vr.get_col(j + 1)));
00165       j++;
00166     }
00167     else {
00168       V.set_col(j, to_cvec(vr.get_col(j)));
00169     }
00170   }
00171 
00172   return (info == 0);
00173 }
00174 
00175 // Non-symmetric matrix
00176 bool eig(const mat &A, cvec &d)
00177 {
00178   it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
00179 
00180   char jobvl = 'N', jobvr = 'N';
00181   int n, lda, ldvl, ldvr, lwork, info;
00182   n = lda = A.rows();
00183   ldvl = 1;
00184   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_assert_debug(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;
00209   ldvr = n;
00210   lwork = std::max(1, 2 * n); // This may be choosen better!
00211 
00212   d.set_size(n, false);
00213   V.set_size(n, n, false);
00214   cvec work(lwork);
00215   vec rwork(std::max(1, 2*n)); // This may be choosen better!
00216   cmat vl;
00217 
00218   cmat B(A); // The routine overwrites input matrix
00219 
00220   zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, V._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
00221 
00222 
00223   return (info == 0);
00224 }
00225 
00226 bool eig(const cmat &A, cvec &d)
00227 {
00228   it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
00229 
00230   char jobvl = 'N', jobvr = 'N';
00231   int n, lda, ldvl, ldvr, lwork, info;
00232   n = lda = A.rows();
00233   ldvl = 1;
00234   ldvr = 1;
00235   lwork = std::max(1, 2 * n); // This may be choosen better!
00236 
00237   d.set_size(n, false);
00238   cvec work(lwork);
00239   vec rwork(std::max(1, 2*n)); // This may be choosen better!
00240   cmat vl, vr;
00241 
00242   cmat B(A); // The routine overwrites input matrix
00243 
00244   zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
00245 
00246 
00247   return (info == 0);
00248 }
00249 
00250 #else
00251 
00252 bool eig_sym(const mat &A, vec &d, mat &V)
00253 {
00254   it_error("LAPACK library is needed to use eig_sym() function");
00255   return false;
00256 }
00257 
00258 bool eig_sym(const mat &A, vec &d)
00259 {
00260   it_error("LAPACK library is needed to use eig_sym() function");
00261   return false;
00262 }
00263 
00264 bool eig_sym(const cmat &A, vec &d, cmat &V)
00265 {
00266   it_error("LAPACK library is needed to use eig_sym() function");
00267   return false;
00268 }
00269 
00270 bool eig_sym(const cmat &A, vec &d)
00271 {
00272   it_error("LAPACK library is needed to use eig_sym() function");
00273   return false;
00274 }
00275 
00276 
00277 bool eig(const mat &A, cvec &d, cmat &V)
00278 {
00279   it_error("LAPACK library is needed to use eig() function");
00280   return false;
00281 }
00282 
00283 bool eig(const mat &A, cvec &d)
00284 {
00285   it_error("LAPACK library is needed to use eig() function");
00286   return false;
00287 }
00288 
00289 bool eig(const cmat &A, cvec &d, cmat &V)
00290 {
00291   it_error("LAPACK library is needed to use eig() function");
00292   return false;
00293 }
00294 
00295 bool eig(const cmat &A, cvec &d)
00296 {
00297   it_error("LAPACK library is needed to use eig() function");
00298   return false;
00299 }
00300 
00301 #endif // HAVE_LAPACK
00302 
00303 vec eig_sym(const mat &A)
00304 {
00305   vec d;
00306   eig_sym(A, d);
00307   return d;
00308 }
00309 
00310 vec eig_sym(const cmat &A)
00311 {
00312   vec d;
00313   eig_sym(A, d);
00314   return d;
00315 }
00316 
00317 
00318 cvec eig(const mat &A)
00319 {
00320   cvec d;
00321   eig(A, d);
00322   return d;
00323 }
00324 
00325 cvec eig(const cmat &A)
00326 {
00327   cvec d;
00328   eig(A, d);
00329   return d;
00330 }
00331 
00332 } //namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4