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 #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; ldvr = n; 00144 lwork = std::max(1,4*n); // This may be choosen better! 00145 00146 vec work(lwork); 00147 vec rwork(std::max(1,2*n)); // This may be choosen better 00148 vec wr(n), wi(n); 00149 mat vl, vr(n,n); 00150 00151 mat B(A); // The routine overwrites input matrix 00152 00153 dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info); 00154 00155 d = to_cvec(wr, wi); 00156 00157 // Fix V 00158 V.set_size(n, n, false); 00159 for (int j=0; j<n; j++) { 00160 // if d(j) and d(j+1) are complex conjugate pairs, treat special 00161 if( (j<n-1) && d(j) == std::conj(d(j+1))) { 00162 V.set_col(j, to_cvec(vr.get_col(j), vr.get_col(j+1)) ); 00163 V.set_col(j+1, to_cvec(vr.get_col(j), -vr.get_col(j+1)) ); 00164 j++; 00165 } else { 00166 V.set_col(j, to_cvec(vr.get_col(j)) ); 00167 } 00168 } 00169 00170 return (info==0); 00171 } 00172 00173 // Non-symmetric matrix 00174 bool eig(const mat &A, cvec &d) 00175 { 00176 it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square"); 00177 00178 char jobvl='N', jobvr='N'; 00179 int n, lda, ldvl, ldvr, lwork, info; 00180 n = lda = A.rows(); 00181 ldvl = 1; ldvr = 1; 00182 lwork = std::max(1,4*n); // This may be choosen better! 00183 00184 vec work(lwork); 00185 vec rwork(std::max(1,2*n)); // This may be choosen better 00186 vec wr(n), wi(n); 00187 mat vl, vr; 00188 00189 mat B(A); // The routine overwrites input matrix 00190 00191 dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info); 00192 00193 d = to_cvec(wr, wi); 00194 00195 return (info==0); 00196 } 00197 00198 bool eig(const cmat &A, cvec &d, cmat &V) 00199 { 00200 it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square"); 00201 00202 char jobvl='N', jobvr='V'; 00203 int n, lda, ldvl, ldvr, lwork, info; 00204 n = lda = A.rows(); 00205 ldvl = 1; ldvr = n; 00206 lwork = std::max(1,2*n); // This may be choosen better! 00207 00208 d.set_size(n, false); 00209 V.set_size(n, n, false); 00210 cvec work(lwork); 00211 vec rwork(std::max(1,2*n)); // This may be choosen better! 00212 cmat vl; 00213 00214 cmat B(A); // The routine overwrites input matrix 00215 00216 zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, V._data(), &ldvr, work._data(), &lwork, rwork._data(), &info); 00217 00218 00219 return (info==0); 00220 } 00221 00222 bool eig(const cmat &A, cvec &d) 00223 { 00224 it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square"); 00225 00226 char jobvl='N', jobvr='N'; 00227 int n, lda, ldvl, ldvr, lwork, info; 00228 n = lda = A.rows(); 00229 ldvl = 1; ldvr = 1; 00230 lwork = std::max(1,2*n); // This may be choosen better! 00231 00232 d.set_size(n, false); 00233 cvec work(lwork); 00234 vec rwork(std::max(1,2*n)); // This may be choosen better! 00235 cmat vl, vr; 00236 00237 cmat B(A); // The routine overwrites input matrix 00238 00239 zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, rwork._data(), &info); 00240 00241 00242 return (info==0); 00243 } 00244 00245 #else 00246 00247 bool eig_sym(const mat &A, vec &d, mat &V) 00248 { 00249 it_error("LAPACK library is needed to use eig_sym() function"); 00250 return false; 00251 } 00252 00253 bool eig_sym(const mat &A, vec &d) 00254 { 00255 it_error("LAPACK library is needed to use eig_sym() function"); 00256 return false; 00257 } 00258 00259 bool eig_sym(const cmat &A, vec &d, cmat &V) 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) 00266 { 00267 it_error("LAPACK library is needed to use eig_sym() function"); 00268 return false; 00269 } 00270 00271 00272 bool eig(const mat &A, cvec &d, cmat &V) 00273 { 00274 it_error("LAPACK library is needed to use eig() function"); 00275 return false; 00276 } 00277 00278 bool eig(const mat &A, cvec &d) 00279 { 00280 it_error("LAPACK library is needed to use eig() function"); 00281 return false; 00282 } 00283 00284 bool eig(const cmat &A, cvec &d, cmat &V) 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) 00291 { 00292 it_error("LAPACK library is needed to use eig() function"); 00293 return false; 00294 } 00295 00296 #endif // HAVE_LAPACK 00297 00298 vec eig_sym(const mat &A) 00299 { 00300 vec d; 00301 eig_sym(A, d); 00302 return d; 00303 } 00304 00305 vec eig_sym(const cmat &A) 00306 { 00307 vec d; 00308 eig_sym(A, d); 00309 return d; 00310 } 00311 00312 00313 cvec eig(const mat &A) 00314 { 00315 cvec d; 00316 eig(A, d); 00317 return d; 00318 } 00319 00320 cvec eig(const cmat &A) 00321 { 00322 cvec d; 00323 eig(A, d); 00324 return d; 00325 } 00326 00327 } //namespace itpp
Generated on Sun Sep 14 18:57:00 2008 for IT++ by Doxygen 1.5.6