IT++ Logo Newcom Logo

qr.cpp

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

Generated on Fri Jun 8 00:27:14 2007 for IT++ by Doxygen 1.5.2