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
Generated on Wed Apr 18 11:45:33 2007 for IT++ by Doxygen 1.5.2