IT++ Logo
svd.cpp
Go to the documentation of this file.
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/svd.h>
00040 
00041 
00042 namespace itpp
00043 {
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   // The theoretical calculation of lwork above results in the minimum size
00062   // needed for dgesvd_ to run to completion without having memory errors.
00063   // For speed improvement it is best to set lwork=-1 and have dgesvd_
00064   // calculate the best workspace requirement.
00065   int lwork_tmp = -1;
00066   dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00067           V._data(), &ldvt, work._data(), &lwork_tmp, &info);
00068   if (info == 0) {
00069     lwork = static_cast<int>(work(0));
00070     work.set_size(lwork, false);
00071   }
00072 
00073   dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00074           V._data(), &ldvt, work._data(), &lwork, &info);
00075 
00076   return (info == 0);
00077 }
00078 
00079 bool svd(const cmat &A, vec &S)
00080 {
00081   char jobu = 'N', jobvt = 'N';
00082   int m, n, lda, ldu, ldvt, lwork, info;
00083   m = lda = ldu = A.rows();
00084   n = ldvt = A.cols();
00085   lwork = 2 * std::min(m, n) + std::max(m, n);
00086 
00087   cvec U, V;
00088   S.set_size(std::min(m, n), false);
00089   cvec work(lwork);
00090   vec rwork(5*std::min(m, n));
00091 
00092   cmat B(A);
00093 
00094   // The theoretical calculation of lwork above results in the minimum size
00095   // needed for zgesvd_ to run to completion without having memory errors.
00096   // For speed improvement it is best to set lwork=-1 and have zgesvd_
00097   // calculate the best workspace requirement.
00098   int lwork_tmp = -1;
00099   zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00100           V._data(), &ldvt, work._data(), &lwork_tmp, rwork._data(), &info);
00101   if (info == 0) {
00102     lwork = static_cast<int>(real(work(0)));
00103     work.set_size(lwork, false);
00104   }
00105 
00106   zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00107           V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
00108 
00109   return (info == 0);
00110 }
00111 
00112 bool svd(const mat &A, mat &U, vec &S, mat &V)
00113 {
00114   char jobu = 'A', jobvt = 'A';
00115   int m, n, lda, ldu, ldvt, lwork, info;
00116   m = lda = ldu = A.rows();
00117   n = ldvt = A.cols();
00118   lwork = std::max(3 * std::min(m, n) + std::max(m, n), 5 * std::min(m, n));
00119 
00120   U.set_size(m, m, false);
00121   V.set_size(n, n, false);
00122   S.set_size(std::min(m, n), false);
00123   vec work(lwork);
00124 
00125   mat B(A);
00126 
00127   // The theoretical calculation of lwork above results in the minimum size
00128   // needed for dgesvd_ to run to completion without having memory errors.
00129   // For speed improvement it is best to set lwork=-1 and have dgesvd_
00130   // calculate the best workspace requirement.
00131   int lwork_tmp = -1;
00132   dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00133           V._data(), &ldvt, work._data(), &lwork_tmp, &info);
00134   if (info == 0) {
00135     lwork = static_cast<int>(work(0));
00136     work.set_size(lwork, false);
00137   }
00138 
00139   dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00140           V._data(), &ldvt, work._data(), &lwork, &info);
00141 
00142   V = V.T(); // This is probably slow!!!
00143 
00144   return (info == 0);
00145 }
00146 
00147 bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
00148 {
00149   char jobu = 'A', jobvt = 'A';
00150   int m, n, lda, ldu, ldvt, lwork, info;
00151   m = lda = ldu = A.rows();
00152   n = ldvt = A.cols();
00153   lwork = 2 * std::min(m, n) + std::max(m, n);
00154 
00155   U.set_size(m, m, false);
00156   V.set_size(n, n, false);
00157   S.set_size(std::min(m, n), false);
00158   cvec work(lwork);
00159   vec rwork(5 * std::min(m, n));
00160 
00161   cmat B(A);
00162 
00163   // The theoretical calculation of lwork above results in the minimum size
00164   // needed for zgesvd_ to run to completion without having memory errors.
00165   // For speed improvement it is best to set lwork=-1 and have zgesvd_
00166   // calculate the best workspace requirement.
00167   int lwork_tmp = -1;
00168   zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00169           V._data(), &ldvt, work._data(), &lwork_tmp, rwork._data(), &info);
00170   if (info == 0) {
00171     lwork = static_cast<int>(real(work(0)));
00172     work.set_size(lwork, false);
00173   }
00174 
00175   zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu,
00176           V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info);
00177 
00178   V = V.H(); // This is slow!!!
00179 
00180   return (info == 0);
00181 }
00182 
00183 #else
00184 
00185 bool svd(const mat &A, vec &S)
00186 {
00187   it_error("LAPACK library is needed to use svd() function");
00188   return false;
00189 }
00190 
00191 bool svd(const cmat &A, vec &S)
00192 {
00193   it_error("LAPACK library is needed to use svd() function");
00194   return false;
00195 }
00196 
00197 bool svd(const mat &A, mat &U, vec &S, mat &V)
00198 {
00199   it_error("LAPACK library is needed to use svd() function");
00200   return false;
00201 }
00202 
00203 bool svd(const cmat &A, cmat &U, vec &S, cmat &V)
00204 {
00205   it_error("LAPACK library is needed to use svd() function");
00206   return false;
00207 }
00208 
00209 #endif // HAVE_LAPACK
00210 
00211 vec svd(const mat &A)
00212 {
00213   vec S;
00214   svd(A, S);
00215   return S;
00216 }
00217 
00218 vec svd(const cmat &A)
00219 {
00220   vec S;
00221   svd(A, S);
00222   return S;
00223 }
00224 
00225 } //namespace itpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4