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
Generated on Thu Apr 19 14:14:57 2007 for IT++ by Doxygen 1.5.1