//! More...
Classes | |
class | auxlib |
wrapper for accessing external functions defined in ATLAS, LAPACK or BLAS libraries More... | |
Functions | |
template<typename eT > | |
static bool | auxlib::inv_noalias (Mat< eT > &out, const Mat< eT > &X) |
immediate matrix inverse | |
template<typename eT > | |
static bool | 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::log_det (eT &out_val, typename get_pod_type< eT >::result &out_sign, const Mat< eT > &X) |
immediate log 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. |
//!
bool auxlib::inv_noalias | ( | Mat< eT > & | out, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
immediate matrix inverse
Definition at line 26 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< eT >::set_size(), Mat< eT >::set_size(), and atlas::tmp_real().
Referenced by op_inv::apply().
{ arma_extra_debug_sigprint(); switch(X.n_rows) { case 1: { out.set_size(1,1); out[0] = eT(1) / X[0]; }; break; case 2: { out.set_size(2,2); const eT a = X.at(0,0); const eT b = X.at(0,1); const eT c = X.at(1,0); const eT d = X.at(1,1); const eT k = eT(1) / (a*d - b*c); out.at(0,0) = d*k; out.at(0,1) = -b*k; out.at(1,0) = -c*k; out.at(1,1) = a*k; }; break; case 3: { out.set_size(3,3); const eT* X_col0 = X.colptr(0); const eT a11 = X_col0[0]; const eT a21 = X_col0[1]; const eT a31 = X_col0[2]; const eT* X_col1 = X.colptr(1); const eT a12 = X_col1[0]; const eT a22 = X_col1[1]; const eT a32 = X_col1[2]; const eT* X_col2 = X.colptr(2); const eT a13 = X_col2[0]; const eT a23 = X_col2[1]; const eT a33 = X_col2[2]; const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); eT* out_col0 = out.colptr(0); out_col0[0] = (a33*a22 - a32*a23) * k; out_col0[1] = -(a33*a21 - a31*a23) * k; out_col0[2] = (a32*a21 - a31*a22) * k; eT* out_col1 = out.colptr(1); out_col1[0] = -(a33*a12 - a32*a13) * k; out_col1[1] = (a33*a11 - a31*a13) * k; out_col1[2] = -(a32*a11 - a31*a12) * k; eT* out_col2 = out.colptr(2); out_col2[0] = (a23*a12 - a22*a13) * k; out_col2[1] = -(a23*a11 - a21*a13) * k; out_col2[2] = (a22*a11 - a21*a12) * k; }; break; case 4: { out.set_size(4,4); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); out /= det(X); }; break; default: { #if defined(ARMA_USE_ATLAS) { out = X; podarray<int> ipiv(out.n_rows); int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr()); if(info == 0) { info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr()); } return (info == 0); } #elif defined(ARMA_USE_LAPACK) { out = X; int n_rows = out.n_rows; int n_cols = out.n_cols; int info = 0; podarray<int> ipiv(out.n_rows); // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6) // based on tests with various matrix types on 32-bit and 64-bit machines // // the "work" array is deliberately long so that a secondary (time-consuming) // memory allocation is avoided, if possible int work_len = (std::max)(1, n_rows*84); podarray<eT> work(work_len); lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info); if(info == 0) { // query for optimum size of work_len int work_len_tmp = -1; lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info); if(info == 0) { int proposed_work_len = static_cast<int>(access::tmp_real(work[0])); // if necessary, allocate more memory if(work_len < proposed_work_len) { work_len = proposed_work_len; work.set_size(work_len); } } lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info); } return (info == 0); } #else { arma_stop("inv(): need ATLAS or LAPACK"); } #endif }; } return true; }
bool auxlib::inv_inplace | ( | Mat< eT > & | X | ) | [inline, static, inherited] |
immediate inplace matrix inverse
Definition at line 202 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), atlas::clapack_getri(), Mat< eT >::colptr(), det(), lapack::getrf_(), lapack::getri_(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< eT >::set_size(), and atlas::tmp_real().
Referenced by op_inv::apply().
{ arma_extra_debug_sigprint(); // for more info, see: // http://www.dr-lex.34sp.com/random/matrix_inv.html // http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html // http://www.euclideanspace.com/maths/algebra/matrix/functions/inverse/fourD/index.htm // http://www.geometrictools.com//LibFoundation/Mathematics/Wm4Matrix4.inl switch(X.n_rows) { case 1: { X[0] = eT(1) / X[0]; }; break; case 2: { const eT a = X.at(0,0); const eT b = X.at(0,1); const eT c = X.at(1,0); const eT d = X.at(1,1); const eT k = eT(1) / (a*d - b*c); X.at(0,0) = d*k; X.at(0,1) = -b*k; X.at(1,0) = -c*k; X.at(1,1) = a*k; }; break; case 3: { eT* X_col0 = X.colptr(0); eT* X_col1 = X.colptr(1); eT* X_col2 = X.colptr(2); const eT a11 = X_col0[0]; const eT a21 = X_col0[1]; const eT a31 = X_col0[2]; const eT a12 = X_col1[0]; const eT a22 = X_col1[1]; const eT a32 = X_col1[2]; const eT a13 = X_col2[0]; const eT a23 = X_col2[1]; const eT a33 = X_col2[2]; const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); X_col0[0] = (a33*a22 - a32*a23) * k; X_col0[1] = -(a33*a21 - a31*a23) * k; X_col0[2] = (a32*a21 - a31*a22) * k; X_col1[0] = -(a33*a12 - a32*a13) * k; X_col1[1] = (a33*a11 - a31*a13) * k; X_col1[2] = -(a32*a11 - a31*a12) * k; X_col2[0] = (a23*a12 - a22*a13) * k; X_col2[1] = -(a23*a11 - a21*a13) * k; X_col2[2] = (a22*a11 - a21*a12) * k; }; break; case 4: { const Mat<eT> A(X); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); 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); X /= det(A); }; break; default: { #if defined(ARMA_USE_ATLAS) { Mat<eT>& out = X; podarray<int> ipiv(out.n_rows); int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr()); if(info == 0) { info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr()); } return (info == 0); } #elif defined(ARMA_USE_LAPACK) { Mat<eT>& out = X; int n_rows = out.n_rows; int n_cols = out.n_cols; int info = 0; podarray<int> ipiv(out.n_rows); // 84 was empirically found -- it is the maximum value suggested by LAPACK (as provided by ATLAS v3.6) // based on tests with various matrix types on 32-bit and 64-bit machines // // the "work" array is deliberately long so that a secondary (time-consuming) // memory allocation is avoided, if possible int work_len = (std::max)(1, n_rows*84); podarray<eT> work(work_len); lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info); if(info == 0) { // query for optimum size of work_len int work_len_tmp = -1; lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info); if(info == 0) { int proposed_work_len = static_cast<int>(access::tmp_real(work[0])); // if necessary, allocate more memory if(work_len < proposed_work_len) { work_len = proposed_work_len; work.set_size(work_len); } } lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info); } return (info == 0); } #else { arma_stop("inv(): need ATLAS or LAPACK"); } #endif } } return true; }
eT auxlib::det | ( | const Mat< eT > & | X | ) | [inline, static, inherited] |
immediate determinant of a matrix using ATLAS or LAPACK
Definition at line 376 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), Mat< eT >::colptr(), lapack::getrf_(), podarray< eT >::mem, podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.
Referenced by det(), inv_inplace(), and inv_noalias().
{ arma_extra_debug_sigprint(); switch(X.n_rows) { case 0: return eT(0); case 1: return X[0]; case 2: return (X.at(0,0)*X.at(1,1) - X.at(0,1)*X.at(1,0)); case 3: { const eT* a_col0 = X.colptr(0); const eT a11 = a_col0[0]; const eT a21 = a_col0[1]; const eT a31 = a_col0[2]; const eT* a_col1 = X.colptr(1); const eT a12 = a_col1[0]; const eT a22 = a_col1[1]; const eT a32 = a_col1[2]; const eT* a_col2 = X.colptr(2); const eT a13 = a_col2[0]; const eT a23 = a_col2[1]; const eT a33 = a_col2[2]; return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) ); // const double tmp1 = X.at(0,0) * X.at(1,1) * X.at(2,2); // const double tmp2 = X.at(0,1) * X.at(1,2) * X.at(2,0); // const double tmp3 = X.at(0,2) * X.at(1,0) * X.at(2,1); // const double tmp4 = X.at(2,0) * X.at(1,1) * X.at(0,2); // const double tmp5 = X.at(2,1) * X.at(1,2) * X.at(0,0); // const double tmp6 = X.at(2,2) * X.at(1,0) * X.at(0,1); // return (tmp1+tmp2+tmp3) - (tmp4+tmp5+tmp6); } case 4: { const eT val = \ X.at(0,3) * X.at(1,2) * X.at(2,1) * X.at(3,0) \ - X.at(0,2) * X.at(1,3) * X.at(2,1) * X.at(3,0) \ - X.at(0,3) * X.at(1,1) * X.at(2,2) * X.at(3,0) \ + X.at(0,1) * X.at(1,3) * X.at(2,2) * X.at(3,0) \ + X.at(0,2) * X.at(1,1) * X.at(2,3) * X.at(3,0) \ - X.at(0,1) * X.at(1,2) * X.at(2,3) * X.at(3,0) \ - X.at(0,3) * X.at(1,2) * X.at(2,0) * X.at(3,1) \ + X.at(0,2) * X.at(1,3) * X.at(2,0) * X.at(3,1) \ + X.at(0,3) * X.at(1,0) * X.at(2,2) * X.at(3,1) \ - X.at(0,0) * X.at(1,3) * X.at(2,2) * X.at(3,1) \ - X.at(0,2) * X.at(1,0) * X.at(2,3) * X.at(3,1) \ + X.at(0,0) * X.at(1,2) * X.at(2,3) * X.at(3,1) \ + X.at(0,3) * X.at(1,1) * X.at(2,0) * X.at(3,2) \ - X.at(0,1) * X.at(1,3) * X.at(2,0) * X.at(3,2) \ - X.at(0,3) * X.at(1,0) * X.at(2,1) * X.at(3,2) \ + X.at(0,0) * X.at(1,3) * X.at(2,1) * X.at(3,2) \ + X.at(0,1) * X.at(1,0) * X.at(2,3) * X.at(3,2) \ - X.at(0,0) * X.at(1,1) * X.at(2,3) * X.at(3,2) \ - X.at(0,2) * X.at(1,1) * X.at(2,0) * X.at(3,3) \ + X.at(0,1) * X.at(1,2) * X.at(2,0) * X.at(3,3) \ + X.at(0,2) * X.at(1,0) * X.at(2,1) * X.at(3,3) \ - X.at(0,0) * X.at(1,2) * X.at(2,1) * X.at(3,3) \ - X.at(0,1) * X.at(1,0) * X.at(2,2) * X.at(3,3) \ + X.at(0,0) * X.at(1,1) * X.at(2,2) * X.at(3,3) \ ; return val; } default: { #if defined(ARMA_USE_ATLAS) { Mat<eT> tmp = X; podarray<int> ipiv(tmp.n_rows); atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr()); // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero eT val = tmp.at(0,0); for(u32 i=1; i < tmp.n_rows; ++i) { val *= tmp.at(i,i); } int sign = +1; for(u32 i=0; i < tmp.n_rows; ++i) { if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0 { sign *= -1; } } return ( (sign < 0) ? -val : val ); } #elif defined(ARMA_USE_LAPACK) { Mat<eT> tmp = X; podarray<int> ipiv(tmp.n_rows); int info = 0; int n_rows = int(tmp.n_rows); int n_cols = int(tmp.n_cols); lapack::getrf_(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info); // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero eT val = tmp.at(0,0); for(u32 i=1; i < tmp.n_rows; ++i) { val *= tmp.at(i,i); } int sign = +1; for(u32 i=0; i < tmp.n_rows; ++i) { if( int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is required as Fortran counts from 1 { sign *= -1; } } return ( (sign < 0) ? -val : val ); } #else { arma_stop("det(): need ATLAS or LAPACK"); return eT(0); } #endif } } }
void auxlib::log_det | ( | eT & | out_val, | |
typename get_pod_type< eT >::result & | out_sign, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
immediate log determinant of a matrix using ATLAS or LAPACK
Definition at line 523 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), lapack::getrf_(), log(), podarray< eT >::mem, podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and access::tmp_real().
Referenced by log_det().
{ arma_extra_debug_sigprint(); typedef typename get_pod_type<eT>::result T; #if defined(ARMA_USE_ATLAS) { Mat<eT> tmp = X; podarray<int> ipiv(tmp.n_rows); atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr()); // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero s32 sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1; eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) ); for(u32 i=1; i < tmp.n_rows; ++i) { const eT x = tmp.at(i,i); sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1; val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x); } for(u32 i=0; i < tmp.n_rows; ++i) { if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the clapack version of getrf() assumes counting from 0 { sign *= -1; } } out_val = val; out_sign = T(sign); } #elif defined(ARMA_USE_LAPACK) { Mat<eT> tmp = X; podarray<int> ipiv(tmp.n_rows); int info = 0; int n_rows = int(tmp.n_rows); int n_cols = int(tmp.n_cols); lapack::getrf_(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info); // on output tmp appears to be L+U_alt, where U_alt is U with the main diagonal set to zero s32 sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? -1 : +1 ) : +1; eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_real( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( tmp.at(0,0) ); for(u32 i=1; i < tmp.n_rows; ++i) { const eT x = tmp.at(i,i); sign *= (is_complex<eT>::value == false) ? ( (access::tmp_real(x) < T(0)) ? -1 : +1 ) : +1; val += (is_complex<eT>::value == false) ? std::log( (access::tmp_real(x) < T(0)) ? x*T(-1) : x ) : std::log(x); } for(u32 i=0; i < tmp.n_rows; ++i) { if( int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is required as Fortran counts from 1 { sign *= -1; } } out_val = val; out_sign = T(sign); } #else { arma_stop("log_det(): need ATLAS or LAPACK"); out_val = eT(0); out_sign = T(0); } #endif }
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 611 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), atlas::clapack_getrf(), lapack::getrf_(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and podarray< eT >::set_size().
{ arma_extra_debug_sigprint(); U = X; #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK) { #if defined(ARMA_USE_ATLAS) { ipiv.set_size(U.n_rows); //int info = atlas::clapack_getrf(atlas::CblasColMajor, U.n_rows, U.n_cols, U.memptr(), U.n_rows, ipiv.memptr()); } #elif defined(ARMA_USE_LAPACK) { ipiv.set_size(U.n_rows); int info = 0; int n_rows = U.n_rows; int n_cols = U.n_cols; lapack::getrf_(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info); // take into account that Fortran counts from 1 for(u32 i=0; i<U.n_rows; ++i) { ipiv[i] -= 1; } } #endif L.set_size(U.n_rows, U.n_rows); for(u32 col=0; col<L.n_cols; ++col) { for(u32 row=0; row<col; ++row) { L.at(row,col) = eT(0); } L.at(col,col) = eT(1); for(u32 row=col+1; row<L.n_rows; ++row) { L.at(row,col) = U.at(row,col); U.at(row,col) = eT(0); } } } #else { arma_stop("lu(): need ATLAS or LAPACK"); } #endif }
void auxlib::lu | ( | Mat< eT > & | L, | |
Mat< eT > & | U, | |||
Mat< eT > & | P, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 680 of file auxlib_meat.hpp.
References lu(), podarray< eT >::n_elem, and Mat< eT >::swap_rows().
{ arma_extra_debug_sigprint(); podarray<int> ipiv; auxlib::lu(L, U, ipiv, X); const u32 n = ipiv.n_elem; Mat<u32> P_tmp(n,n); Mat<u32> ident = eye< Mat<u32> >(n,n); for(u32 i=n; i>0; --i) { const u32 j = i-1; const u32 k = ipiv[j]; ident.swap_rows(j,k); if(i == n) { P_tmp = ident; } else { P_tmp *= ident; } ident.swap_rows(j,k); } P = conv_to< Mat<eT> >::from(P_tmp); }
void auxlib::lu | ( | Mat< eT > & | L, | |
Mat< eT > & | U, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 719 of file auxlib_meat.hpp.
References lu().
{ arma_extra_debug_sigprint(); podarray<int> ipiv; auxlib::lu(L, U, ipiv, X); }
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 733 of file auxlib_meat.hpp.
References arma_stop(), unwrap_check< T1 >::M, podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().
Referenced by eig_sym().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { const unwrap_check<Mat<eT> > tmp(A_orig, eigval); const Mat<eT>& A = tmp.M; arma_debug_check( (A.n_rows != A.n_cols), "eig_sym(): given matrix is not square"); // rudimentary "better-than-nothing" test for symmetry //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" ); char jobz = 'N'; char uplo = 'U'; int n_rows = A.n_rows; int lwork = (std::max)(1,3*n_rows-1); eigval.set_size(n_rows); podarray<eT> work(lwork); Mat<eT> A_copy = A; int info; arma_extra_debug_print("lapack::syev_()"); lapack::syev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info); } #else { arma_stop("eig_sym(): need LAPACK"); } #endif }
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 775 of file auxlib_meat.hpp.
References arma_stop(), lapack::heev_(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), and Col< eT >::set_size().
{ arma_extra_debug_sigprint(); typedef typename std::complex<T> eT; #if defined(ARMA_USE_LAPACK) { arma_debug_check( (A.n_rows != A.n_cols), "eig_sym(): given matrix is not hermitian"); char jobz = 'N'; char uplo = 'U'; int n_rows = A.n_rows; int lda = A.n_rows; int lwork = (std::max)(1, 2*n_rows - 1); // TODO: automatically find best size of lwork eigval.set_size(n_rows); podarray<eT> work(lwork); podarray<T> rwork( (std::max)(1, 3*n_rows - 2) ); Mat<eT> A_copy = A; int info; arma_extra_debug_print("lapack::heev_()"); lapack::heev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info); } #else { arma_stop("eig_sym(): need LAPACK"); } #endif }
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 816 of file auxlib_meat.hpp.
References arma_stop(), unwrap_check< T1 >::M, podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Col< eT >::set_size(), and lapack::syev_().
{ arma_extra_debug_sigprint(); // TODO: check for aliasing #if defined(ARMA_USE_LAPACK) { const unwrap_check< Mat<eT> > tmp1(A_orig, eigval); const Mat<eT>& A_tmp = tmp1.M; const unwrap_check< Mat<eT> > tmp2(A_tmp, eigvec); const Mat<eT>& A = tmp2.M; arma_debug_check( (A.n_rows != A.n_cols), "eig_sym(): given matrix is not square" ); // rudimentary "better-than-nothing" test for symmetry //arma_debug_check( (A.at(A.n_rows-1, 0) != A.at(0, A.n_cols-1)), "auxlib::eig(): given matrix is not symmetric" ); char jobz = 'V'; char uplo = 'U'; int n_rows = A.n_rows; int lwork = (std::max)(1, 3*n_rows-1); eigval.set_size(n_rows); podarray<eT> work(lwork); eigvec = A; int info; arma_extra_debug_print("lapack::syev_()"); lapack::syev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info); } #else { arma_stop("eig_sym(): need LAPACK"); } #endif }
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 865 of file auxlib_meat.hpp.
References arma_stop(), lapack::heev_(), unwrap_check< T1 >::M, max(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Col< eT >::set_size().
{ arma_extra_debug_sigprint(); typedef typename std::complex<T> eT; #if defined(ARMA_USE_LAPACK) { const unwrap_check< Mat<eT> > tmp(A_orig, eigvec); const Mat<eT>& A = tmp.M; arma_debug_check( (A.n_rows != A.n_cols), "eig_sym(): given matrix is not hermitian" ); char jobz = 'V'; char uplo = 'U'; int n_rows = A.n_rows; int lda = A.n_rows; int lwork = (std::max)(1, 2*n_rows - 1); // TODO: automatically find best size of lwork eigval.set_size(n_rows); podarray<eT> work(lwork); podarray<T> rwork( (std::max)(1, 3*n_rows - 2) ); eigvec = A; int info; arma_extra_debug_print("lapack::heev_()"); lapack::heev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info); } #else { arma_stop("eig_sym(): need LAPACK"); } #endif }
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 913 of file auxlib_meat.hpp.
References arma_stop(), lapack::geev_(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().
{ arma_extra_debug_sigprint(); // TODO: check for aliasing #if defined(ARMA_USE_LAPACK) { arma_debug_check( (A.n_rows != A.n_cols), "eig_gen(): given matrix is not square" ); char jobvl; char jobvr; switch(side) { case 'l': // left jobvl = 'V'; jobvr = 'N'; break; case 'r': // right jobvl = 'N'; jobvr = 'V'; break; case 'b': // both jobvl = 'V'; jobvr = 'V'; break; case 'n': // neither jobvl = 'N'; jobvr = 'N'; break; default: arma_stop("eig_gen(): parameter 'side' is invalid"); } int n_rows = A.n_rows; int lda = A.n_rows; int lwork = (std::max)(1, 4*n_rows); // TODO: automatically find best size of lwork eigval.set_size(n_rows); l_eigvec.set_size(n_rows, n_rows); r_eigvec.set_size(n_rows, n_rows); podarray<T> work(lwork); podarray<T> rwork( (std::max)(1, 3*n_rows) ); podarray<T> wr(n_rows); podarray<T> wi(n_rows); Mat<T> A_copy = A; int info; arma_extra_debug_print("lapack::cx_geev_()"); 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); eigval.set_size(n_rows); for(u32 i=0; i<u32(n_rows); ++i) { eigval[i] = std::complex<T>(wr[i], wi[i]); } } #else { arma_stop("eig_gen(): need LAPACK"); } #endif }
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 1006 of file auxlib_meat.hpp.
References arma_stop(), lapack::cx_geev_(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, and Mat< eT >::set_size().
{ arma_extra_debug_sigprint(); // TODO: check for aliasing typedef typename std::complex<T> eT; #if defined(ARMA_USE_LAPACK) { arma_debug_check( (A.n_rows != A.n_cols), "eig_gen(): given matrix is not square" ); char jobvl; char jobvr; switch(side) { case 'l': // left jobvl = 'V'; jobvr = 'N'; break; case 'r': // right jobvl = 'N'; jobvr = 'V'; break; case 'b': // both jobvl = 'V'; jobvr = 'V'; break; case 'n': // neither jobvl = 'N'; jobvr = 'N'; break; default: arma_stop("eig_gen(): parameter 'side' is invalid"); } int n_rows = A.n_rows; int lda = A.n_rows; int lwork = (std::max)(1, 4*n_rows); // TODO: automatically find best size of lwork eigval.set_size(n_rows); l_eigvec.set_size(n_rows, n_rows); r_eigvec.set_size(n_rows, n_rows); podarray<eT> work(lwork); podarray<T> rwork( (std::max)(1, 3*n_rows) ); // was 2,3 Mat<eT> A_copy = A; int info; arma_extra_debug_print("lapack::cx_geev_()"); 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); } #else { arma_stop("eig_gen(): need LAPACK"); } #endif }
bool auxlib::chol | ( | Mat< eT > & | out, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1084 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), Mat< eT >::memptr(), Mat< eT >::n_rows, and lapack::potrf_().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { char uplo = 'U'; int n = X.n_rows; int lda = X.n_rows; int info; out = X; lapack::potrf_(&uplo, &n, out.memptr(), &lda, &info); for(u32 col=0; col<X.n_rows; ++col) { eT* colptr = out.colptr(col); for(u32 row=col+1; row<X.n_rows; ++row) { colptr[row] = eT(0); } } return (info == 0); } #else { arma_stop("chol(): need LAPACK"); return false; } #endif }
bool auxlib::qr | ( | Mat< eT > & | Q, | |
Mat< eT > & | R, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1122 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::at(), lapack::geqrf_(), Mat< eT >::mem, podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_elem, Mat< eT >::n_rows, lapack::orgqr_(), Mat< eT >::set_size(), podarray< eT >::set_size(), atlas::tmp_real(), and lapack::ungqr_().
Referenced by qr().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { int m = static_cast<int>(X.n_rows); int n = static_cast<int>(X.n_cols); int work_len = (std::max)(1,n); int work_len_tmp; int k = (std::min)(m,n); int info; podarray<eT> tau(k); podarray<eT> work(work_len); R = X; // query for the optimum value of work_len work_len_tmp = -1; lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); if(info == 0) { work_len = static_cast<int>(access::tmp_real(work[0])); work.set_size(work_len); } lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); Q.set_size(X.n_rows, X.n_rows); eT* Q_mem = Q.memptr(); const eT* R_mem = R.mem; const u32 n_elem_copy = (std::min)(Q.n_elem, R.n_elem); for(u32 i=0; i < n_elem_copy; ++i) { Q_mem[i] = R_mem[i]; } // construct R for(u32 row=0; row < R.n_rows; ++row) { const u32 n_elem_tmp = (std::min)(row, R.n_cols); for(u32 col=0; col < n_elem_tmp; ++col) { R.at(row,col) = eT(0); } } if( (is_float<eT>::value == true) || (is_double<eT>::value == true) ) { // query for the optimum value of work_len work_len_tmp = -1; lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); if(info == 0) { work_len = static_cast<int>(access::tmp_real(work[0])); work.set_size(work_len); } lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); } else if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) ) { // query for the optimum value of work_len work_len_tmp = -1; lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info); if(info == 0) { work_len = static_cast<int>(access::tmp_real(work[0])); work.set_size(work_len); } lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info); } return (info == 0); } #else { arma_stop("qr(): need LAPACK"); return false; } #endif }
bool auxlib::svd | ( | Col< eT > & | S, | |
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1220 of file auxlib_meat.hpp.
References arma_stop(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< eT >::set_size(), and Col< eT >::set_size().
Referenced by rank(), and svd().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { Mat<eT> A = X; Mat<eT> U(1, 1); Mat<eT> V(1, A.n_cols); char jobu = 'N'; char jobvt = 'N'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldu = U.n_rows; int ldvt = V.n_rows; int lwork = 2 * (std::max)(1, (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) ); int info; S.set_size( (std::min)(m, n) ); podarray<eT> work(lwork); // let gesvd_() calculate the optimum size of the workspace int lwork_tmp = -1; lapack::gesvd_<eT> ( &jobu, &jobvt, &m,&n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork_tmp, &info ); if(info == 0) { int proposed_lwork = static_cast<int>(work[0]); if(proposed_lwork > lwork) { lwork = proposed_lwork; work.set_size(lwork); } lapack::gesvd_<eT> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, &info ); } return (info == 0); } #else { arma_stop("svd(): need LAPACK"); return false; } #endif }
bool auxlib::svd | ( | Col< T > & | S, | |
const Mat< std::complex< T > > & | X | |||
) | [inline, static, inherited] |
Definition at line 1300 of file auxlib_meat.hpp.
References arma_stop(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< eT >::set_size(), and Col< eT >::set_size().
{ arma_extra_debug_sigprint(); typedef std::complex<T> eT; #if defined(ARMA_USE_LAPACK) { Mat<eT> A = X; Mat<eT> U(1, 1); Mat<eT> V(1, A.n_cols); char jobu = 'N'; char jobvt = 'N'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldu = U.n_rows; int ldvt = V.n_rows; int lwork = 2 * (std::max)(1, 2*(std::min)(m,n)+(std::max)(m,n) ); int info; S.set_size( (std::min)(m,n) ); podarray<eT> work(lwork); podarray<T> rwork( 5*(std::min)(m,n) ); // let gesvd_() calculate the optimum size of the workspace int lwork_tmp = -1; lapack::cx_gesvd_<T> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork_tmp, rwork.memptr(), &info ); if(info == 0) { int proposed_lwork = static_cast<int>(real(work[0])); if(proposed_lwork > lwork) { lwork = proposed_lwork; work.set_size(lwork); } lapack::cx_gesvd_<T> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, rwork.memptr(), &info ); } return (info == 0); } #else { arma_stop("svd(): need LAPACK"); return false; } #endif }
bool auxlib::svd | ( | Mat< eT > & | U, | |
Col< eT > & | S, | |||
Mat< eT > & | V, | |||
const Mat< eT > & | X | |||
) | [inline, static, inherited] |
Definition at line 1383 of file auxlib_meat.hpp.
References op_trans::apply(), arma_stop(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, podarray< eT >::set_size(), Col< eT >::set_size(), and Mat< eT >::set_size().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { Mat<eT> A = X; U.set_size(A.n_rows, A.n_rows); V.set_size(A.n_cols, A.n_cols); char jobu = 'A'; char jobvt = 'A'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldu = U.n_rows; int ldvt = V.n_rows; int lwork = 2 * (std::max)(1, (std::max)( (3*(std::min)(m,n) + (std::max)(m,n)), 5*(std::min)(m,n) ) ); int info; S.set_size( (std::min)(m,n) ); podarray<eT> work(lwork); // let gesvd_() calculate the optimum size of the workspace int lwork_tmp = -1; lapack::gesvd_<eT> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork_tmp, &info ); if(info == 0) { int proposed_lwork = static_cast<int>(work[0]); if(proposed_lwork > lwork) { lwork = proposed_lwork; work.set_size(lwork); } lapack::gesvd_<eT> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, &info ); op_trans::apply(V,V); // op_trans will work out that an in-place transpose can be done } return (info == 0); } #else { arma_stop("svd(): need LAPACK"); return false; } #endif }
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 1463 of file auxlib_meat.hpp.
References op_htrans::apply(), arma_stop(), max(), podarray< eT >::memptr(), Mat< eT >::memptr(), min(), Mat< eT >::n_cols, Mat< eT >::n_rows, real(), podarray< eT >::set_size(), Col< eT >::set_size(), and Mat< eT >::set_size().
{ arma_extra_debug_sigprint(); typedef std::complex<T> eT; #if defined(ARMA_USE_LAPACK) { Mat<eT> A = X; U.set_size(A.n_rows, A.n_rows); V.set_size(A.n_cols, A.n_cols); char jobu = 'A'; char jobvt = 'A'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldu = U.n_rows; int ldvt = V.n_rows; int lwork = 2 * (std::max)(1, 2*(std::min)(m,n)+(std::max)(m,n) ); int info; S.set_size( (std::min)(m,n) ); podarray<eT> work(lwork); podarray<T> rwork( 5*(std::min)(m,n) ); // let gesvd_() calculate the optimum size of the workspace int lwork_tmp = -1; lapack::cx_gesvd_<T> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork_tmp, rwork.memptr(), &info ); if(info == 0) { int proposed_lwork = static_cast<int>(real(work[0])); if(proposed_lwork > lwork) { lwork = proposed_lwork; work.set_size(lwork); } lapack::cx_gesvd_<T> ( &jobu, &jobvt, &m, &n, A.memptr(), &lda, S.memptr(), U.memptr(), &ldu, V.memptr(), &ldvt, work.memptr(), &lwork, rwork.memptr(), &info ); op_htrans::apply(V,V); // op_htrans will work out that an in-place transpose can be done } return (info == 0); } #else { arma_stop("svd(): need LAPACK"); return false; } #endif }
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 1550 of file auxlib_meat.hpp.
References arma_stop(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, and Mat< eT >::n_rows.
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { int n = A.n_rows; int lda = A.n_rows; int ldb = A.n_rows; int nrhs = B.n_cols; int info; podarray<int> ipiv(n); out = B; Mat<eT> A_copy = A; lapack::gesv_<eT>(&n, &nrhs, A_copy.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info); return (info == 0); } #else { arma_stop("solve(): need LAPACK"); return false; } #endif }
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 1587 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), and trans().
Referenced by glue_solve::apply().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { char trans = 'N'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldb = A.n_rows; int nrhs = B.n_cols; int lwork = n + (std::max)(n, nrhs); int info; Mat<eT> A_copy = A; Mat<eT> tmp = B; podarray<eT> work(lwork); arma_extra_debug_print("lapack::gels_()"); // NOTE: // the dgels() function in the lapack library supplied by ATLAS 3.6 // seems to have problems lapack::gels_<eT> ( &trans, &m, &n, &nrhs, A_copy.memptr(), &lda, tmp.memptr(), &ldb, work.memptr(), &lwork, &info ); arma_extra_debug_print("lapack::gels_() -- finished"); out.set_size(A.n_cols, B.n_cols); for(u32 col=0; col<B.n_cols; ++col) { syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols ); } return (info == 0); } #else { arma_stop("auxlib::solve_od(): need LAPACK"); return false; } #endif }
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 1651 of file auxlib_meat.hpp.
References arma_stop(), Mat< eT >::colptr(), syslib::copy_elem(), podarray< eT >::memptr(), Mat< eT >::memptr(), Mat< eT >::n_cols, Mat< eT >::n_rows, Mat< eT >::set_size(), trans(), and Mat< eT >::zeros().
Referenced by glue_solve::apply().
{ arma_extra_debug_sigprint(); #if defined(ARMA_USE_LAPACK) { char trans = 'N'; int m = A.n_rows; int n = A.n_cols; int lda = A.n_rows; int ldb = A.n_cols; int nrhs = B.n_cols; int lwork = m + (std::max)(m,nrhs); int info; Mat<eT> A_copy = A; Mat<eT> tmp; tmp.zeros(A.n_cols, B.n_cols); for(u32 col=0; col<B.n_cols; ++col) { eT* tmp_colmem = tmp.colptr(col); syslib::copy_elem( tmp_colmem, B.colptr(col), B.n_rows ); for(u32 row=B.n_rows; row<A.n_cols; ++row) { tmp_colmem[row] = eT(0); } } podarray<eT> work(lwork); arma_extra_debug_print("lapack::gels_()"); // NOTE: // the dgels() function in the lapack library supplied by ATLAS 3.6 // seems to have problems lapack::gels_<eT> ( &trans, &m, &n, &nrhs, A_copy.memptr(), &lda, tmp.memptr(), &ldb, work.memptr(), &lwork, &info ); arma_extra_debug_print("lapack::gels_() -- finished"); out.set_size(A.n_cols, B.n_cols); for(u32 col=0; col<B.n_cols; ++col) { syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols ); } return (info == 0); } #else { arma_stop("auxlib::solve_ud(): need LAPACK"); return false; } #endif }