Classes | |
class | auxlib |
wrapper for accessing external functions defined in ATLAS, LAPACK or BLAS libraries More... | |
Functions | |
template<typename eT > | |
static void | auxlib::inv_noalias (Mat< eT > &out, const Mat< eT > &X) |
immediate matrix inverse | |
template<typename eT > | |
static void | auxlib::inv_inplace (Mat< eT > &X) |
immediate inplace matrix inverse | |
template<typename eT > | |
static eT | auxlib::det (const Mat< eT > &X) |
immediate determinant of a matrix using ATLAS or LAPACK | |
template<typename eT > | |
static void | auxlib::lu (Mat< eT > &L, Mat< eT > &U, podarray< int > &ipiv, const Mat< eT > &X_orig) |
immediate LU decomposition of a matrix using ATLAS or LAPACK | |
template<typename eT > | |
static void | auxlib::lu (Mat< eT > &L, Mat< eT > &U, Mat< eT > &P, const Mat< eT > &X) |
template<typename eT > | |
static void | auxlib::lu (Mat< eT > &L, Mat< eT > &U, const Mat< eT > &X) |
template<typename eT > | |
static void | auxlib::eig_sym (Col< eT > &eigval, const Mat< eT > &A) |
immediate eigenvalues of a symmetric real matrix using LAPACK | |
template<typename T > | |
static void | auxlib::eig_sym (Col< T > &eigval, const Mat< std::complex< T > > &A) |
immediate eigenvalues of a hermitian complex matrix using LAPACK | |
template<typename eT > | |
static void | auxlib::eig_sym (Col< eT > &eigval, Mat< eT > &eigvec, const Mat< eT > &A) |
immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK | |
template<typename T > | |
static void | auxlib::eig_sym (Col< T > &eigval, Mat< std::complex< T > > &eigvec, const Mat< std::complex< T > > &A) |
immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK | |
template<typename T > | |
void | auxlib::eig_gen (Col< std::complex< T > > &eigval, Mat< T > &l_eigvec, Mat< T > &r_eigvec, const Mat< T > &A, const char side) |
Eigenvalues and eigenvectors of a general square real matrix using LAPACK. The argument 'side' specifies which eigenvectors should be calculated (see code for mode details). | |
template<typename T > | |
static void | auxlib::eig_gen (Col< std::complex< T > > &eigval, Mat< std::complex< T > > &l_eigvec, Mat< std::complex< T > > &r_eigvec, const Mat< std::complex< T > > &A, const char side) |
Eigenvalues and eigenvectors of a general square complex matrix using LAPACK The argument 'side' specifies which eigenvectors should be calculated (see code for mode details). | |
template<typename eT > | |
static bool | auxlib::chol (Mat< eT > &out, const Mat< eT > &X) |
template<typename eT > | |
static bool | auxlib::qr (Mat< eT > &Q, Mat< eT > &R, const Mat< eT > &X) |
template<typename eT > | |
static bool | auxlib::svd (Col< eT > &S, const Mat< eT > &X) |
template<typename T > | |
static bool | auxlib::svd (Col< T > &S, const Mat< std::complex< T > > &X) |
template<typename eT > | |
static bool | auxlib::svd (Mat< eT > &U, Col< eT > &S, Mat< eT > &V, const Mat< eT > &X) |
template<typename T > | |
static bool | auxlib::svd (Mat< std::complex< T > > &U, Col< T > &S, Mat< std::complex< T > > &V, const Mat< std::complex< T > > &X) |
template<typename eT > | |
static bool | auxlib::solve (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B) |
Solve a system of linear equations Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows. | |
template<typename eT > | |
static bool | auxlib::solve_od (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B) |
Solve an over-determined system. Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows. | |
template<typename eT > | |
static bool | auxlib::solve_ud (Mat< eT > &out, const Mat< eT > &A, const Mat< eT > &B) |
Solve an under-determined system. Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows. |
void auxlib::inv_noalias | ( | Mat< eT > & | out, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
immediate matrix inverse
Definition at line 25 of file auxlib_meat.hpp.
References arma_check(), arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), Mat< eT >::set_size(), and tmp_real().
Referenced by op_inv::apply().
00026 { 00027 arma_extra_debug_sigprint(); 00028 00029 switch(X.n_rows) 00030 { 00031 case 1: 00032 { 00033 out.set_size(1,1); 00034 out[0] = eT(1) / X[0]; 00035 }; 00036 break; 00037 00038 case 2: 00039 { 00040 out.set_size(2,2); 00041 00042 const eT a = X.at(0,0); 00043 const eT b = X.at(0,1); 00044 const eT c = X.at(1,0); 00045 const eT d = X.at(1,1); 00046 00047 const eT k = eT(1) / (a*d - b*c); 00048 00049 out.at(0,0) = d*k; 00050 out.at(0,1) = -b*k; 00051 out.at(1,0) = -c*k; 00052 out.at(1,1) = a*k; 00053 }; 00054 break; 00055 00056 case 3: 00057 { 00058 out.set_size(3,3); 00059 00060 const eT* X_col0 = X.colptr(0); 00061 const eT a11 = X_col0[0]; 00062 const eT a21 = X_col0[1]; 00063 const eT a31 = X_col0[2]; 00064 00065 const eT* X_col1 = X.colptr(1); 00066 const eT a12 = X_col1[0]; 00067 const eT a22 = X_col1[1]; 00068 const eT a32 = X_col1[2]; 00069 00070 const eT* X_col2 = X.colptr(2); 00071 const eT a13 = X_col2[0]; 00072 const eT a23 = X_col2[1]; 00073 const eT a33 = X_col2[2]; 00074 00075 const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); 00076 00077 00078 eT* out_col0 = out.colptr(0); 00079 out_col0[0] = (a33*a22 - a32*a23) * k; 00080 out_col0[1] = -(a33*a21 - a31*a23) * k; 00081 out_col0[2] = (a32*a21 - a31*a22) * k; 00082 00083 eT* out_col1 = out.colptr(1); 00084 out_col1[0] = -(a33*a12 - a32*a13) * k; 00085 out_col1[1] = (a33*a11 - a31*a13) * k; 00086 out_col1[2] = -(a32*a11 - a31*a12) * k; 00087 00088 eT* out_col2 = out.colptr(2); 00089 out_col2[0] = (a23*a12 - a22*a13) * k; 00090 out_col2[1] = -(a23*a11 - a21*a13) * k; 00091 out_col2[2] = (a22*a11 - a21*a12) * k; 00092 }; 00093 break; 00094 00095 case 4: 00096 { 00097 out.set_size(4,4); 00098 00099 out.at(0,0) = X.at(1,2)*X.at(2,3)*X.at(3,1) - X.at(1,3)*X.at(2,2)*X.at(3,1) + X.at(1,3)*X.at(2,1)*X.at(3,2) - X.at(1,1)*X.at(2,3)*X.at(3,2) - X.at(1,2)*X.at(2,1)*X.at(3,3) + X.at(1,1)*X.at(2,2)*X.at(3,3); 00100 out.at(1,0) = X.at(1,3)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,0)*X.at(3,2) + X.at(1,0)*X.at(2,3)*X.at(3,2) + X.at(1,2)*X.at(2,0)*X.at(3,3) - X.at(1,0)*X.at(2,2)*X.at(3,3); 00101 out.at(2,0) = X.at(1,1)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,1)*X.at(3,0) + X.at(1,3)*X.at(2,0)*X.at(3,1) - X.at(1,0)*X.at(2,3)*X.at(3,1) - X.at(1,1)*X.at(2,0)*X.at(3,3) + X.at(1,0)*X.at(2,1)*X.at(3,3); 00102 out.at(3,0) = X.at(1,2)*X.at(2,1)*X.at(3,0) - X.at(1,1)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,0)*X.at(3,1) + X.at(1,0)*X.at(2,2)*X.at(3,1) + X.at(1,1)*X.at(2,0)*X.at(3,2) - X.at(1,0)*X.at(2,1)*X.at(3,2); 00103 00104 out.at(0,1) = X.at(0,3)*X.at(2,2)*X.at(3,1) - X.at(0,2)*X.at(2,3)*X.at(3,1) - X.at(0,3)*X.at(2,1)*X.at(3,2) + X.at(0,1)*X.at(2,3)*X.at(3,2) + X.at(0,2)*X.at(2,1)*X.at(3,3) - X.at(0,1)*X.at(2,2)*X.at(3,3); 00105 out.at(1,1) = X.at(0,2)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,2)*X.at(3,0) + X.at(0,3)*X.at(2,0)*X.at(3,2) - X.at(0,0)*X.at(2,3)*X.at(3,2) - X.at(0,2)*X.at(2,0)*X.at(3,3) + X.at(0,0)*X.at(2,2)*X.at(3,3); 00106 out.at(2,1) = X.at(0,3)*X.at(2,1)*X.at(3,0) - X.at(0,1)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,0)*X.at(3,1) + X.at(0,0)*X.at(2,3)*X.at(3,1) + X.at(0,1)*X.at(2,0)*X.at(3,3) - X.at(0,0)*X.at(2,1)*X.at(3,3); 00107 out.at(3,1) = X.at(0,1)*X.at(2,2)*X.at(3,0) - X.at(0,2)*X.at(2,1)*X.at(3,0) + X.at(0,2)*X.at(2,0)*X.at(3,1) - X.at(0,0)*X.at(2,2)*X.at(3,1) - X.at(0,1)*X.at(2,0)*X.at(3,2) + X.at(0,0)*X.at(2,1)*X.at(3,2); 00108 00109 out.at(0,2) = X.at(0,2)*X.at(1,3)*X.at(3,1) - X.at(0,3)*X.at(1,2)*X.at(3,1) + X.at(0,3)*X.at(1,1)*X.at(3,2) - X.at(0,1)*X.at(1,3)*X.at(3,2) - X.at(0,2)*X.at(1,1)*X.at(3,3) + X.at(0,1)*X.at(1,2)*X.at(3,3); 00110 out.at(1,2) = X.at(0,3)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,0)*X.at(3,2) + X.at(0,0)*X.at(1,3)*X.at(3,2) + X.at(0,2)*X.at(1,0)*X.at(3,3) - X.at(0,0)*X.at(1,2)*X.at(3,3); 00111 out.at(2,2) = X.at(0,1)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,1)*X.at(3,0) + X.at(0,3)*X.at(1,0)*X.at(3,1) - X.at(0,0)*X.at(1,3)*X.at(3,1) - X.at(0,1)*X.at(1,0)*X.at(3,3) + X.at(0,0)*X.at(1,1)*X.at(3,3); 00112 out.at(3,2) = X.at(0,2)*X.at(1,1)*X.at(3,0) - X.at(0,1)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,0)*X.at(3,1) + X.at(0,0)*X.at(1,2)*X.at(3,1) + X.at(0,1)*X.at(1,0)*X.at(3,2) - X.at(0,0)*X.at(1,1)*X.at(3,2); 00113 00114 out.at(0,3) = X.at(0,3)*X.at(1,2)*X.at(2,1) - X.at(0,2)*X.at(1,3)*X.at(2,1) - X.at(0,3)*X.at(1,1)*X.at(2,2) + X.at(0,1)*X.at(1,3)*X.at(2,2) + X.at(0,2)*X.at(1,1)*X.at(2,3) - X.at(0,1)*X.at(1,2)*X.at(2,3); 00115 out.at(1,3) = X.at(0,2)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,2)*X.at(2,0) + X.at(0,3)*X.at(1,0)*X.at(2,2) - X.at(0,0)*X.at(1,3)*X.at(2,2) - X.at(0,2)*X.at(1,0)*X.at(2,3) + X.at(0,0)*X.at(1,2)*X.at(2,3); 00116 out.at(2,3) = X.at(0,3)*X.at(1,1)*X.at(2,0) - X.at(0,1)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,0)*X.at(2,1) + X.at(0,0)*X.at(1,3)*X.at(2,1) + X.at(0,1)*X.at(1,0)*X.at(2,3) - X.at(0,0)*X.at(1,1)*X.at(2,3); 00117 out.at(3,3) = X.at(0,1)*X.at(1,2)*X.at(2,0) - X.at(0,2)*X.at(1,1)*X.at(2,0) + X.at(0,2)*X.at(1,0)*X.at(2,1) - X.at(0,0)*X.at(1,2)*X.at(2,1) - X.at(0,1)*X.at(1,0)*X.at(2,2) + X.at(0,0)*X.at(1,1)*X.at(2,2); 00118 00119 out /= det(X); 00120 }; 00121 break; 00122 00123 default: 00124 { 00125 #if defined(ARMA_USE_ATLAS) 00126 { 00127 out = X; 00128 podarray<int> ipiv(out.n_rows); 00129 00130 int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr()); 00131 00132 if(info == 0) 00133 { 00134 info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr()); 00135 } 00136 00137 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" ); 00138 } 00139 #elif defined(ARMA_USE_LAPACK) 00140 { 00141 out = X; 00142 00143 int n_rows = out.n_rows; 00144 int n_cols = out.n_cols; 00145 int info = 0; 00146 00147 podarray<int> ipiv(out.n_rows); 00148 00149 // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6) 00150 // based on tests with various matrix types on 32-bit and 64-bit machines 00151 // 00152 // the "work" array is deliberately long so that a secondary (time-consuming) 00153 // memory allocation is avoided, if possible 00154 00155 int work_len = (std::max)(1, n_rows*84); 00156 podarray<eT> work(work_len); 00157 00158 lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info); 00159 00160 if(info == 0) 00161 { 00162 // query for optimum size of work_len 00163 00164 int work_len_tmp = -1; 00165 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info); 00166 00167 if(info == 0) 00168 { 00169 int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0])); 00170 00171 // if necessary, allocate more memory 00172 if(work_len < proposed_work_len) 00173 { 00174 work_len = proposed_work_len; 00175 work.set_size(work_len); 00176 } 00177 } 00178 00179 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info); 00180 } 00181 00182 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" ); 00183 } 00184 #else 00185 { 00186 arma_stop("auxlib::inv_noalias(): need ATLAS or LAPACK library"); 00187 } 00188 #endif 00189 }; 00190 } 00191 }
void auxlib::inv_inplace | ( | Mat< eT > & | X | ) | [inline, static, inherited] |
immediate inplace matrix inverse
Definition at line 199 of file auxlib_meat.hpp.
References arma_check(), arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), and tmp_real().
Referenced by op_inv::apply().
00200 { 00201 arma_extra_debug_sigprint(); 00202 00203 // for more info, see: 00204 // http://www.dr-lex.34sp.com/random/matrix_inv.html 00205 // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html 00206 // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm 00207 // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl 00208 00209 switch(X.n_rows) 00210 { 00211 case 1: 00212 { 00213 X[0] = eT(1) / X[0]; 00214 }; 00215 break; 00216 00217 case 2: 00218 { 00219 const eT a = X.at(0,0); 00220 const eT b = X.at(0,1); 00221 const eT c = X.at(1,0); 00222 const eT d = X.at(1,1); 00223 00224 const eT k = eT(1) / (a*d - b*c); 00225 00226 X.at(0,0) = d*k; 00227 X.at(0,1) = -b*k; 00228 X.at(1,0) = -c*k; 00229 X.at(1,1) = a*k; 00230 }; 00231 break; 00232 00233 case 3: 00234 { 00235 eT* X_col0 = X.colptr(0); 00236 eT* X_col1 = X.colptr(1); 00237 eT* X_col2 = X.colptr(2); 00238 00239 const eT a11 = X_col0[0]; 00240 const eT a21 = X_col0[1]; 00241 const eT a31 = X_col0[2]; 00242 00243 const eT a12 = X_col1[0]; 00244 const eT a22 = X_col1[1]; 00245 const eT a32 = X_col1[2]; 00246 00247 const eT a13 = X_col2[0]; 00248 const eT a23 = X_col2[1]; 00249 const eT a33 = X_col2[2]; 00250 00251 const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); 00252 00253 X_col0[0] = (a33*a22 - a32*a23) * k; 00254 X_col0[1] = -(a33*a21 - a31*a23) * k; 00255 X_col0[2] = (a32*a21 - a31*a22) * k; 00256 00257 X_col1[0] = -(a33*a12 - a32*a13) * k; 00258 X_col1[1] = (a33*a11 - a31*a13) * k; 00259 X_col1[2] = -(a32*a11 - a31*a12) * k; 00260 00261 X_col2[0] = (a23*a12 - a22*a13) * k; 00262 X_col2[1] = -(a23*a11 - a21*a13) * k; 00263 X_col2[2] = (a22*a11 - a21*a12) * k; 00264 }; 00265 break; 00266 00267 case 4: 00268 { 00269 const Mat<eT> A(X); 00270 00271 X.at(0,0) = A.at(1,2)*A.at(2,3)*A.at(3,1) - A.at(1,3)*A.at(2,2)*A.at(3,1) + A.at(1,3)*A.at(2,1)*A.at(3,2) - A.at(1,1)*A.at(2,3)*A.at(3,2) - A.at(1,2)*A.at(2,1)*A.at(3,3) + A.at(1,1)*A.at(2,2)*A.at(3,3); 00272 X.at(1,0) = A.at(1,3)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,0)*A.at(3,2) + A.at(1,0)*A.at(2,3)*A.at(3,2) + A.at(1,2)*A.at(2,0)*A.at(3,3) - A.at(1,0)*A.at(2,2)*A.at(3,3); 00273 X.at(2,0) = A.at(1,1)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,1)*A.at(3,0) + A.at(1,3)*A.at(2,0)*A.at(3,1) - A.at(1,0)*A.at(2,3)*A.at(3,1) - A.at(1,1)*A.at(2,0)*A.at(3,3) + A.at(1,0)*A.at(2,1)*A.at(3,3); 00274 X.at(3,0) = A.at(1,2)*A.at(2,1)*A.at(3,0) - A.at(1,1)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,0)*A.at(3,1) + A.at(1,0)*A.at(2,2)*A.at(3,1) + A.at(1,1)*A.at(2,0)*A.at(3,2) - A.at(1,0)*A.at(2,1)*A.at(3,2); 00275 00276 X.at(0,1) = A.at(0,3)*A.at(2,2)*A.at(3,1) - A.at(0,2)*A.at(2,3)*A.at(3,1) - A.at(0,3)*A.at(2,1)*A.at(3,2) + A.at(0,1)*A.at(2,3)*A.at(3,2) + A.at(0,2)*A.at(2,1)*A.at(3,3) - A.at(0,1)*A.at(2,2)*A.at(3,3); 00277 X.at(1,1) = A.at(0,2)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,2)*A.at(3,0) + A.at(0,3)*A.at(2,0)*A.at(3,2) - A.at(0,0)*A.at(2,3)*A.at(3,2) - A.at(0,2)*A.at(2,0)*A.at(3,3) + A.at(0,0)*A.at(2,2)*A.at(3,3); 00278 X.at(2,1) = A.at(0,3)*A.at(2,1)*A.at(3,0) - A.at(0,1)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,0)*A.at(3,1) + A.at(0,0)*A.at(2,3)*A.at(3,1) + A.at(0,1)*A.at(2,0)*A.at(3,3) - A.at(0,0)*A.at(2,1)*A.at(3,3); 00279 X.at(3,1) = A.at(0,1)*A.at(2,2)*A.at(3,0) - A.at(0,2)*A.at(2,1)*A.at(3,0) + A.at(0,2)*A.at(2,0)*A.at(3,1) - A.at(0,0)*A.at(2,2)*A.at(3,1) - A.at(0,1)*A.at(2,0)*A.at(3,2) + A.at(0,0)*A.at(2,1)*A.at(3,2); 00280 00281 X.at(0,2) = A.at(0,2)*A.at(1,3)*A.at(3,1) - A.at(0,3)*A.at(1,2)*A.at(3,1) + A.at(0,3)*A.at(1,1)*A.at(3,2) - A.at(0,1)*A.at(1,3)*A.at(3,2) - A.at(0,2)*A.at(1,1)*A.at(3,3) + A.at(0,1)*A.at(1,2)*A.at(3,3); 00282 X.at(1,2) = A.at(0,3)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,0)*A.at(3,2) + A.at(0,0)*A.at(1,3)*A.at(3,2) + A.at(0,2)*A.at(1,0)*A.at(3,3) - A.at(0,0)*A.at(1,2)*A.at(3,3); 00283 X.at(2,2) = A.at(0,1)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,1)*A.at(3,0) + A.at(0,3)*A.at(1,0)*A.at(3,1) - A.at(0,0)*A.at(1,3)*A.at(3,1) - A.at(0,1)*A.at(1,0)*A.at(3,3) + A.at(0,0)*A.at(1,1)*A.at(3,3); 00284 X.at(3,2) = A.at(0,2)*A.at(1,1)*A.at(3,0) - A.at(0,1)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,0)*A.at(3,1) + A.at(0,0)*A.at(1,2)*A.at(3,1) + A.at(0,1)*A.at(1,0)*A.at(3,2) - A.at(0,0)*A.at(1,1)*A.at(3,2); 00285 00286 X.at(0,3) = A.at(0,3)*A.at(1,2)*A.at(2,1) - A.at(0,2)*A.at(1,3)*A.at(2,1) - A.at(0,3)*A.at(1,1)*A.at(2,2) + A.at(0,1)*A.at(1,3)*A.at(2,2) + A.at(0,2)*A.at(1,1)*A.at(2,3) - A.at(0,1)*A.at(1,2)*A.at(2,3); 00287 X.at(1,3) = A.at(0,2)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,2)*A.at(2,0) + A.at(0,3)*A.at(1,0)*A.at(2,2) - A.at(0,0)*A.at(1,3)*A.at(2,2) - A.at(0,2)*A.at(1,0)*A.at(2,3) + A.at(0,0)*A.at(1,2)*A.at(2,3); 00288 X.at(2,3) = A.at(0,3)*A.at(1,1)*A.at(2,0) - A.at(0,1)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,0)*A.at(2,1) + A.at(0,0)*A.at(1,3)*A.at(2,1) + A.at(0,1)*A.at(1,0)*A.at(2,3) - A.at(0,0)*A.at(1,1)*A.at(2,3); 00289 X.at(3,3) = A.at(0,1)*A.at(1,2)*A.at(2,0) - A.at(0,2)*A.at(1,1)*A.at(2,0) + A.at(0,2)*A.at(1,0)*A.at(2,1) - A.at(0,0)*A.at(1,2)*A.at(2,1) - A.at(0,1)*A.at(1,0)*A.at(2,2) + A.at(0,0)*A.at(1,1)*A.at(2,2); 00290 00291 X /= det(A); 00292 }; 00293 break; 00294 00295 default: 00296 { 00297 #if defined(ARMA_USE_ATLAS) 00298 { 00299 Mat<eT>& out = X; 00300 podarray<int> ipiv(out.n_rows); 00301 00302 int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr()); 00303 00304 if(info == 0) 00305 { 00306 info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr()); 00307 } 00308 00309 arma_check( (info > 0), "auxlib::inv_inplace(): matrix appears to be singular" ); 00310 } 00311 #elif defined(ARMA_USE_LAPACK) 00312 { 00313 Mat<eT>& out = X; 00314 00315 int n_rows = out.n_rows; 00316 int n_cols = out.n_cols; 00317 int info = 0; 00318 00319 podarray<int> ipiv(out.n_rows); 00320 00321 // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6) 00322 // based on tests with various matrix types on 32-bit and 64-bit machines 00323 // 00324 // the "work" array is deliberately long so that a secondary (time-consuming) 00325 // memory allocation is avoided, if possible 00326 00327 int work_len = (std::max)(1, n_rows*84); 00328 podarray<eT> work(work_len); 00329 00330 lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info); 00331 00332 if(info == 0) 00333 { 00334 // query for optimum size of work_len 00335 00336 int work_len_tmp = -1; 00337 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info); 00338 00339 if(info == 0) 00340 { 00341 int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0])); 00342 00343 // if necessary, allocate more memory 00344 if(work_len < proposed_work_len) 00345 { 00346 work_len = proposed_work_len; 00347 work.set_size(work_len); 00348 } 00349 } 00350 00351 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info); 00352 } 00353 00354 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" ); 00355 } 00356 #else 00357 { 00358 arma_stop("auxlib::inv_inplace(): need ATLAS or LAPACK"); 00359 } 00360 #endif 00361 } 00362 00363 } 00364 00365 }
eT auxlib::det | ( | const Mat< eT > & | X | ) | [inline, static, inherited] |
immediate determinant of a matrix using ATLAS or LAPACK
Definition at line 372 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), Mat< eT >::colptr(), lapack::getrf_(), podarray< T1 >::mem, podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.
Referenced by det(), inv_inplace(), and inv_noalias().
00373 { 00374 arma_extra_debug_sigprint(); 00375 00376 switch(X.n_rows) 00377 { 00378 case 0: 00379 return 0.0; 00380 00381 case 1: 00382 return X[0]; 00383 00384 case 2: 00385 return (X.at(0,0)*X.at(1,1) - X.at(0,1)*X.at(1,0)); 00386 00387 case 3: 00388 { 00389 const eT* a_col0 = X.colptr(0); 00390 const eT a11 = a_col0[0]; 00391 const eT a21 = a_col0[1]; 00392 const eT a31 = a_col0[2]; 00393 00394 const eT* a_col1 = X.colptr(1); 00395 const eT a12 = a_col1[0]; 00396 const eT a22 = a_col1[1]; 00397 const eT a32 = a_col1[2]; 00398 00399 const eT* a_col2 = X.colptr(2); 00400 const eT a13 = a_col2[0]; 00401 const eT a23 = a_col2[1]; 00402 const eT a33 = a_col2[2]; 00403 00404 return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); 00405 00406 // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2); 00407 // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0); 00408 // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1); 00409 // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2); 00410 // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0); 00411 // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1); 00412 // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6); 00413 } 00414 00415 case 4: 00416 { 00417 const eT val = \ 00418 X.at(0,3) * X.at(1,2) * X.at(2,1) * X.at(3,0) \ 00419 - X.at(0,2) * X.at(1,3) * X.at(2,1) * X.at(3,0) \ 00420 - X.at(0,3) * X.at(1,1) * X.at(2,2) * X.at(3,0) \ 00421 + X.at(0,1) * X.at(1,3) * X.at(2,2) * X.at(3,0) \ 00422 + X.at(0,2) * X.at(1,1) * X.at(2,3) * X.at(3,0) \ 00423 - X.at(0,1) * X.at(1,2) * X.at(2,3) * X.at(3,0) \ 00424 - X.at(0,3) * X.at(1,2) * X.at(2,0) * X.at(3,1) \ 00425 + X.at(0,2) * X.at(1,3) * X.at(2,0) * X.at(3,1) \ 00426 + X.at(0,3) * X.at(1,0) * X.at(2,2) * X.at(3,1) \ 00427 - X.at(0,0) * X.at(1,3) * X.at(2,2) * X.at(3,1) \ 00428 - X.at(0,2) * X.at(1,0) * X.at(2,3) * X.at(3,1) \ 00429 + X.at(0,0) * X.at(1,2) * X.at(2,3) * X.at(3,1) \ 00430 + X.at(0,3) * X.at(1,1) * X.at(2,0) * X.at(3,2) \ 00431 - X.at(0,1) * X.at(1,3) * X.at(2,0) * X.at(3,2) \ 00432 - X.at(0,3) * X.at(1,0) * X.at(2,1) * X.at(3,2) \ 00433 + X.at(0,0) * X.at(1,3) * X.at(2,1) * X.at(3,2) \ 00434 + X.at(0,1) * X.at(1,0) * X.at(2,3) * X.at(3,2) \ 00435 - X.at(0,0) * X.at(1,1) * X.at(2,3) * X.at(3,2) \ 00436 - X.at(0,2) * X.at(1,1) * X.at(2,0) * X.at(3,3) \ 00437 + X.at(0,1) * X.at(1,2) * X.at(2,0) * X.at(3,3) \ 00438 + X.at(0,2) * X.at(1,0) * X.at(2,1) * X.at(3,3) \ 00439 - X.at(0,0) * X.at(1,2) * X.at(2,1) * X.at(3,3) \ 00440 - X.at(0,1) * X.at(1,0) * X.at(2,2) * X.at(3,3) \ 00441 + X.at(0,0) * X.at(1,1) * X.at(2,2) * X.at(3,3) \ 00442 ; 00443 00444 return val; 00445 } 00446 00447 default: 00448 { 00449 #if defined(ARMA_USE_ATLAS) 00450 { 00451 Mat<eT> tmp = X; 00452 podarray<int> ipiv(tmp.n_rows); 00453 00454 atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr()); 00455 00456 // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero 00457 eT val = tmp.at(0,0); 00458 for(u32 i=1; i < tmp.n_rows; ++i) 00459 { 00460 val *= tmp.at(i,i); 00461 } 00462 00463 eT sign = eT(1); 00464 for(u32 i=0; i < tmp.n_rows; ++i) 00465 { 00466 if(s32(i) != ipiv.mem[i]) 00467 { 00468 sign *= eT(-1); 00469 } 00470 } 00471 00472 return sign*val; 00473 } 00474 #elif defined(ARMA_USE_LAPACK) 00475 { 00476 Mat<eT> tmp = X; 00477 podarray<int> ipiv(tmp.n_rows); 00478 00479 int info = 0; 00480 int n_rows = tmp.n_rows; 00481 int n_cols = tmp.n_cols; 00482 00483 lapack::getrf_(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info); 00484 00485 // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero 00486 eT val = tmp.at(0,0); 00487 for(u32 i=1; i < tmp.n_rows; ++i) 00488 { 00489 val *= tmp.at(i,i); 00490 } 00491 00492 eT sign = eT(1); 00493 for(u32 i=0; i < tmp.n_rows; ++i) 00494 { 00495 if(s32(i) != ipiv.mem[i]) 00496 { 00497 sign *= eT(-1); 00498 } 00499 } 00500 00501 return sign*val; 00502 } 00503 #else 00504 { 00505 arma_stop("auxlib::det(): need ATLAS or LAPACK library"); 00506 return eT(0); 00507 } 00508 #endif 00509 } 00510 } 00511 }
void auxlib::lu | ( | Mat< eT > & | L, | |
Mat< eT > & | U, | |||
podarray< int > & | ipiv, | |||
const Mat< eT > & | X_orig | |||
) | [inline, static, inherited] |
immediate LU decomposition of a matrix using ATLAS or LAPACK
Definition at line 519 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), lapack::getrf_(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and podarray< T1 >::set_size().
00520 { 00521 arma_extra_debug_sigprint(); 00522 00523 arma_debug_check( (&L == &U), "auxlib::lu(): L and U are the same object"); 00524 00525 unwrap_check< Mat<eT> > tmp1(X_orig, U); 00526 const Mat<eT>& X_tmp = tmp1.M; 00527 00528 unwrap_check< Mat<eT> > tmp2(X_tmp, L); 00529 const Mat<eT>& X = tmp2.M; 00530 00531 U = X; 00532 00533 #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK) 00534 { 00535 00536 #if defined(ARMA_USE_ATLAS) 00537 { 00538 ipiv.set_size(U.n_rows); 00539 00540 //int info = 00541 atlas::clapack_getrf(atlas::CblasColMajor, U.n_rows, U.n_cols, U.memptr(), U.n_rows, ipiv.memptr()); 00542 } 00543 #elif defined(ARMA_USE_LAPACK) 00544 { 00545 ipiv.set_size(U.n_rows); 00546 int info = 0; 00547 00548 int n_rows = U.n_rows; 00549 int n_cols = U.n_cols; 00550 00551 lapack::getrf_(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info); 00552 } 00553 #endif 00554 00555 00556 L.set_size(U.n_rows, U.n_rows); 00557 00558 for(u32 col=0; col<L.n_cols; ++col) 00559 { 00560 00561 for(u32 row=0; row<col; ++row) 00562 { 00563 L.at(row,col) = eT(0); 00564 } 00565 00566 L.at(col,col) = eT(1); 00567 00568 for(u32 row=col+1; row<L.n_rows; ++row) 00569 { 00570 L.at(row,col) = U.at(row,col); 00571 U.at(row,col) = eT(0); 00572 } 00573 00574 } 00575 } 00576 #else 00577 { 00578 arma_stop("auxlib::lu(): need ATLAS or LAPACK library"); 00579 } 00580 #endif 00581 00582 }
void auxlib::lu | ( | Mat< eT > & | L, | |
Mat< eT > & | U, | |||
Mat< eT > & | P, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 589 of file auxlib_meat.hpp.
References eye(), lu(), podarray< T1 >::n_elem, Mat< eT >::set_size(), and Mat< eT >::swap_rows().
00590 { 00591 arma_extra_debug_sigprint(); 00592 00593 arma_debug_check( ( (&P == &L) || (&P == &U) ), "auxlib::lu(): P is the same object as L or U" ); 00594 00595 podarray<int> ipiv; 00596 auxlib::lu(L, U, ipiv, X); 00597 00598 const u32 n = ipiv.n_elem; 00599 00600 P.set_size(n,n); 00601 Mat<eT> ident = eye(n,n); 00602 00603 for(u32 i=n; i>0; --i) 00604 { 00605 const u32 j = i-1; 00606 const u32 k = ipiv[j]-1; 00607 00608 ident.swap_rows(j,k); 00609 00610 if(i == n) 00611 P = ident; 00612 else 00613 P *= ident; 00614 00615 ident.swap_rows(j,k); 00616 } 00617 00618 }
void auxlib::lu | ( | Mat< eT > & | L, | |
Mat< eT > & | U, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 625 of file auxlib_meat.hpp.
References lu().
00626 { 00627 arma_extra_debug_sigprint(); 00628 00629 podarray<int> ipiv; 00630 auxlib::lu(L, U, ipiv, X); 00631 }
void auxlib::eig_sym | ( | Col< eT > & | eigval, | |
const Mat< eT > & | A | |||
) | [inline, static, inherited] |
immediate eigenvalues of a symmetric real matrix using LAPACK
Definition at line 639 of file auxlib_meat.hpp.
References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().
Referenced by eig_sym().
00640 { 00641 arma_extra_debug_sigprint(); 00642 00643 #if defined(ARMA_USE_LAPACK) 00644 { 00645 const unwrap_check<Mat<eT> > tmp(A_orig, eigval); 00646 const Mat<eT>& A = tmp.M; 00647 00648 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square"); 00649 00650 // rudimentary "better-than-nothing" test for symmetry 00651 //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" ); 00652 00653 char jobz = 'N'; 00654 char uplo = 'U'; 00655 00656 int n_rows = A.n_rows; 00657 int lwork = (std::max)(1,3*n_rows-1); 00658 00659 eigval.set_size(n_rows); 00660 podarray<eT> work(lwork); 00661 00662 Mat<eT> A_copy = A; 00663 int info; 00664 00665 arma_extra_debug_print("lapack::syev_()"); 00666 lapack::syev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info); 00667 } 00668 #else 00669 { 00670 arma_stop("auxlib::eig_sym(): need LAPACK library"); 00671 } 00672 #endif 00673 }
void auxlib::eig_sym | ( | Col< T > & | eigval, | |
const Mat< std::complex< T > > & | A | |||
) | [inline, static, inherited] |
immediate eigenvalues of a hermitian complex matrix using LAPACK
Definition at line 681 of file auxlib_meat.hpp.
References arma_stop(), lapack::heev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), and Col< eT >::set_size().
00682 { 00683 arma_extra_debug_sigprint(); 00684 00685 typedef typename std::complex<T> eT; 00686 00687 #if defined(ARMA_USE_LAPACK) 00688 { 00689 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian"); 00690 00691 char jobz = 'N'; 00692 char uplo = 'U'; 00693 00694 int n_rows = A.n_rows; 00695 int lda = A.n_rows; 00696 int lwork = (std::max)(1, 2*n_rows - 1); // TODO: automatically find best size of lwork 00697 00698 eigval.set_size(n_rows); 00699 00700 podarray<eT> work(lwork); 00701 podarray<T> rwork( (std::max)(1, 3*n_rows - 2) ); 00702 00703 Mat<eT> A_copy = A; 00704 int info; 00705 00706 arma_extra_debug_print("lapack::heev_()"); 00707 lapack::heev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info); 00708 } 00709 #else 00710 { 00711 arma_stop("auxlib::eig_sym(): need LAPACK library"); 00712 } 00713 #endif 00714 }
void auxlib::eig_sym | ( | Col< eT > & | eigval, | |
Mat< eT > & | eigvec, | |||
const Mat< eT > & | A | |||
) | [inline, static, inherited] |
immediate eigenvalues and eigenvectors of a symmetric real matrix using LAPACK
Definition at line 722 of file auxlib_meat.hpp.
References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().
00723 { 00724 arma_extra_debug_sigprint(); 00725 00726 // TODO: check for aliasing 00727 00728 #if defined(ARMA_USE_LAPACK) 00729 { 00730 const unwrap_check< Mat<eT> > tmp1(A_orig, eigval); 00731 const Mat<eT>& A_tmp = tmp1.M; 00732 00733 const unwrap_check< Mat<eT> > tmp2(A_tmp, eigvec); 00734 const Mat<eT>& A = tmp2.M; 00735 00736 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square" ); 00737 00738 // rudimentary "better-than-nothing" test for symmetry 00739 //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" ); 00740 00741 00742 char jobz = 'V'; 00743 char uplo = 'U'; 00744 00745 int n_rows = A.n_rows; 00746 int lwork = (std::max)(1, 3*n_rows-1); 00747 00748 eigval.set_size(n_rows); 00749 podarray<eT> work(lwork); 00750 00751 eigvec = A; 00752 int info; 00753 00754 arma_extra_debug_print("lapack::syev_()"); 00755 lapack::syev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info); 00756 } 00757 #else 00758 { 00759 arma_stop("auxlib::eig_sym(): need LAPACK library"); 00760 } 00761 #endif 00762 00763 }
void auxlib::eig_sym | ( | Col< T > & | eigval, | |
Mat< std::complex< T > > & | eigvec, | |||
const Mat< std::complex< T > > & | A | |||
) | [inline, static, inherited] |
immediate eigenvalues and eigenvectors of a hermitian complex matrix using LAPACK
Definition at line 771 of file auxlib_meat.hpp.
References arma_stop(), lapack::heev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Col< eT >::set_size().
00772 { 00773 arma_extra_debug_sigprint(); 00774 00775 typedef typename std::complex<T> eT; 00776 00777 #if defined(ARMA_USE_LAPACK) 00778 { 00779 const unwrap_check< Mat<eT> > tmp(A_orig, eigvec); 00780 const Mat<eT>& A = tmp.M; 00781 00782 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian" ); 00783 00784 char jobz = 'V'; 00785 char uplo = 'U'; 00786 00787 int n_rows = A.n_rows; 00788 int lda = A.n_rows; 00789 int lwork = (std::max)(1, 2*n_rows - 1); // TODO: automatically find best size of lwork 00790 00791 eigval.set_size(n_rows); 00792 00793 podarray<eT> work(lwork); 00794 podarray<T> rwork( (std::max)(1, 3*n_rows - 2) ); 00795 00796 eigvec = A; 00797 int info; 00798 00799 arma_extra_debug_print("lapack::heev_()"); 00800 lapack::heev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info); 00801 } 00802 #else 00803 { 00804 arma_stop("auxlib::eig_sym(): need LAPACK library"); 00805 } 00806 #endif 00807 00808 }
void auxlib::eig_gen | ( | Col< std::complex< T > > & | eigval, | |
Mat< T > & | l_eigvec, | |||
Mat< T > & | r_eigvec, | |||
const Mat< T > & | A, | |||
const char | side | |||
) | [inline, inherited] |
Eigenvalues and eigenvectors of a general square real matrix using LAPACK. The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).
Definition at line 819 of file auxlib_meat.hpp.
References arma_stop(), lapack::geev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().
00826 { 00827 arma_extra_debug_sigprint(); 00828 00829 // TODO: check for aliasing 00830 00831 #if defined(ARMA_USE_LAPACK) 00832 { 00833 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" ); 00834 00835 char jobvl; 00836 char jobvr; 00837 00838 switch(side) 00839 { 00840 case 'l': // left 00841 jobvl = 'V'; 00842 jobvr = 'N'; 00843 break; 00844 00845 case 'r': // right 00846 jobvl = 'N'; 00847 jobvr = 'V'; 00848 break; 00849 00850 case 'b': // both 00851 jobvl = 'V'; 00852 jobvr = 'V'; 00853 break; 00854 00855 case 'n': // neither 00856 jobvl = 'N'; 00857 jobvr = 'N'; 00858 break; 00859 00860 default: 00861 arma_stop("auxlib::eig_gen(): parameter 'side' is invalid"); 00862 } 00863 00864 00865 int n_rows = A.n_rows; 00866 int lda = A.n_rows; 00867 int lwork = (std::max)(1, 4*n_rows); // TODO: automatically find best size of lwork 00868 00869 eigval.set_size(n_rows); 00870 l_eigvec.set_size(n_rows, n_rows); 00871 r_eigvec.set_size(n_rows, n_rows); 00872 00873 podarray<T> work(lwork); 00874 podarray<T> rwork( (std::max)(1, 3*n_rows) ); 00875 00876 podarray<T> wr(n_rows); 00877 podarray<T> wi(n_rows); 00878 00879 Mat<T> A_copy = A; 00880 int info; 00881 00882 arma_extra_debug_print("lapack::cx_geev_()"); 00883 lapack::geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, wr.memptr(), wi.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, &info); 00884 00885 00886 eigval.set_size(n_rows); 00887 for(u32 i=0; i<u32(n_rows); ++i) 00888 { 00889 eigval[i] = std::complex<T>(wr[i], wi[i]); 00890 } 00891 00892 } 00893 #else 00894 { 00895 arma_stop("auxlib::eig_gen(): need LAPACK library"); 00896 } 00897 #endif 00898 00899 }
void auxlib::eig_gen | ( | Col< std::complex< T > > & | eigval, | |
Mat< std::complex< T > > & | l_eigvec, | |||
Mat< std::complex< T > > & | r_eigvec, | |||
const Mat< std::complex< T > > & | A, | |||
const char | side | |||
) | [inline, static, inherited] |
Eigenvalues and eigenvectors of a general square complex matrix using LAPACK The argument 'side' specifies which eigenvectors should be calculated (see code for mode details).
Definition at line 912 of file auxlib_meat.hpp.
References arma_stop(), lapack::cx_geev_(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().
00919 { 00920 arma_extra_debug_sigprint(); 00921 00922 // TODO: check for aliasing 00923 00924 typedef typename std::complex<T> eT; 00925 00926 #if defined(ARMA_USE_LAPACK) 00927 { 00928 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" ); 00929 00930 char jobvl; 00931 char jobvr; 00932 00933 switch(side) 00934 { 00935 case 'l': // left 00936 jobvl = 'V'; 00937 jobvr = 'N'; 00938 break; 00939 00940 case 'r': // right 00941 jobvl = 'N'; 00942 jobvr = 'V'; 00943 break; 00944 00945 case 'b': // both 00946 jobvl = 'V'; 00947 jobvr = 'V'; 00948 break; 00949 00950 case 'n': // neither 00951 jobvl = 'N'; 00952 jobvr = 'N'; 00953 break; 00954 00955 default: 00956 arma_stop("auxlib::eig_gen(): parameter 'side' is invalid"); 00957 } 00958 00959 00960 int n_rows = A.n_rows; 00961 int lda = A.n_rows; 00962 int lwork = (std::max)(1, 4*n_rows); // TODO: automatically find best size of lwork 00963 00964 eigval.set_size(n_rows); 00965 l_eigvec.set_size(n_rows, n_rows); 00966 r_eigvec.set_size(n_rows, n_rows); 00967 00968 podarray<eT> work(lwork); 00969 podarray<T> rwork( (std::max)(1, 3*n_rows) ); // was 2,3 00970 00971 Mat<eT> A_copy = A; 00972 int info; 00973 00974 arma_extra_debug_print("lapack::cx_geev_()"); 00975 lapack::cx_geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, rwork.memptr(), &info); 00976 } 00977 #else 00978 { 00979 arma_stop("auxlib::eig_gen(): need LAPACK library"); 00980 } 00981 #endif 00982 00983 }
bool auxlib::chol | ( | Mat< eT > & | out, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 990 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), Mat< eT >::memptr(), Mat< eT >::n_rows, and lapack::potrf_().
Referenced by chol().
00991 { 00992 arma_extra_debug_sigprint(); 00993 00994 // TODO: check for aliasing 00995 00996 00997 #if defined(ARMA_USE_LAPACK) 00998 { 00999 char uplo = 'U'; 01000 int n = X.n_rows; 01001 int lda = X.n_rows; 01002 int info; 01003 01004 out = X; 01005 lapack::potrf_(&uplo, &n, out.memptr(), &lda, &info); 01006 01007 for(u32 col=0; col<X.n_rows; ++col) 01008 { 01009 eT* colptr = out.colptr(col); 01010 for(u32 row=col+1; row<X.n_rows; ++row) 01011 { 01012 colptr[row] = eT(0); 01013 } 01014 } 01015 01016 return (info == 0); 01017 } 01018 #else 01019 { 01020 arma_stop("auxlib::chol(): need LAPACK library"); 01021 return false; 01022 } 01023 #endif 01024 }
bool auxlib::qr | ( | Mat< eT > & | Q, | |
Mat< eT > & | R, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1031 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), lapack::geqrf_(), max(), Mat< eT >::mem, podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_elem, Mat< eT >::n_rows, lapack::orgqr_(), Mat< eT >::set_size(), podarray< T1 >::set_size(), tmp_real(), and lapack::ungqr_().
Referenced by qr().
01032 { 01033 arma_extra_debug_sigprint(); 01034 01035 // TODO: check for aliasing 01036 01037 #if defined(ARMA_USE_LAPACK) 01038 { 01039 int m = static_cast<int>(X.n_rows); 01040 int n = static_cast<int>(X.n_cols); 01041 int work_len = (std::max)(1,n); 01042 int work_len_tmp; 01043 int k = (std::min)(m,n); 01044 int info; 01045 01046 podarray<eT> tau(k); 01047 podarray<eT> work(work_len); 01048 01049 R = X; 01050 01051 // query for the optimum value of work_len 01052 work_len_tmp = -1; 01053 lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); 01054 01055 if(info == 0) 01056 { 01057 work_len = static_cast<int>(auxlib::tmp_real(work[0])); 01058 work.set_size(work_len); 01059 } 01060 01061 lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); 01062 01063 Q.set_size(X.n_rows,X.n_rows); 01064 01065 eT* Q_mem = Q.memptr(); 01066 const eT* R_mem = R.mem; 01067 01068 const u32 n_elem_copy = (std::min)(Q.n_elem, R.n_elem); 01069 for(u32 i=0; i < n_elem_copy; ++i) 01070 { 01071 Q_mem[i] = R_mem[i]; 01072 } 01073 01074 01075 // construct R 01076 for(u32 row=0; row < R.n_rows; ++row) 01077 { 01078 const u32 n_elem_tmp = (std::min)(row, R.n_cols); 01079 for(u32 col=0; col < n_elem_tmp; ++col) 01080 { 01081 R.at(row,col) = eT(0); 01082 } 01083 } 01084 01085 01086 if( (is_float<eT>::value == true) || (is_double<eT>::value == true) ) 01087 { 01088 // query for the optimum value of work_len 01089 work_len_tmp = -1; 01090 lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); 01091 01092 if(info == 0) 01093 { 01094 work_len = static_cast<int>(auxlib::tmp_real(work[0])); 01095 work.set_size(work_len); 01096 } 01097 01098 lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); 01099 } 01100 else 01101 if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) ) 01102 { 01103 // query for the optimum value of work_len 01104 work_len_tmp = -1; 01105 lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); 01106 01107 if(info == 0) 01108 { 01109 work_len = static_cast<int>(auxlib::tmp_real(work[0])); 01110 work.set_size(work_len); 01111 } 01112 01113 lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); 01114 } 01115 01116 return (info == 0); 01117 } 01118 #else 01119 { 01120 arma_stop("auxlib::qr(): need LAPACK library"); 01121 return false; 01122 } 01123 #endif 01124 }
bool auxlib::svd | ( | Col< eT > & | S, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1131 of file auxlib_meat.hpp.
References arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), and Col< eT >::set_size().
Referenced by svd().
01132 { 01133 arma_extra_debug_sigprint(); 01134 01135 // TODO: check for aliasing 01136 01137 #if defined(ARMA_USE_LAPACK) 01138 { 01139 Mat<eT> A = X; 01140 01141 char jobu = 'N'; 01142 char jobvt = 'N'; 01143 01144 int m = A.n_rows; 01145 int n = A.n_cols; 01146 int lda = A.n_rows; 01147 int ldu = A.n_rows; 01148 int ldvt = A.n_cols; 01149 int lwork = 2; 01150 int info; 01151 01152 Mat<eT> U(1,1); 01153 Mat<eT> V(1,1); 01154 01155 S.set_size( (std::min)(m, n) ); 01156 01157 podarray<eT> work(lwork); 01158 01159 01160 // let gesvd_() calculate the optimum size of the workspace 01161 int lwork_tmp = -1; 01162 01163 lapack::gesvd_<eT> 01164 ( 01165 &jobu, &jobvt, 01166 &m,&n, 01167 A.memptr(), &lda, 01168 S.memptr(), 01169 U.memptr(), &ldu, 01170 V.memptr(), &ldvt, 01171 work.memptr(), &lwork_tmp, 01172 &info 01173 ); 01174 01175 if(info == 0) 01176 { 01177 lwork = static_cast<int>(work[0]); 01178 work.set_size(lwork); 01179 01180 lapack::gesvd_<eT> 01181 ( 01182 &jobu, &jobvt, 01183 &m, &n, 01184 A.memptr(), &lda, 01185 S.memptr(), 01186 U.memptr(), &ldu, 01187 V.memptr(), &ldvt, 01188 work.memptr(), &lwork, 01189 &info 01190 ); 01191 01192 return (info == 0); 01193 } 01194 else 01195 { 01196 return false; 01197 } 01198 } 01199 #else 01200 { 01201 arma_stop("auxlib::svd(): need LAPACK library"); 01202 return false; 01203 } 01204 #endif 01205 }
bool auxlib::svd | ( | Col< T > & | S, | |
const Mat< std::complex< T > > & | X | |||
) | [inline, static, inherited] |
Definition at line 1212 of file auxlib_meat.hpp.
References arma_stop(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< T1 >::set_size(), and Col< eT >::set_size().
01213 { 01214 arma_extra_debug_sigprint(); 01215 01216 typedef std::complex<T> eT; 01217 01218 // TODO: check for aliasing 01219 01220 #if defined(ARMA_USE_LAPACK) 01221 { 01222 Mat<eT> A = X; 01223 01224 char jobu = 'N'; 01225 char jobvt = 'N'; 01226 01227 int m = A.n_rows; 01228 int n = A.n_cols; 01229 int lda = A.n_rows; 01230 int ldu = A.n_rows; 01231 int ldvt = A.n_cols; 01232 int lwork = 2 * (std::min)(m,n) + (std::max)(m,n); 01233 int info; 01234 01235 Mat<eT> U(1,1); 01236 Mat<eT> V(1,1); 01237 01238 S.set_size( (std::min)(m,n) ); 01239 01240 podarray<eT> work(lwork); 01241 podarray<T> rwork( 5*(std::min)(m,n) ); 01242 01243 // let gesvd_() calculate the optimum size of the workspace 01244 int lwork_tmp = -1; 01245 01246 lapack::cx_gesvd_<T> 01247 ( 01248 &jobu, &jobvt, 01249 &m, &n, 01250 A.memptr(), &lda, 01251 S.memptr(), 01252 U.memptr(), &ldu, 01253 V.memptr(), &ldvt, 01254 work.memptr(), &lwork_tmp, 01255 rwork.memptr(), 01256 &info 01257 ); 01258 01259 if(info == 0) 01260 { 01261 int proposed_lwork = static_cast<int>(real(work[0])); 01262 if(proposed_lwork > lwork) 01263 { 01264 lwork = proposed_lwork; 01265 work.set_size(lwork); 01266 } 01267 01268 lapack::cx_gesvd_<T> 01269 ( 01270 &jobu, &jobvt, 01271 &m, &n, 01272 A.memptr(), &lda, 01273 S.memptr(), 01274 U.memptr(), &ldu, 01275 V.memptr(), &ldvt, 01276 work.memptr(), &lwork, 01277 rwork.memptr(), 01278 &info 01279 ); 01280 01281 return (info == 0); 01282 } 01283 else 01284 { 01285 return false; 01286 } 01287 } 01288 #else 01289 { 01290 arma_stop("auxlib::svd(): need LAPACK library"); 01291 return false; 01292 } 01293 #endif 01294 }
bool auxlib::svd | ( | Mat< eT > & | U, | |
Col< eT > & | S, | |||
Mat< eT > & | V, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1301 of file auxlib_meat.hpp.
References op_trans::apply(), arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< T1 >::set_size(), Col< eT >::set_size(), and Mat< eT >::set_size().
01302 { 01303 arma_extra_debug_sigprint(); 01304 01305 // TODO: check for aliasing 01306 01307 #if defined(ARMA_USE_LAPACK) 01308 { 01309 Mat<eT> A = X; 01310 01311 char jobu = 'A'; 01312 char jobvt = 'A'; 01313 01314 int m = A.n_rows; 01315 int n = A.n_cols; 01316 int lda = A.n_rows; 01317 int ldu = A.n_rows; 01318 int ldvt = A.n_cols; 01319 int lwork = 2; 01320 int info; 01321 01322 U.set_size(m,m); 01323 V.set_size(n,n); 01324 01325 S.set_size( (std::min)(m,n) ); 01326 podarray<eT> work(lwork); 01327 01328 // let gesvd_() calculate the optimum size of the workspace 01329 int lwork_tmp = -1; 01330 01331 lapack::gesvd_<eT> 01332 ( 01333 &jobu, &jobvt, 01334 &m, &n, 01335 A.memptr(), &lda, 01336 S.memptr(), 01337 U.memptr(), &ldu, 01338 V.memptr(), &ldvt, 01339 work.memptr(), &lwork_tmp, 01340 &info 01341 ); 01342 01343 if(info == 0) 01344 { 01345 lwork = static_cast<int>(work[0]); 01346 work.set_size(lwork); 01347 01348 lapack::gesvd_<eT> 01349 ( 01350 &jobu, &jobvt, 01351 &m, &n, 01352 A.memptr(), &lda, 01353 S.memptr(), 01354 U.memptr(), &ldu, 01355 V.memptr(), &ldvt, 01356 work.memptr(), &lwork, 01357 &info 01358 ); 01359 01360 op_trans::apply(V,V); // op_trans will work out that an in-place transpose can be done 01361 01362 return (info == 0); 01363 } 01364 else 01365 { 01366 return false; 01367 } 01368 } 01369 #else 01370 { 01371 arma_stop("auxlib::svd(): need LAPACK library"); 01372 return false; 01373 } 01374 #endif 01375 }
bool auxlib::svd | ( | Mat< std::complex< T > > & | U, | |
Col< T > & | S, | |||
Mat< std::complex< T > > & | V, | |||
const Mat< std::complex< T > > & | X | |||
) | [inline, static, inherited] |
Definition at line 1382 of file auxlib_meat.hpp.
References op_htrans::apply(), arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< T1 >::set_size(), and Col< eT >::set_size().
01383 { 01384 arma_extra_debug_sigprint(); 01385 01386 typedef std::complex<T> eT; 01387 01388 // TODO: check for aliasing 01389 01390 #if defined(ARMA_USE_LAPACK) 01391 { 01392 Mat<eT> A = X; 01393 01394 char jobu = 'A'; 01395 char jobvt = 'A'; 01396 01397 int m = A.n_rows; 01398 int n = A.n_cols; 01399 int lda = A.n_rows; 01400 int ldu = A.n_rows; 01401 int ldvt = A.n_cols; 01402 int lwork = 2; 01403 int info; 01404 01405 U.set_size(m,m); 01406 V.set_size(n,n); 01407 01408 S.set_size( (std::min)(m,n) ); 01409 01410 podarray<eT> work(lwork); 01411 podarray<T> rwork( 5*(std::min)(m,n) ); 01412 01413 // let gesvd_() calculate the optimum size of the workspace 01414 int lwork_tmp = -1; 01415 lapack::cx_gesvd_<T> 01416 ( 01417 &jobu, &jobvt, 01418 &m, &n, 01419 A.memptr(), &lda, 01420 S.memptr(), 01421 U.memptr(), &ldu, 01422 V.memptr(), &ldvt, 01423 work.memptr(), &lwork_tmp, 01424 rwork.memptr(), 01425 &info 01426 ); 01427 01428 if(info == 0) 01429 { 01430 lwork = static_cast<int>(real(work[0])); 01431 work.set_size(lwork); 01432 01433 lapack::cx_gesvd_<T> 01434 ( 01435 &jobu, &jobvt, 01436 &m, &n, 01437 A.memptr(), &lda, 01438 S.memptr(), 01439 U.memptr(), &ldu, 01440 V.memptr(), &ldvt, 01441 work.memptr(), &lwork, 01442 rwork.memptr(), 01443 &info 01444 ); 01445 01446 op_htrans::apply(V,V); // op_htrans will work out that an in-place transpose can be done 01447 01448 return (info == 0); 01449 } 01450 else 01451 { 01452 return false; 01453 } 01454 } 01455 #else 01456 { 01457 arma_stop("auxlib::svd(): need LAPACK library"); 01458 return false; 01459 } 01460 #endif 01461 01462 }
bool auxlib::solve | ( | Mat< eT > & | out, | |
const Mat< eT > & | A, | |||
const Mat< eT > & | B | |||
) | [inline, static, inherited] |
Solve a system of linear equations Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows.
Definition at line 1471 of file auxlib_meat.hpp.
References arma_stop(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.
Referenced by solve().
01472 { 01473 arma_extra_debug_sigprint(); 01474 01475 #if defined(ARMA_USE_LAPACK) 01476 { 01477 int n = A.n_rows; 01478 int lda = A.n_rows; 01479 int ldb = A.n_rows; 01480 int nrhs = B.n_cols; 01481 int info; 01482 01483 podarray<int> ipiv(n); 01484 01485 out = B; 01486 Mat<eT> A_copy = A; 01487 01488 lapack::gesv_<eT>(&n, &nrhs, A_copy.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info); 01489 01490 return (info == 0); 01491 } 01492 #else 01493 { 01494 arma_stop("auxlib::solve(): need LAPACK library"); 01495 return false; 01496 } 01497 #endif 01498 }
bool auxlib::solve_od | ( | Mat< eT > & | out, | |
const Mat< eT > & | A, | |||
const Mat< eT > & | B | |||
) | [inline, static, inherited] |
Solve an over-determined system. Assumes that A.n_rows > A.n_cols and B.n_rows = A.n_rows.
Definition at line 1508 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and trans().
Referenced by solve().
01509 { 01510 arma_extra_debug_sigprint(); 01511 01512 #if defined(ARMA_USE_LAPACK) 01513 { 01514 char trans = 'N'; 01515 01516 int m = A.n_rows; 01517 int n = A.n_cols; 01518 int lda = A.n_rows; 01519 int ldb = A.n_rows; 01520 int nrhs = B.n_cols; 01521 int lwork = n + (std::max)(n, nrhs); 01522 int info; 01523 01524 Mat<eT> A_copy = A; 01525 Mat<eT> tmp = B; 01526 01527 01528 podarray<eT> work(lwork); 01529 01530 arma_extra_debug_print("lapack::gels_()"); 01531 01532 // NOTE: 01533 // the dgels() function in the lapack library supplied by ATLAS 3.6 01534 // seems to have problems 01535 01536 lapack::gels_<eT> 01537 ( 01538 &trans, &m, &n, &nrhs, 01539 A_copy.memptr(), &lda, 01540 tmp.memptr(), &ldb, 01541 work.memptr(), &lwork, 01542 &info 01543 ); 01544 01545 arma_extra_debug_print("lapack::gels_() -- finished"); 01546 01547 out.set_size(A.n_cols, B.n_cols); 01548 01549 for(u32 col=0; col<B.n_cols; ++col) 01550 { 01551 syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols ); 01552 } 01553 01554 return (info == 0); 01555 } 01556 #else 01557 { 01558 arma_stop("auxlib::solve_od(): need LAPACK library"); 01559 return false; 01560 } 01561 #endif 01562 }
bool auxlib::solve_ud | ( | Mat< eT > & | out, | |
const Mat< eT > & | A, | |||
const Mat< eT > & | B | |||
) | [inline, static, inherited] |
Solve an under-determined system. Assumes that A.n_rows < A.n_cols and B.n_rows = A.n_rows.
Definition at line 1572 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), max(), podarray< T1 >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), trans(), and Mat< eT >::zeros().
Referenced by solve().
01573 { 01574 arma_extra_debug_sigprint(); 01575 01576 #if defined(ARMA_USE_LAPACK) 01577 { 01578 char trans = 'N'; 01579 01580 int m = A.n_rows; 01581 int n = A.n_cols; 01582 int lda = A.n_rows; 01583 int ldb = A.n_cols; 01584 int nrhs = B.n_cols; 01585 int lwork = m + (std::max)(m,nrhs); 01586 int info; 01587 01588 01589 Mat<eT> A_copy = A; 01590 01591 Mat<eT> tmp; 01592 tmp.zeros(A.n_cols, B.n_cols); 01593 01594 for(u32 col=0; col<B.n_cols; ++col) 01595 { 01596 eT* tmp_colmem = tmp.colptr(col); 01597 01598 syslib::copy_elem( tmp_colmem, B.colptr(col), B.n_rows ); 01599 01600 for(u32 row=B.n_rows; row<A.n_cols; ++row) 01601 { 01602 tmp_colmem[row] = eT(0); 01603 } 01604 } 01605 01606 podarray<eT> work(lwork); 01607 01608 arma_extra_debug_print("lapack::gels_()"); 01609 01610 // NOTE: 01611 // the dgels() function in the lapack library supplied by ATLAS 3.6 01612 // seems to have problems 01613 01614 lapack::gels_<eT> 01615 ( 01616 &trans, &m, &n, &nrhs, 01617 A_copy.memptr(), &lda, 01618 tmp.memptr(), &ldb, 01619 work.memptr(), &lwork, 01620 &info 01621 ); 01622 01623 arma_extra_debug_print("lapack::gels_() -- finished"); 01624 01625 out.set_size(A.n_cols, B.n_cols); 01626 01627 for(u32 col=0; col<B.n_cols; ++col) 01628 { 01629 syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols ); 01630 } 01631 01632 return (info == 0); 01633 } 01634 #else 01635 { 01636 arma_stop("auxlib::solve_ud(): need LAPACK library"); 01637 return false; 01638 } 01639 #endif 01640 }