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 #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 &Q, mat &R, bmat &P)
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     ivec jpvt(n);
00101     jpvt.zeros();
00102 
00103     R = A;
00104 
00105     // perform workspace query for optimum lwork value
00106     int lwork_tmp = -1;
00107     dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00108             &lwork_tmp, &info);
00109     if (info == 0) {
00110       lwork = static_cast<int>(work(0));
00111       work.set_size(lwork, false);
00112     }
00113     dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00114             &lwork, &info);
00115     Q = R;
00116     Q.set_size(m, m, true);
00117 
00118     // construct permutation matrix
00119     P = zeros_b(n, n);
00120     for (int j = 0; j < n; j++)
00121       P(jpvt(j)-1, j) = 1;
00122 
00123     // construct R
00124     for (int i = 0; i < m; i++)
00125       for (int j = 0; j < std::min(i, n); j++)
00126         R(i, j) = 0;
00127 
00128     // perform workspace query for optimum lwork value
00129     lwork_tmp = -1;
00130     dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00131             &info);
00132     if (info == 0) {
00133       lwork = static_cast<int>(work(0));
00134       work.set_size(lwork, false);
00135     }
00136     dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00137             &info);
00138 
00139     return (info==0);
00140   }
00141 
00142 
00143 
00144   bool qr(const cmat &A, cmat &Q, cmat &R)
00145   {
00146     int info;
00147     int m = A.rows();
00148     int n = A.cols();
00149     int lwork = n;
00150     int k = std::min(m, n);
00151     cvec tau(k);
00152     cvec work(lwork);
00153 
00154     R = A;
00155 
00156     // perform workspace query for optimum lwork value
00157     int lwork_tmp = -1;
00158     zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork_tmp,
00159             &info);
00160     if (info == 0) {
00161       lwork = static_cast<int>(real(work(0)));
00162       work.set_size(lwork, false);
00163     }
00164     zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
00165 
00166     Q = R;
00167     Q.set_size(m, m, true);
00168 
00169     // construct R
00170     for (int i = 0; i < m; i++)
00171       for (int j = 0; j < std::min(i, n); j++)
00172         R(i, j) = 0;
00173 
00174     // perform workspace query for optimum lwork value
00175     lwork_tmp = -1;
00176     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00177             &info);
00178     if (info == 0) {
00179       lwork = static_cast<int>(real(work(0)));
00180       work.set_size(lwork, false);
00181     }
00182     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00183             &info);
00184 
00185     return (info==0);
00186   }
00187 
00188   bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00189   {
00190     int info;
00191     int m = A.rows();
00192     int n = A.cols();
00193     int lwork = n;
00194     int k = std::min(m, n);
00195     cvec tau(k);
00196     cvec work(lwork);
00197     vec rwork(std::max(1, 2*n));
00198     ivec jpvt(n);
00199     jpvt.zeros();
00200 
00201     R = A;
00202 
00203     // perform workspace query for optimum lwork value
00204     int lwork_tmp = -1;
00205     zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00206             &lwork_tmp, rwork._data(), &info);
00207     if (info == 0) {
00208       lwork = static_cast<int>(real(work(0)));
00209       work.set_size(lwork, false);
00210     }
00211     zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(),
00212             &lwork, rwork._data(), &info);
00213 
00214     Q = R;
00215     Q.set_size(m, m, true);
00216 
00217     // construct permutation matrix
00218     P = zeros_b(n, n);
00219     for (int j = 0; j < n; j++)
00220       P(jpvt(j)-1, j) = 1;
00221 
00222     // construct R
00223     for (int i = 0; i < m; i++)
00224       for (int j = 0; j < std::min(i, n); j++)
00225         R(i, j) = 0;
00226 
00227     // perform workspace query for optimum lwork value
00228     lwork_tmp = -1;
00229     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork_tmp,
00230             &info);
00231     if (info == 0) {
00232       lwork = static_cast<int>(real(work(0)));
00233       work.set_size(lwork, false);
00234     }
00235     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork,
00236             &info);
00237 
00238     return (info==0);
00239   }
00240 
00241 #else
00242 
00243   bool qr(const mat &A, mat &Q, mat &R)
00244   {
00245     it_error("LAPACK library is needed to use qr() function");
00246     return false;
00247   }
00248 
00249   bool qr(const mat &A, mat &Q, mat &R, bmat &P)
00250   {
00251     it_error("LAPACK library is needed to use qr() function");
00252     return false;
00253   }
00254 
00255   bool qr(const cmat &A, cmat &Q, cmat &R)
00256   {
00257     it_error("LAPACK library is needed to use qr() function");
00258     return false;
00259   }
00260 
00261   bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00262   {
00263     it_error("LAPACK library is needed to use qr() function");
00264     return false;
00265   }
00266 
00267 #endif // HAVE_LAPACK
00268 
00269 } // namespace itpp
SourceForge Logo

Generated on Sat Apr 19 10:43:50 2008 for IT++ by Doxygen 1.5.5