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 
00042 
00043 namespace itpp {
00044 
00045 #if defined(HAVE_LAPACK)
00046 
00047   bool qr(const mat &A, mat &Q, mat &R)
00048   {
00049     int m, n, k, info, lwork, i, j;
00050 
00051     m = A.rows(); n = A.cols();
00052     lwork = 16*n;
00053     k = std::min(m,n);
00054     vec tau(k);
00055     vec work(lwork);
00056 
00057     R = A;
00058     dgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
00059     Q = R;
00060 
00061     // construct R
00062     for (i=0; i<m; i++)
00063       for (j=0; j<std::min(i,n); j++)
00064   R(i,j) = 0;
00065 
00066     Q.set_size(m, m, true);
00067     dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, &info);
00068 
00069     return (info==0);
00070   }
00071 
00072   bool qr(const mat &A, mat &Q, mat &R, bmat &P)
00073   {
00074     int m, n, k, info, lwork, i, j;
00075 
00076     m = A.rows(); n = A.cols();
00077     lwork = 16*n;
00078     k = std::min(m,n);
00079     vec tau(k);
00080     vec work(lwork);
00081     ivec jpvt(n); jpvt.zeros();
00082 
00083     R = A;
00084     P.set_size(n, n, false); P.zeros();
00085     dgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), &lwork, &info);
00086     Q = R;
00087 
00088     // construct permutation matrix
00089     for (j=0; j<n; j++)
00090       P(jpvt(j)-1, j) = 1;
00091 
00092     // construct R
00093     for (i=0; i<m; i++)
00094       for (j=0; j<std::min(i,n); j++)
00095   R(i,j) = 0;
00096 
00097     Q.set_size(m, m, true);
00098     dorgqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, &info);
00099 
00100     return (info==0);
00101   }
00102 
00103 
00104 
00105   bool qr(const cmat &A, cmat &Q, cmat &R)
00106   {
00107     int m, n, k, info, lwork, i, j;
00108 
00109     m = A.rows(); n = A.cols();
00110     lwork = 16*n;
00111     k = std::min(m,n);
00112     cvec tau(k);
00113     cvec work(lwork);
00114 
00115     R = A;
00116 
00117     zgeqrf_(&m, &n, R._data(), &m, tau._data(), work._data(), &lwork, &info);
00118     Q = R;
00119 
00120     // construct R
00121     for (i=0; i<m; i++)
00122       for (j=0; j<std::min(i,n); j++)
00123   R(i,j) = 0;
00124 
00125 
00126     Q.set_size(m, m, true);
00127     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, &info);
00128 
00129     return (info==0);
00130   }
00131 
00132   bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00133   {
00134     int m, n, k, info, lwork, i, j;
00135 
00136     m = A.rows(); n = A.cols();
00137     lwork = 16*n;
00138     k = std::min(m,n);
00139     cvec tau(k);
00140     cvec work(lwork);
00141     vec rwork(std::max(1, 2*n));
00142     ivec jpvt(n);
00143     jpvt.zeros();
00144 
00145     R = A;
00146     P.set_size(n, n, false); P.zeros();
00147 
00148     zgeqp3_(&m, &n, R._data(), &m, jpvt._data(), tau._data(), work._data(), &lwork, rwork._data(), &info);
00149     Q = R;
00150 
00151     // construct permutation matrix
00152     for (j=0; j<n; j++)
00153       P(jpvt(j)-1, j) = 1;
00154 
00155     // construct R
00156     for (i=0; i<m; i++)
00157       for (j=0; j<std::min(i,n); j++)
00158   R(i,j) = 0;
00159 
00160 
00161     Q.set_size(m, m, true);
00162     zungqr_(&m, &m, &k, Q._data(), &m, tau._data(), work._data(), &lwork, &info);
00163 
00164     return (info==0);
00165   }
00166 
00167 #else
00168 
00169   bool qr(const mat &A, mat &Q, mat &R)
00170   {
00171     it_error("LAPACK library is needed to use qr() function");
00172     return false;
00173   }
00174 
00175   bool qr(const mat &A, mat &Q, mat &R, bmat &P)
00176   {
00177     it_error("LAPACK library is needed to use qr() function");
00178     return false;
00179   }
00180 
00181   bool qr(const cmat &A, cmat &Q, cmat &R)
00182   {
00183     it_error("LAPACK library is needed to use qr() function");
00184     return false;
00185   }
00186 
00187   bool qr(const cmat &A, cmat &Q, cmat &R, bmat &P)
00188   {
00189     it_error("LAPACK library is needed to use qr() function");
00190     return false;
00191   }
00192 
00193 #endif // HAVE_LAPACK
00194 
00195 } // namespace itpp
SourceForge Logo

Generated on Sun Sep 14 18:52:33 2008 for IT++ by Doxygen 1.5.6