IT++ Logo Newcom Logo

svd.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/svd.h>
00044 
00045 
00046 namespace itpp { 
00047 
00048 #if defined(HAVE_LAPACK)
00049 
00050   bool svd(const mat &A, vec &S)
00051   {
00052     char jobu='N', jobvt='N';
00053     int m, n, lda, ldu, ldvt, lwork, info;
00054     m = lda = ldu = A.rows();
00055     n = ldvt = A.cols();
00056     lwork = std::max(3*std::min(m,n)+std::max(m,n), 5*std::min(m,n));
00057 
00058     mat U, V;
00059     S.set_size(std::min(m,n), false);
00060     vec work(lwork);
00061 
00062     mat B(A);
00063 
00064     dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, &info);
00065 
00066     return (info==0);
00067   }
00068 
00069   bool svd(const cmat &A, vec &S)
00070   {
00071     char jobu='N', jobvt='N';
00072     int m, n, lda, ldu, ldvt, lwork, info;
00073     m = lda = ldu = A.rows();
00074     n = ldvt = A.cols();
00075     lwork = 2*std::min(m,n)+std::max(m,n);
00076 
00077     cvec U, V;
00078     //U.set_size(m,m, false);
00079     //V.set_size(n,n, false);
00080     S.set_size(std::min(m,n), false);
00081     cvec work(lwork);
00082     vec rwork(std::max(1, 5*std::min(m, n)));
00083 
00084     cmat B(A);
00085 
00086     zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
00087 
00088     return (info==0);
00089   }
00090 
00091   bool svd(const mat &A, mat &U, vec &S, mat &V)
00092   {
00093     char jobu='A', jobvt='A';
00094     int m, n, lda, ldu, ldvt, lwork, info;
00095     m = lda = ldu = A.rows();
00096     n = ldvt = A.cols();
00097     lwork = std::max(3*std::min(m,n)+std::max(m,n), 5*std::min(m,n));
00098 
00099     U.set_size(m,m, false);
00100     V.set_size(n,n, false);
00101     S.set_size(std::min(m,n), false);
00102     vec work(lwork);
00103 
00104     mat B(A);
00105 
00106     dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, &info);
00107 
00108     V = V.T(); // This is probably slow!!!
00109 
00110     return (info==0);
00111   }
00112 
00113   bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
00114   {
00115     char jobu='A', jobvt='A';
00116     int m, n, lda, ldu, ldvt, lwork, info;
00117     m = lda = ldu = A.rows();
00118     n = ldvt = A.cols();
00119     lwork = 2*std::min(m,n)+std::max(m,n);
00120 
00121     U.set_size(m,m, false);
00122     V.set_size(n,n, false);
00123     S.set_size(std::min(m,n), false);
00124     cvec work(lwork);
00125     vec rwork(std::max(1, 5*std::min(m, n)));
00126 
00127     cmat B(A);
00128 
00129     zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
00130 
00131     V = V.H(); // This is slow!!!
00132 
00133     return (info==0);
00134   }
00135 
00136 #else
00137 
00138   bool svd(const mat &A, vec &S)
00139   {
00140     it_error("LAPACK library is needed to use svd() function");
00141     return false;
00142   }
00143 
00144   bool svd(const cmat &A, vec &S)
00145   {
00146     it_error("LAPACK library is needed to use svd() function");
00147     return false;
00148   }
00149 
00150   bool svd(const mat &A, mat &U, vec &S, mat &V)
00151   {   
00152     it_error("LAPACK library is needed to use svd() function");
00153     return false;
00154   }
00155 
00156   bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
00157   {
00158     it_error("LAPACK library is needed to use svd() function");
00159     return false;   
00160   }
00161 
00162 #endif // HAVE_LAPACK
00163 
00164   vec svd(const mat &A)
00165   {
00166     vec S;
00167     svd(A, S);
00168     return S;
00169   }
00170 
00171   vec svd(const cmat &A)
00172   {
00173     vec S;
00174     svd(A, S);
00175     return S;
00176   }
00177 
00178 } //namespace itpp
SourceForge Logo

Generated on Fri Jan 11 08:51:37 2008 for IT++ by Doxygen 1.3.9.1