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
Generated on Thu Apr 23 20:04:03 2009 for IT++ by Doxygen 1.5.8