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/schur.h> 00044 00045 00046 namespace itpp { 00047 00048 #if defined(HAVE_LAPACK) 00049 00050 bool schur(const mat &A, mat &U, mat &T) { 00051 it_assert1(A.rows() == A.cols(), "schur(): Matrix is not square"); 00052 00053 char jobvs = 'V'; 00054 char sort = 'N'; 00055 int info; 00056 int n = A.rows(); 00057 int lda = n; 00058 int ldvs = n; 00059 int lwork = 3 * n; // This may be choosen better! 00060 int sdim = 0; 00061 vec wr(n); 00062 vec wi(n); 00063 vec work(lwork); 00064 00065 T.set_size(lda, n, false); 00066 U.set_size(ldvs, n, false); 00067 00068 T = A; // The routine overwrites input matrix with eigenvectors 00069 00070 dgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, wr._data(), wi._data(), 00071 U._data(), &ldvs, work._data(), &lwork, 0, &info); 00072 00073 return (info == 0); 00074 } 00075 00076 00077 bool schur(const cmat &A, cmat &U, cmat &T) { 00078 it_assert1(A.rows() == A.cols(), "schur(): Matrix is not square"); 00079 00080 char jobvs = 'V'; 00081 char sort = 'N'; 00082 int info; 00083 int n = A.rows(); 00084 int lda = n; 00085 int ldvs = n; 00086 int lwork = 2 * n; // This may be choosen better! 00087 int sdim = 0; 00088 vec rwork(n); 00089 cvec w(n); 00090 cvec work(lwork); 00091 00092 T.set_size(lda, n, false); 00093 U.set_size(ldvs, n, false); 00094 00095 T = A; // The routine overwrites input matrix with eigenvectors 00096 00097 zgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, w._data(), U._data(), 00098 &ldvs, work._data(), &lwork, rwork._data(), 0, &info); 00099 00100 return (info == 0); 00101 } 00102 00103 #else 00104 00105 bool schur(const mat &A, mat &U, mat &T) { 00106 it_error("LAPACK library is needed to use schur() function"); 00107 return false; 00108 } 00109 00110 00111 bool schur(const cmat &A, cmat &U, cmat &T) { 00112 it_error("LAPACK library is needed to use schur() function"); 00113 return false; 00114 } 00115 00116 #endif // HAVE_LAPACK 00117 00118 mat schur(const mat &A) { 00119 mat U, T; 00120 schur(A, U, T); 00121 return T; 00122 } 00123 00124 00125 cmat schur(const cmat &A) { 00126 cmat U, T; 00127 schur(A, U, T); 00128 return T; 00129 } 00130 00131 } // namespace itpp
Generated on Fri Jun 8 01:07:10 2007 for IT++ by Doxygen 1.5.2