00001 00029 #ifndef _MSC_VER 00030 # include <itpp/config.h> 00031 #else 00032 # include <itpp/config_msvc.h> 00033 #endif 00034 00035 #if defined(HAVE_LAPACK) 00036 # include <itpp/base/algebra/lapack.h> 00037 #endif 00038 00039 #include <itpp/base/algebra/schur.h> 00040 00041 00042 namespace itpp 00043 { 00044 00045 #if defined(HAVE_LAPACK) 00046 00047 bool schur(const mat &A, mat &U, mat &T) 00048 { 00049 it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square"); 00050 00051 char jobvs = 'V'; 00052 char sort = 'N'; 00053 int info; 00054 int n = A.rows(); 00055 int lda = n; 00056 int ldvs = n; 00057 int lwork = 3 * n; // This may be choosen better! 00058 int sdim = 0; 00059 vec wr(n); 00060 vec wi(n); 00061 vec work(lwork); 00062 00063 T.set_size(lda, n, false); 00064 U.set_size(ldvs, n, false); 00065 00066 T = A; // The routine overwrites input matrix with eigenvectors 00067 00068 dgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, wr._data(), wi._data(), 00069 U._data(), &ldvs, work._data(), &lwork, 0, &info); 00070 00071 return (info == 0); 00072 } 00073 00074 00075 bool schur(const cmat &A, cmat &U, cmat &T) 00076 { 00077 it_assert_debug(A.rows() == A.cols(), "schur(): Matrix is not square"); 00078 00079 char jobvs = 'V'; 00080 char sort = 'N'; 00081 int info; 00082 int n = A.rows(); 00083 int lda = n; 00084 int ldvs = n; 00085 int lwork = 2 * n; // This may be choosen better! 00086 int sdim = 0; 00087 vec rwork(n); 00088 cvec w(n); 00089 cvec work(lwork); 00090 00091 T.set_size(lda, n, false); 00092 U.set_size(ldvs, n, false); 00093 00094 T = A; // The routine overwrites input matrix with eigenvectors 00095 00096 zgees_(&jobvs, &sort, 0, &n, T._data(), &lda, &sdim, w._data(), U._data(), 00097 &ldvs, work._data(), &lwork, rwork._data(), 0, &info); 00098 00099 return (info == 0); 00100 } 00101 00102 #else 00103 00104 bool schur(const mat &A, mat &U, mat &T) 00105 { 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 { 00113 it_error("LAPACK library is needed to use schur() function"); 00114 return false; 00115 } 00116 00117 #endif // HAVE_LAPACK 00118 00119 mat schur(const mat &A) 00120 { 00121 mat U, T; 00122 schur(A, U, T); 00123 return T; 00124 } 00125 00126 00127 cmat schur(const cmat &A) 00128 { 00129 cmat U, T; 00130 schur(A, U, T); 00131 return T; 00132 } 00133 00134 } // namespace itpp
Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4