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 00046 #if defined(HAVE_LAPACK) 00047 00048 bool schur(const mat &A, mat &U, mat &T) 00049 { 00050 it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square"); 00051 00052 char jobvs = 'V'; 00053 char sort = 'N'; 00054 int info; 00055 int n = A.rows(); 00056 int lda = n; 00057 int ldvs = n; 00058 int lwork = 3 * n; // This may be choosen better! 00059 int sdim = 0; 00060 vec wr(n); 00061 vec wi(n); 00062 vec work(lwork); 00063 00064 T.set_size(lda, n, false); 00065 U.set_size(ldvs, n, false); 00066 00067 T = A; // The routine overwrites input matrix with eigenvectors 00068 00069 dgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, wr._data(), wi._data(), 00070 U._data(), &ldvs, work._data(), &lwork, 0, &info); 00071 00072 return (info == 0); 00073 } 00074 00075 00076 bool schur(const cmat &A, cmat &U, cmat &T) 00077 { 00078 it_assert_debug(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 { 00107 it_error("LAPACK library is needed to use schur() function"); 00108 return false; 00109 } 00110 00111 00112 bool schur(const cmat &A, cmat &U, cmat &T) 00113 { 00114 it_error("LAPACK library is needed to use schur() function"); 00115 return false; 00116 } 00117 00118 #endif // HAVE_LAPACK 00119 00120 mat schur(const mat &A) 00121 { 00122 mat U, T; 00123 schur(A, U, T); 00124 return T; 00125 } 00126 00127 00128 cmat schur(const cmat &A) 00129 { 00130 cmat U, T; 00131 schur(A, U, T); 00132 return T; 00133 } 00134 00135 } // namespace itpp
Generated on Fri Aug 14 15:28:04 2009 for IT++ by Doxygen 1.5.9