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
Generated on Sat Apr 19 10:57:50 2008 for IT++ by Doxygen 1.5.5