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 dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, &info); 00062 00063 return (info==0); 00064 } 00065 00066 bool svd(const cmat &A, vec &S) 00067 { 00068 char jobu='N', jobvt='N'; 00069 int m, n, lda, ldu, ldvt, lwork, info; 00070 m = lda = ldu = A.rows(); 00071 n = ldvt = A.cols(); 00072 lwork = 2*std::min(m,n)+std::max(m,n); 00073 00074 cvec U, V; 00075 //U.set_size(m,m, false); 00076 //V.set_size(n,n, false); 00077 S.set_size(std::min(m,n), false); 00078 cvec work(lwork); 00079 vec rwork(std::max(1, 5*std::min(m, n))); 00080 00081 cmat B(A); 00082 00083 zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info); 00084 00085 return (info==0); 00086 } 00087 00088 bool svd(const mat &A, mat &U, vec &S, mat &V) 00089 { 00090 char jobu='A', jobvt='A'; 00091 int m, n, lda, ldu, ldvt, lwork, info; 00092 m = lda = ldu = A.rows(); 00093 n = ldvt = A.cols(); 00094 lwork = std::max(3*std::min(m,n)+std::max(m,n), 5*std::min(m,n)); 00095 00096 U.set_size(m,m, false); 00097 V.set_size(n,n, false); 00098 S.set_size(std::min(m,n), false); 00099 vec work(lwork); 00100 00101 mat B(A); 00102 00103 dgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, &info); 00104 00105 V = V.T(); // This is probably slow!!! 00106 00107 return (info==0); 00108 } 00109 00110 bool svd(const cmat &A, cmat &U, vec &S, cmat &V) 00111 { 00112 char jobu='A', jobvt='A'; 00113 int m, n, lda, ldu, ldvt, lwork, info; 00114 m = lda = ldu = A.rows(); 00115 n = ldvt = A.cols(); 00116 lwork = 2*std::min(m,n)+std::max(m,n); 00117 00118 U.set_size(m,m, false); 00119 V.set_size(n,n, false); 00120 S.set_size(std::min(m,n), false); 00121 cvec work(lwork); 00122 vec rwork(std::max(1, 5*std::min(m, n))); 00123 00124 cmat B(A); 00125 00126 zgesvd_(&jobu, &jobvt, &m, &n, B._data(), &lda, S._data(), U._data(), &ldu, V._data(), &ldvt, work._data(), &lwork, rwork._data(), &info); 00127 00128 V = V.H(); // This is slow!!! 00129 00130 return (info==0); 00131 } 00132 00133 #else 00134 00135 bool svd(const mat &A, vec &S) 00136 { 00137 it_error("LAPACK library is needed to use svd() function"); 00138 return false; 00139 } 00140 00141 bool svd(const cmat &A, vec &S) 00142 { 00143 it_error("LAPACK library is needed to use svd() function"); 00144 return false; 00145 } 00146 00147 bool svd(const mat &A, mat &U, vec &S, mat &V) 00148 { 00149 it_error("LAPACK library is needed to use svd() function"); 00150 return false; 00151 } 00152 00153 bool svd(const cmat &A, cmat &U, vec &S, cmat &V) 00154 { 00155 it_error("LAPACK library is needed to use svd() function"); 00156 return false; 00157 } 00158 00159 #endif // HAVE_LAPACK 00160 00161 vec svd(const mat &A) 00162 { 00163 vec S; 00164 svd(A, S); 00165 return S; 00166 } 00167 00168 vec svd(const cmat &A) 00169 { 00170 vec S; 00171 svd(A, S); 00172 return S; 00173 } 00174 00175 } //namespace itpp
Generated on Sun Dec 9 17:30:58 2007 for IT++ by Doxygen 1.5.4