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
Generated on Thu Apr 23 20:06:46 2009 for IT++ by Doxygen 1.5.8