IT++ Logo
qr.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:29 for IT++ by Doxygen 1.7.4