00001 00029 #ifndef _MSC_VER 00030 # include <itpp/config.h> 00031 #else 00032 # include <itpp/config_msvc.h> 00033 #endif 00034 00035 #if defined(HAVE_LAPACK) 00036 # include <itpp/base/algebra/lapack.h> 00037 #endif 00038 00039 #include <itpp/base/algebra/qr.h> 00040 #include <itpp/base/specmat.h> 00041 00042 00043 namespace itpp 00044 { 00045 00046 #if defined(HAVE_LAPACK) 00047 00048 bool qr(const mat &A, mat &Q, mat &R) 00049 { 00050 int info; 00051 int m = A.rows(); 00052 int n = A.cols(); 00053 int lwork = n; 00054 int k = std::min(m, n); 00055 vec tau(k); 00056 vec work(lwork); 00057 00058 R = A; 00059 00060 // perform workspace query for optimum lwork value 00061 int lwork_tmp = -1; 00062 dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 00063 &info); 00064 if (info == 0) { 00065 lwork = static_cast<int>(work(0)); 00066 work.set_size(lwork, false); 00067 } 00068 dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); 00069 Q = R; 00070 Q.set_size(m, m, true); 00071 00072 // construct R 00073 for (int i = 0; i < m; i++) 00074 for (int j = 0; j < std::min(i, n); j++) 00075 R(i, j) = 0; 00076 00077 // perform workspace query for optimum lwork value 00078 lwork_tmp = -1; 00079 dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp, 00080 &info); 00081 if (info == 0) { 00082 lwork = static_cast<int>(work(0)); 00083 work.set_size(lwork, false); 00084 } 00085 dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, 00086 &info); 00087 00088 return (info == 0); 00089 } 00090 00091 bool qr(const mat &A, mat &R) 00092 { 00093 int info; 00094 int m = A.rows(); 00095 int n = A.cols(); 00096 int lwork = n; 00097 int k = std::min(m, n); 00098 vec tau(k); 00099 vec work(lwork); 00100 00101 R = A; 00102 00103 // perform workspace query for optimum lwork value 00104 int lwork_tmp = -1; 00105 dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 00106 &info); 00107 if (info == 0) { 00108 lwork = static_cast<int>(work(0)); 00109 work.set_size(lwork, false); 00110 } 00111 dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); 00112 00113 // construct R 00114 for (int i = 0; i < m; i++) 00115 for (int j = 0; j < std::min(i, n); j++) 00116 R(i, j) = 0; 00117 00118 return (info == 0); 00119 } 00120 00121 bool qr(const mat &A, mat &Q, mat &R, bmat &P) 00122 { 00123 int info; 00124 int m = A.rows(); 00125 int n = A.cols(); 00126 int lwork = n; 00127 int k = std::min(m, n); 00128 vec tau(k); 00129 vec work(lwork); 00130 ivec jpvt(n); 00131 jpvt.zeros(); 00132 00133 R = A; 00134 00135 // perform workspace query for optimum lwork value 00136 int lwork_tmp = -1; 00137 dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), 00138 &lwork_tmp, &info); 00139 if (info == 0) { 00140 lwork = static_cast<int>(work(0)); 00141 work.set_size(lwork, false); 00142 } 00143 dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), 00144 &lwork, &info); 00145 Q = R; 00146 Q.set_size(m, m, true); 00147 00148 // construct permutation matrix 00149 P = zeros_b(n, n); 00150 for (int j = 0; j < n; j++) 00151 P(jpvt(j) - 1, j) = 1; 00152 00153 // construct R 00154 for (int i = 0; i < m; i++) 00155 for (int j = 0; j < std::min(i, n); j++) 00156 R(i, j) = 0; 00157 00158 // perform workspace query for optimum lwork value 00159 lwork_tmp = -1; 00160 dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp, 00161 &info); 00162 if (info == 0) { 00163 lwork = static_cast<int>(work(0)); 00164 work.set_size(lwork, false); 00165 } 00166 dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, 00167 &info); 00168 00169 return (info == 0); 00170 } 00171 00172 00173 00174 bool qr(const cmat &A, cmat &Q, cmat &R) 00175 { 00176 int info; 00177 int m = A.rows(); 00178 int n = A.cols(); 00179 int lwork = n; 00180 int k = std::min(m, n); 00181 cvec tau(k); 00182 cvec work(lwork); 00183 00184 R = A; 00185 00186 // perform workspace query for optimum lwork value 00187 int lwork_tmp = -1; 00188 zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 00189 &info); 00190 if (info == 0) { 00191 lwork = static_cast<int>(real(work(0))); 00192 work.set_size(lwork, false); 00193 } 00194 zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); 00195 00196 Q = R; 00197 Q.set_size(m, m, true); 00198 00199 // construct R 00200 for (int i = 0; i < m; i++) 00201 for (int j = 0; j < std::min(i, n); j++) 00202 R(i, j) = 0; 00203 00204 // perform workspace query for optimum lwork value 00205 lwork_tmp = -1; 00206 zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp, 00207 &info); 00208 if (info == 0) { 00209 lwork = static_cast<int>(real(work(0))); 00210 work.set_size(lwork, false); 00211 } 00212 zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, 00213 &info); 00214 00215 return (info == 0); 00216 } 00217 00218 bool qr(const cmat &A, cmat &R) 00219 { 00220 int info; 00221 int m = A.rows(); 00222 int n = A.cols(); 00223 int lwork = n; 00224 int k = std::min(m, n); 00225 cvec tau(k); 00226 cvec work(lwork); 00227 00228 R = A; 00229 00230 // perform workspace query for optimum lwork value 00231 int lwork_tmp = -1; 00232 zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp, 00233 &info); 00234 if (info == 0) { 00235 lwork = static_cast<int>(real(work(0))); 00236 work.set_size(lwork, false); 00237 } 00238 zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info); 00239 00240 // construct R 00241 for (int i = 0; i < m; i++) 00242 for (int j = 0; j < std::min(i, n); j++) 00243 R(i, j) = 0; 00244 00245 return (info == 0); 00246 } 00247 00248 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P) 00249 { 00250 int info; 00251 int m = A.rows(); 00252 int n = A.cols(); 00253 int lwork = n; 00254 int k = std::min(m, n); 00255 cvec tau(k); 00256 cvec work(lwork); 00257 vec rwork(std::max(1, 2*n)); 00258 ivec jpvt(n); 00259 jpvt.zeros(); 00260 00261 R = A; 00262 00263 // perform workspace query for optimum lwork value 00264 int lwork_tmp = -1; 00265 zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), 00266 &lwork_tmp, rwork._data(), &info); 00267 if (info == 0) { 00268 lwork = static_cast<int>(real(work(0))); 00269 work.set_size(lwork, false); 00270 } 00271 zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), 00272 &lwork, rwork._data(), &info); 00273 00274 Q = R; 00275 Q.set_size(m, m, true); 00276 00277 // construct permutation matrix 00278 P = zeros_b(n, n); 00279 for (int j = 0; j < n; j++) 00280 P(jpvt(j) - 1, j) = 1; 00281 00282 // construct R 00283 for (int i = 0; i < m; i++) 00284 for (int j = 0; j < std::min(i, n); j++) 00285 R(i, j) = 0; 00286 00287 // perform workspace query for optimum lwork value 00288 lwork_tmp = -1; 00289 zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp, 00290 &info); 00291 if (info == 0) { 00292 lwork = static_cast<int>(real(work(0))); 00293 work.set_size(lwork, false); 00294 } 00295 zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, 00296 &info); 00297 00298 return (info == 0); 00299 } 00300 00301 #else 00302 00303 bool qr(const mat &A, mat &Q, mat &R) 00304 { 00305 it_error("LAPACK library is needed to use qr() function"); 00306 return false; 00307 } 00308 00309 bool qr(const mat &A, mat &R) 00310 { 00311 it_error("LAPACK library is needed to use qr() function"); 00312 return false; 00313 } 00314 00315 bool qr(const mat &A, mat &Q, mat &R, bmat &P) 00316 { 00317 it_error("LAPACK library is needed to use qr() function"); 00318 return false; 00319 } 00320 00321 bool qr(const cmat &A, cmat &Q, cmat &R) 00322 { 00323 it_error("LAPACK library is needed to use qr() function"); 00324 return false; 00325 } 00326 00327 bool qr(const cmat &A, cmat &R) 00328 { 00329 it_error("LAPACK library is needed to use qr() function"); 00330 return false; 00331 } 00332 00333 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P) 00334 { 00335 it_error("LAPACK library is needed to use qr() function"); 00336 return false; 00337 } 00338 00339 #endif // HAVE_LAPACK 00340 00341 } // namespace itpp
Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4