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
Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4