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
Generated on Sun Sep 14 18:52:25 2008 for IT++ by Doxygen 1.5.6