IT++ Logo

qr.cpp

Go to the documentation of this file.
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/qr.h>
00041 #include <itpp/base/specmat.h>
00042 
00043 
00044 namespace itpp
00045 {
00046 
00047 #if defined(HAVE_LAPACK)
00048 
00049 bool qr(const mat &A, mat &Q, mat &R)
00050 {
00051   int info;
00052   int m = A.rows();
00053   int n = A.cols();
00054   int lwork = n;
00055   int k = std::min(m, n);
00056   vec tau(k);
00057   vec work(lwork);
00058 
00059   R = A;
00060 
00061   // perform workspace query for optimum lwork value
00062   int lwork_tmp = -1;
00063   dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
00064           &info);
00065   if (info == 0) {
00066     lwork = static_cast<int>(work(0));
00067     work.set_size(lwork, false);
00068   }
00069   dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
00070   Q = R;
00071   Q.set_size(m, m, true);
00072 
00073   // construct R
00074   for (int i = 0; i < m; i++)
00075     for (int j = 0; j < std::min(i, n); j++)
00076       R(i, j) = 0;
00077 
00078   // perform workspace query for optimum lwork value
00079   lwork_tmp = -1;
00080   dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00081           &info);
00082   if (info == 0) {
00083     lwork = static_cast<int>(work(0));
00084     work.set_size(lwork, false);
00085   }
00086   dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00087           &info);
00088 
00089   return (info == 0);
00090 }
00091 
00092 bool qr(const mat &A, mat &Q, mat &R, bmat &P)
00093 {
00094   int info;
00095   int m = A.rows();
00096   int n = A.cols();
00097   int lwork = n;
00098   int k = std::min(m, n);
00099   vec tau(k);
00100   vec work(lwork);
00101   ivec jpvt(n);
00102   jpvt.zeros();
00103 
00104   R = A;
00105 
00106   // perform workspace query for optimum lwork value
00107   int lwork_tmp = -1;
00108   dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00109           &lwork_tmp, &info);
00110   if (info == 0) {
00111     lwork = static_cast<int>(work(0));
00112     work.set_size(lwork, false);
00113   }
00114   dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00115           &lwork, &info);
00116   Q = R;
00117   Q.set_size(m, m, true);
00118 
00119   // construct permutation matrix
00120   P = zeros_b(n, n);
00121   for (int j = 0; j < n; j++)
00122     P(jpvt(j) - 1, j) = 1;
00123 
00124   // construct R
00125   for (int i = 0; i < m; i++)
00126     for (int j = 0; j < std::min(i, n); j++)
00127       R(i, j) = 0;
00128 
00129   // perform workspace query for optimum lwork value
00130   lwork_tmp = -1;
00131   dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00132           &info);
00133   if (info == 0) {
00134     lwork = static_cast<int>(work(0));
00135     work.set_size(lwork, false);
00136   }
00137   dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00138           &info);
00139 
00140   return (info == 0);
00141 }
00142 
00143 
00144 
00145 bool qr(const cmat &A, cmat &Q, cmat &R)
00146 {
00147   int info;
00148   int m = A.rows();
00149   int n = A.cols();
00150   int lwork = n;
00151   int k = std::min(m, n);
00152   cvec tau(k);
00153   cvec work(lwork);
00154 
00155   R = A;
00156 
00157   // perform workspace query for optimum lwork value
00158   int lwork_tmp = -1;
00159   zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
00160           &info);
00161   if (info == 0) {
00162     lwork = static_cast<int>(real(work(0)));
00163     work.set_size(lwork, false);
00164   }
00165   zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
00166 
00167   Q = R;
00168   Q.set_size(m, m, true);
00169 
00170   // construct R
00171   for (int i = 0; i < m; i++)
00172     for (int j = 0; j < std::min(i, n); j++)
00173       R(i, j) = 0;
00174 
00175   // perform workspace query for optimum lwork value
00176   lwork_tmp = -1;
00177   zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00178           &info);
00179   if (info == 0) {
00180     lwork = static_cast<int>(real(work(0)));
00181     work.set_size(lwork, false);
00182   }
00183   zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00184           &info);
00185 
00186   return (info == 0);
00187 }
00188 
00189 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00190 {
00191   int info;
00192   int m = A.rows();
00193   int n = A.cols();
00194   int lwork = n;
00195   int k = std::min(m, n);
00196   cvec tau(k);
00197   cvec work(lwork);
00198   vec rwork(std::max(1, 2*n));
00199   ivec jpvt(n);
00200   jpvt.zeros();
00201 
00202   R = A;
00203 
00204   // perform workspace query for optimum lwork value
00205   int lwork_tmp = -1;
00206   zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00207           &lwork_tmp, rwork._data(), &info);
00208   if (info == 0) {
00209     lwork = static_cast<int>(real(work(0)));
00210     work.set_size(lwork, false);
00211   }
00212   zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00213           &lwork, rwork._data(), &info);
00214 
00215   Q = R;
00216   Q.set_size(m, m, true);
00217 
00218   // construct permutation matrix
00219   P = zeros_b(n, n);
00220   for (int j = 0; j < n; j++)
00221     P(jpvt(j) - 1, j) = 1;
00222 
00223   // construct R
00224   for (int i = 0; i < m; i++)
00225     for (int j = 0; j < std::min(i, n); j++)
00226       R(i, j) = 0;
00227 
00228   // perform workspace query for optimum lwork value
00229   lwork_tmp = -1;
00230   zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00231           &info);
00232   if (info == 0) {
00233     lwork = static_cast<int>(real(work(0)));
00234     work.set_size(lwork, false);
00235   }
00236   zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00237           &info);
00238 
00239   return (info == 0);
00240 }
00241 
00242 #else
00243 
00244 bool qr(const mat &A, mat &Q, mat &R)
00245 {
00246   it_error("LAPACK library is needed to use qr() function");
00247   return false;
00248 }
00249 
00250 bool qr(const mat &A, mat &Q, mat &R, bmat &P)
00251 {
00252   it_error("LAPACK library is needed to use qr() function");
00253   return false;
00254 }
00255 
00256 bool qr(const cmat &A, cmat &Q, cmat &R)
00257 {
00258   it_error("LAPACK library is needed to use qr() function");
00259   return false;
00260 }
00261 
00262 bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00263 {
00264   it_error("LAPACK library is needed to use qr() function");
00265   return false;
00266 }
00267 
00268 #endif // HAVE_LAPACK
00269 
00270 } // namespace itpp
SourceForge Logo

Generated on Sun Jul 26 08:54:53 2009 for IT++ by Doxygen 1.5.9