IT++ Logo

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

Generated on Sun Sep 14 18:54:50 2008 for IT++ by Doxygen 1.5.6