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
Generated on Sat Apr 19 10:41:54 2008 for IT++ by Doxygen 1.5.5