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