IT++ Logo

schur.cpp

Go to the documentation of this file.
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
SourceForge Logo

Generated on Sat Apr 19 10:41:54 2008 for IT++ by Doxygen 1.5.5