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.
Please see the following code.
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;
}


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.

Hi ggal,

Thank you for your reply.
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
};