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

Generated on Thu Apr 23 20:06:46 2009 for IT++ by Doxygen 1.5.8