## Cholesky decomposition problem with Intel MKL

kimsoohwan
Registered Member
Posts
7
Karma
0

### Cholesky decomposition problem with Intel MKL

Mon Aug 18, 2014 11:46 am
Hi,

First of all, thank you very much for this great library.
I have a problem when I do cholesky decomposition with Intel MKL. (Without Intel MKL, everything is fine.)
Even though the matrix is not singular, LLT.info() returns Eigen::NumericalIssue, but the calculation is still correct.
How can I check whether the matrix is singular or not with Intel MKL.
Thank you very much for your help.

Cheers,
S.

Code: Select all
`#include <iostream>#define EIGEN_NO_DEBUG      // to speed up#define EIGEN_USE_MKL_ALL   // to use Intel MKL#include <Eigen/Dense>int main(){   // dim   const int D = 4;   // matrix   Eigen::MatrixXf non_singular_mat(D, D);   non_singular_mat << 2.250000000000000f, 1.172075789902540f, 0.079002611998964f, 0.300544591865789f,                       1.172075789902540f, 2.250000000000000f, 0.499333922441860f, 0.352442509452200f,                       0.079002611998964f, 0.499333922441860f, 2.250000000000000f, 0.455259868374894f,                       0.300544591865789f, 0.352442509452200f, 0.455259868374894f, 2.250000000000000f;   // cholesky decomposition   Eigen::LLT<Eigen::MatrixXf> L(non_singular_mat);   if(L.info() != Eigen::Success)   {      switch(L.info())      {         case Eigen::NumericalIssue :         {            std::cout << "NumericalIssue" << std::endl; // With Intel MKL (EIGEN_USE_MKL_ALL), this flag is on!!!            break;         }         case Eigen::NoConvergence :         {            std::cout << "NoConvergence" << std::endl;            break;         }         case Eigen::InvalidInput :         {            std::cout << "InvalidInput" << std::endl;            break;         }      }   }   // inverse   Eigen::MatrixXf inv_non_singular_mat = L.solve(Eigen::MatrixXf::Identity(D, D));   // (LL')*inv(Cov) = I   // I   if((non_singular_mat * inv_non_singular_mat).isApprox(Eigen::MatrixXf::Identity(D, D)))       std::cout << "Non Singular: OK" << std::endl; // The calculation is correct.   else      std::cout << "Non Singular: Error!!!" << std::endl;    return 0;}`
ggael
Moderator
Posts
3075
Karma
18
OS

### Re: Cholesky decomposition problem with Intel MKL

Thu Aug 21, 2014 8:55 pm
I don't have MKL right now to debug that, but you can easily explore what's going one by looking at the file Eigen/src/Cholesky/LLT_MKL.h, and printing out the value of 'info' as returned by spotrf. See here for its meaning: http://www.netlib.org/lapack/explore-3. ... trf.f.html.
kimsoohwan
Registered Member
Posts
7
Karma
0

### Re: Cholesky decomposition problem with Intel MKL

Sat Aug 23, 2014 8:16 am
Hi ggal,

I've already checked the info variable in the Eigen::LLT class.
But it was returned as a non-zero value (which means the matrix is singular) even for a non-singular matrix.
I tried to parse the Eigen::LLT class, but that was not easy to debug. (Probably it is not a bug.)
So, I implemented my LLT class and would like to leave here for others as below.
It will work was Eigen::LLT if the flags for Intel MKL are not on, but replace the Eigen::compute function otherwise.

Cheers

Code: Select all
`/** @class   LLT_MKL  * @brief   Wrapper of LLT for using Intel MKL  * @note   Refer to Eigen::LLT  */template<typename _MatrixType, int _UpLo = Eigen::Lower>class LLT_MKL : public Eigen::LLT<_MatrixType, _UpLo>{public:   LLT_MKL(const MatrixType &A)      : Eigen::LLT<_MatrixType, _UpLo>()   {      compute(A);   }   /**@ brief Cholesky decomposition */   Eigen::LLT<_MatrixType, _UpLo>& compute(const MatrixType& A)   {#if defined(EIGEN_USE_LAPACKE) || defined(EIGEN_USE_LAPACKE_STRICT)      // resize      eigen_assert(A.rows() == A.cols());      const Index size = A.rows();      m_matrix.resize(size, size);      m_matrix = A;      m_isInitialized = true;      // set up parameters for ?potrf      const int StorageOrder = MatrixType::Flags & Eigen::RowMajorBit ? Eigen::RowMajor : Eigen::ColMajor;      const int matrix_order = (StorageOrder==Eigen::RowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR;      // lhs      const lapack_int n = m_matrix.rows();      const char uplo = (_UpLo == Eigen::Lower)   ? 'L' : 'U';      MatrixType::Scalar* a = &(m_matrix.coeffRef(0,0));      const lapack_int lda = m_matrix.outerStride();      // potrf      const lapack_int info = potrf<MatrixType::Scalar>(matrix_order, uplo, n, a, lda);      m_info = (info==0) ? Eigen::Success : Eigen::NumericalIssue;      return static_cast<Eigen::LLT<_MatrixType, _UpLo> >(*this);#else      return Eigen::LLT<_MatrixType, _UpLo>::compute(A);#endif   }protected:#if defined(EIGEN_USE_LAPACKE) || defined(EIGEN_USE_LAPACKE_STRICT)   template <typename T>   inline lapack_int potrf(int matrix_order, char uplo, lapack_int n, T* a, lapack_int lda);   template <>   inline lapack_int potrf<float>(int matrix_order, char uplo, lapack_int n, float* a, lapack_int lda)   {      return LAPACKE_spotrf(matrix_order, uplo, n, a, lda);   }   template <>   inline lapack_int potrf<double>(int matrix_order, char uplo, lapack_int n, double* a, lapack_int lda)   {      return LAPACKE_dpotrf(matrix_order, uplo, n, a, lda);   }#endif};`

## Who is online

Registered users: Baidu [Spider], Bing [Bot], Exabot [Bot], Google [Bot], Sogou [Bot], Yahoo [Bot]