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     // 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
SourceForge Logo

Generated on Sat Apr 19 10:57:50 2008 for IT++ by Doxygen 1.5.5