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