Eigen::PardisoLLT< MatrixType, _UpLo > Class Template Reference

A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library. More...

#include <src/eigen/Eigen/src/PardisoSupport/PardisoSupport.h>

Inheritance diagram for Eigen::PardisoLLT< MatrixType, _UpLo >:

Collaboration diagram for Eigen::PardisoLLT< MatrixType, _UpLo >:

Public Types
enum	{ UpLo = _UpLo }
typedef Base::StorageIndex	StorageIndex
enum
typedef Traits::MatrixType	MatrixType
typedef SparseMatrix< Scalar, RowMajor, StorageIndex >	SparseMatrixType
typedef Matrix< Scalar, Dynamic, 1 >	VectorType
typedef Matrix< StorageIndex, 1, MatrixType::ColsAtCompileTime >	IntRowVectorType
typedef Matrix< StorageIndex, MatrixType::RowsAtCompileTime, 1 >	IntColVectorType
typedef Array< StorageIndex, 64, 1, DontAlign >	ParameterType

Public Member Functions
	PardisoLLT ()
	PardisoLLT (const MatrixType &matrix)
Index	cols () const
Index	rows () const
ComputationInfo	info () const
	Reports whether previous computation was successful.
ParameterType &	pardisoParameterArray ()
PardisoLLT< MatrixType, _UpLo > &	analyzePattern (const MatrixType &matrix)
PardisoLLT< MatrixType, _UpLo > &	factorize (const MatrixType &matrix)
PardisoLLT< MatrixType, _UpLo > &	compute (const MatrixType &matrix)
void	_solve_impl (const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
void	_solve_impl (const MatrixBase< BDerived > &b, MatrixBase< XDerived > &x) const
void	_solve_impl (const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
template<typename Rhs >
const Solve< Derived, Rhs >	solve (const MatrixBase< Rhs > &b) const
template<typename Rhs >
const Solve< Derived, Rhs >	solve (const SparseMatrixBase< Rhs > &b) const

Protected Types
typedef PardisoImpl< PardisoLLT< MatrixType, _UpLo > >	Base
typedef Base::Scalar	Scalar
typedef Base::RealScalar	RealScalar
typedef internal::pardiso_traits< PardisoLLT< MatrixType, _UpLo > >	Traits

Protected Member Functions
void	getMatrix (const MatrixType &matrix)
void	pardisoRelease ()
void	pardisoInit (int type)
void	manageErrorCode (Index error) const
PardisoLLT< MatrixType, _UpLo > &	derived ()
const PardisoLLT< MatrixType, _UpLo > &	derived () const

Protected Attributes
SparseMatrixType	m_matrix
ComputationInfo	m_info
bool	m_analysisIsOk
bool	m_factorizationIsOk
StorageIndex	m_type
StorageIndex	m_msglvl
void *	m_pt [64]
ParameterType	m_iparm
IntColVectorType	m_perm
Index	m_size
bool	m_isInitialized

Friends
class	PardisoImpl< PardisoLLT< MatrixType, _UpLo > >

Detailed Description

template<typename MatrixType, int _UpLo>
class Eigen::PardisoLLT< MatrixType, _UpLo >

A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library.

This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite. The vectors or matrices X and B can be either dense or sparse.

By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set:

solver.pardisoParameterArray()[59] = 1; 

Template Parameters

MatrixType	the type of the sparse matrix A, it must be a SparseMatrix<>
UpLo	can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used. Upper|Lower can be used to tell both triangular parts can be used as input.

\implsparsesolverconcept

See also

TutorialSparseSolverConcept, class SimplicialLLT

Member Typedef Documentation

Base

template<typename MatrixType , int _UpLo>

typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Eigen::PardisoLLT< MatrixType, _UpLo >::Base

protected

IntColVectorType

typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::IntColVectorType

inherited

IntRowVectorType

typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::IntRowVectorType

inherited

MatrixType

typedef Traits::MatrixType Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::MatrixType

inherited

ParameterType

typedef Array<StorageIndex,64,1,DontAlign> Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::ParameterType

inherited

RealScalar

template<typename MatrixType , int _UpLo>

typedef Base::RealScalar Eigen::PardisoLLT< MatrixType, _UpLo >::RealScalar

protected

Scalar

template<typename MatrixType , int _UpLo>

typedef Base::Scalar Eigen::PardisoLLT< MatrixType, _UpLo >::Scalar

protected

SparseMatrixType

typedef SparseMatrix<Scalar,RowMajor,StorageIndex> Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::SparseMatrixType

inherited

StorageIndex

template<typename MatrixType , int _UpLo>

typedef Base::StorageIndex Eigen::PardisoLLT< MatrixType, _UpLo >::StorageIndex

Traits

typedef internal::pardiso_traits<PardisoLLT< MatrixType, _UpLo > > Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::Traits

protectedinherited

VectorType

typedef Matrix<Scalar,Dynamic,1> Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::VectorType

inherited

Member Enumeration Documentation

anonymous enum

anonymous enum

inherited

         {
      ScalarIsComplex = NumTraits<Scalar>::IsComplex,
      ColsAtCompileTime = Dynamic,
      MaxColsAtCompileTime = Dynamic
    };

anonymous enum

template<typename MatrixType , int _UpLo>

anonymous enum

Enumerator
UpLo

{ UpLo = _UpLo };

Constructor & Destructor Documentation

PardisoLLT() [1/2]

template<typename MatrixType , int _UpLo>

Eigen::PardisoLLT< MatrixType, _UpLo >::PardisoLLT ( )

inline

      : Base()
    {
      pardisoInit(Base::ScalarIsComplex ? 4 : 2);
    }

References Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::pardisoInit(), and Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::ScalarIsComplex.

Here is the call graph for this function:

PardisoLLT() [2/2]

template<typename MatrixType , int _UpLo>

Eigen::PardisoLLT< MatrixType, _UpLo >::PardisoLLT ( const MatrixType &  matrix )

inlineexplicit

      : Base()
    {
      pardisoInit(Base::ScalarIsComplex ? 4 : 2);
      compute(matrix);
    }

References Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::compute(), Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::pardisoInit(), and Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::ScalarIsComplex.

Here is the call graph for this function:

Member Function Documentation

_solve_impl() [1/3]

void Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::_solve_impl ( const MatrixBase< BDerived > &  b, MatrixBase< XDerived > &  x  ) const

inherited

{
  if(m_iparm[0] == 0) // Factorization was not computed
  {
    m_info = InvalidInput;
    return;
  }
 
  //Index n = m_matrix.rows();
  Index nrhs = Index(b.cols());
  eigen_assert(m_size==b.rows());
  eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
  eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
  eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
 
 
//  switch (transposed) {
//    case SvNoTrans    : m_iparm[11] = 0 ; break;
//    case SvTranspose  : m_iparm[11] = 2 ; break;
//    case SvAdjoint    : m_iparm[11] = 1 ; break;
//    default:
//      //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the PARDISO backend\n";
//      m_iparm[11] = 0;
//  }
 
  Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
  Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
  
  // Pardiso cannot solve in-place
  if(rhs_ptr == x.derived().data())
  {
    tmp = b;
    rhs_ptr = tmp.data();
  }
  
  Index error;
  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 33, internal::convert_index<StorageIndex>(m_size),
                                                            m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
                                                            m_perm.data(), internal::convert_index<StorageIndex>(nrhs), m_iparm.data(), m_msglvl,
                                                            rhs_ptr, x.derived().data());
 
  manageErrorCode(error);
}

_solve_impl() [2/3]

void Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::_solve_impl ( const MatrixBase< Rhs > &  b, MatrixBase< Dest > &  dest  ) const

inherited

_solve_impl() [3/3]

void Eigen::SparseSolverBase< PardisoLLT< MatrixType, _UpLo > >::_solve_impl ( const SparseMatrixBase< Rhs > &  b, SparseMatrixBase< Dest > &  dest  ) const

inlineinherited

    {
      internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived());
    }

analyzePattern()

PardisoLLT< MatrixType, _UpLo > & Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::analyzePattern ( const MatrixType &  matrix )

inherited

Performs a symbolic decomposition on the sparcity of matrix.

This function is particularly useful when solving for several problems having the same structure.

See also

factorize()

{
  m_size = a.rows();
  eigen_assert(m_size == a.cols());
 
  pardisoRelease();
  m_perm.setZero(m_size);
  derived().getMatrix(a);
  
  Index error;
  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 11, internal::convert_index<StorageIndex>(m_size),
                                                            m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
                                                            m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
  
  manageErrorCode(error);
  m_analysisIsOk = true;
  m_factorizationIsOk = false;
  m_isInitialized = true;
  return derived();
}

cols()

Index Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::cols ( ) const

inlineinherited

{ return m_size; }

compute()

PardisoLLT< MatrixType, _UpLo > & Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::compute ( const MatrixType &  matrix )

inherited

{
  m_size = a.rows();
  eigen_assert(a.rows() == a.cols());
 
  pardisoRelease();
  m_perm.setZero(m_size);
  derived().getMatrix(a);
  
  Index error;
  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 12, internal::convert_index<StorageIndex>(m_size),
                                                            m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
                                                            m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
  manageErrorCode(error);
  m_analysisIsOk = true;
  m_factorizationIsOk = true;
  m_isInitialized = true;
  return derived();
}

derived() [1/2]

PardisoLLT< MatrixType, _UpLo > & Eigen::SparseSolverBase< PardisoLLT< MatrixType, _UpLo > >::derived ( )

inlineprotectedinherited

{ return *static_cast<Derived*>(this); }

derived() [2/2]

const PardisoLLT< MatrixType, _UpLo > & Eigen::SparseSolverBase< PardisoLLT< MatrixType, _UpLo > >::derived ( ) const

inlineprotectedinherited

{ return *static_cast<const Derived*>(this); }

factorize()

PardisoLLT< MatrixType, _UpLo > & Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::factorize ( const MatrixType &  matrix )

inherited

Performs a numeric decomposition of matrix

The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.

See also

analyzePattern()

{
  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
  eigen_assert(m_size == a.rows() && m_size == a.cols());
  
  derived().getMatrix(a);
 
  Index error;
  error = internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, 22, internal::convert_index<StorageIndex>(m_size),
                                                            m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
                                                            m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
  
  manageErrorCode(error);
  m_factorizationIsOk = true;
  return derived();
}

getMatrix()

template<typename MatrixType , int _UpLo>

void Eigen::PardisoLLT< MatrixType, _UpLo >::getMatrix ( const MatrixType &  matrix )

inlineprotected

    {
      // PARDISO supports only upper, row-major matrices
      PermutationMatrix<Dynamic,Dynamic,StorageIndex> p_null;
      m_matrix.resize(matrix.rows(), matrix.cols());
      m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
      m_matrix.makeCompressed();
    }

References Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_matrix, Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::makeCompressed(), Eigen::SparseMatrix< _Scalar, _Options, _StorageIndex >::resize(), and Eigen::SparseMatrixBase< Derived >::twistedBy().

Here is the call graph for this function:

info()

ComputationInfo Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::info ( ) const

inlineinherited

Reports whether previous computation was successful.

Returns

Success if computation was succesful, NumericalIssue if the matrix appears to be negative.

    {
      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
      return m_info;
    }

manageErrorCode()

void Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::manageErrorCode ( Index  error ) const

inlineprotectedinherited

    {
      switch(error)
      {
        case 0:
          m_info = Success;
          break;
        case -4:
        case -7:
          m_info = NumericalIssue;
          break;
        default:
          m_info = InvalidInput;
      }
    }

pardisoInit()

void Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::pardisoInit ( int  type )

inlineprotectedinherited

    {
      m_type = type;
      bool symmetric = std::abs(m_type) < 10;
      m_iparm[0] = 1;   // No solver default
      m_iparm[1] = 2;   // use Metis for the ordering
      m_iparm[2] = 0;   // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
      m_iparm[3] = 0;   // No iterative-direct algorithm
      m_iparm[4] = 0;   // No user fill-in reducing permutation
      m_iparm[5] = 0;   // Write solution into x, b is left unchanged
      m_iparm[6] = 0;   // Not in use
      m_iparm[7] = 2;   // Max numbers of iterative refinement steps
      m_iparm[8] = 0;   // Not in use
      m_iparm[9] = 13;  // Perturb the pivot elements with 1E-13
      m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
      m_iparm[11] = 0;  // Not in use
      m_iparm[12] = symmetric ? 0 : 1;  // Maximum weighted matching algorithm is switched-off (default for symmetric).
                                        // Try m_iparm[12] = 1 in case of inappropriate accuracy
      m_iparm[13] = 0;  // Output: Number of perturbed pivots
      m_iparm[14] = 0;  // Not in use
      m_iparm[15] = 0;  // Not in use
      m_iparm[16] = 0;  // Not in use
      m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
      m_iparm[18] = -1; // Output: Mflops for LU factorization
      m_iparm[19] = 0;  // Output: Numbers of CG Iterations
      
      m_iparm[20] = 0;  // 1x1 pivoting
      m_iparm[26] = 0;  // No matrix checker
      m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
      m_iparm[34] = 1;  // C indexing
      m_iparm[36] = 0;  // CSR
      m_iparm[59] = 0;  // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
      
      memset(m_pt, 0, sizeof(m_pt));
    }

pardisoParameterArray()

ParameterType & Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::pardisoParameterArray ( )

inlineinherited

Warning

for advanced usage only.

Returns

a reference to the parameter array controlling PARDISO. See the PARDISO manual to know how to use it.

    {
      return m_iparm;
    }

pardisoRelease()

void Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::pardisoRelease ( )

inlineprotectedinherited

    {
      if(m_isInitialized) // Factorization ran at least once
      {
        internal::pardiso_run_selector<StorageIndex>::run(m_pt, 1, 1, m_type, -1, internal::convert_index<StorageIndex>(m_size),0, 0, 0, m_perm.data(), 0,
                                                          m_iparm.data(), m_msglvl, NULL, NULL);
        m_isInitialized = false;
      }
    }

rows()

Index Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::rows ( ) const

inlineinherited

{ return m_size; }

solve() [1/2]

template<typename Derived >

template<typename Rhs >

const Solve< Derived, Rhs > Eigen::SparseSolverBase< Derived >::solve ( const MatrixBase< Rhs > &  b ) const

inlineinherited

Returns

an expression of the solution x of using the current decomposition of A.

See also

compute()

    {
      eigen_assert(m_isInitialized && "Solver is not initialized.");
      eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
      return Solve<Derived, Rhs>(derived(), b.derived());
    }

References Eigen::SparseSolverBase< Derived >::derived(), eigen_assert, and Eigen::SparseSolverBase< Derived >::m_isInitialized.

Referenced by igl::embree::bone_heat(), igl::Frame_field_deformer::compute_optimal_positions(), igl::eigs(), igl::copyleft::comiso::FrameInterpolator::interpolateSymmetric(), igl::slim::solve_weighted_arap(), and igl::copyleft::comiso::NRosyField::solveNoRoundings().

Here is the call graph for this function:

Here is the caller graph for this function:

solve() [2/2]

template<typename Derived >

template<typename Rhs >

const Solve< Derived, Rhs > Eigen::SparseSolverBase< Derived >::solve ( const SparseMatrixBase< Rhs > &  b ) const

inlineinherited

Returns

an expression of the solution x of using the current decomposition of A.

See also

compute()

    {
      eigen_assert(m_isInitialized && "Solver is not initialized.");
      eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b");
      return Solve<Derived, Rhs>(derived(), b.derived());
    }

Friends And Related Symbol Documentation

PardisoImpl< PardisoLLT< MatrixType, _UpLo > >

template<typename MatrixType , int _UpLo>

friend class PardisoImpl< PardisoLLT< MatrixType, _UpLo > >

friend

Member Data Documentation

m_analysisIsOk

bool Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_analysisIsOk

protectedinherited

m_factorizationIsOk

bool Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_factorizationIsOk

protectedinherited

m_info

ComputationInfo Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_info

mutableprotectedinherited

m_iparm

ParameterType Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_iparm

mutableprotectedinherited

m_isInitialized

bool Eigen::SparseSolverBase< PardisoLLT< MatrixType, _UpLo > >::m_isInitialized

mutableprotectedinherited

m_matrix

SparseMatrixType Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_matrix

mutableprotectedinherited

m_msglvl

StorageIndex Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_msglvl

protectedinherited

m_perm

IntColVectorType Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_perm

mutableprotectedinherited

m_pt

void* Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_pt[64]

mutableprotectedinherited

m_size

Index Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_size

protectedinherited

m_type

StorageIndex Eigen::PardisoImpl< PardisoLLT< MatrixType, _UpLo > >::m_type

protectedinherited

---

The documentation for this class was generated from the following file:

• src/eigen/Eigen/src/PardisoSupport/PardisoSupport.h