IT++ Logo Newcom Logo

schur.cpp

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

Generated on Fri Jun 8 02:08:53 2007 for IT++ by Doxygen 1.5.2