Eigen::PardisoLDLT< MatrixType, Options > Class Template Reference

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

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

Inheritance diagram for Eigen::PardisoLDLT< MatrixType, Options >:

Collaboration diagram for Eigen::PardisoLDLT< MatrixType, Options >:

Public Types
enum	{ UpLo = Options&(Upper|Lower) }
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
	PardisoLDLT ()
	PardisoLDLT (const MatrixType &matrix)
void	getMatrix (const MatrixType &matrix)
Index	cols () const
Index	rows () const
ComputationInfo	info () const
	Reports whether previous computation was successful.
ParameterType &	pardisoParameterArray ()
PardisoLDLT< MatrixType, Options > &	analyzePattern (const MatrixType &matrix)
PardisoLDLT< MatrixType, Options > &	factorize (const MatrixType &matrix)
PardisoLDLT< MatrixType, Options > &	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< PardisoLDLT< MatrixType, Options > >	Base
typedef Base::Scalar	Scalar
typedef Base::RealScalar	RealScalar
typedef internal::pardiso_traits< PardisoLDLT< MatrixType, Options > >	Traits

Protected Member Functions
void	pardisoRelease ()
void	pardisoInit (int type)
void	manageErrorCode (Index error) const
PardisoLDLT< MatrixType, Options > &	derived ()
const PardisoLDLT< MatrixType, Options > &	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< PardisoLDLT< MatrixType, Options > >

Detailed Description

template<typename MatrixType, int Options>
class Eigen::PardisoLDLT< MatrixType, Options >

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

This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite. For complex matrices, A can also be symmetric only, see the Options template parameter. 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<>
Options	can be any bitwise combination of Upper, Lower, and Symmetric. The default is Upper, meaning only the upper triangular part has to be used. Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix. Upper|Lower can be used to tell both triangular parts can be used as input.

\implsparsesolverconcept

See also

TutorialSparseSolverConcept, class SimplicialLDLT

Member Typedef Documentation

Base

template<typename MatrixType , int Options>

typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Eigen::PardisoLDLT< MatrixType, Options >::Base

protected

IntColVectorType

typedef Matrix<StorageIndex, MatrixType::RowsAtCompileTime, 1> Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::IntColVectorType

inherited

IntRowVectorType

typedef Matrix<StorageIndex, 1, MatrixType::ColsAtCompileTime> Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::IntRowVectorType

inherited

MatrixType

typedef Traits::MatrixType Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::MatrixType

inherited

ParameterType

typedef Array<StorageIndex,64,1,DontAlign> Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::ParameterType

inherited

RealScalar

template<typename MatrixType , int Options>

typedef Base::RealScalar Eigen::PardisoLDLT< MatrixType, Options >::RealScalar

protected

Scalar

template<typename MatrixType , int Options>

typedef Base::Scalar Eigen::PardisoLDLT< MatrixType, Options >::Scalar

protected

SparseMatrixType

typedef SparseMatrix<Scalar,RowMajor,StorageIndex> Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::SparseMatrixType

inherited

StorageIndex

template<typename MatrixType , int Options>

typedef Base::StorageIndex Eigen::PardisoLDLT< MatrixType, Options >::StorageIndex

Traits

typedef internal::pardiso_traits<PardisoLDLT< MatrixType, Options > > Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::Traits

protectedinherited

VectorType

typedef Matrix<Scalar,Dynamic,1> Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::VectorType

inherited

Member Enumeration Documentation

anonymous enum

anonymous enum

inherited

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

anonymous enum

template<typename MatrixType , int Options>

anonymous enum

Enumerator
UpLo

{ UpLo = Options&(Upper|Lower) };

Constructor & Destructor Documentation

PardisoLDLT() [1/2]

template<typename MatrixType , int Options>

Eigen::PardisoLDLT< MatrixType, Options >::PardisoLDLT ( )

inline

      : Base()
    {
      pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
    }

References Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::pardisoInit(), Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::ScalarIsComplex, and Eigen::Symmetric.

Here is the call graph for this function:

PardisoLDLT() [2/2]

template<typename MatrixType , int Options>

Eigen::PardisoLDLT< MatrixType, Options >::PardisoLDLT ( const MatrixType &  matrix )

inlineexplicit

      : Base()
    {
      pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
      compute(matrix);
    }

References Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::compute(), Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::pardisoInit(), Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::ScalarIsComplex, and Eigen::Symmetric.

Here is the call graph for this function:

Member Function Documentation

_solve_impl() [1/3]

void Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::_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< PardisoLDLT< MatrixType, Options > >::_solve_impl ( const MatrixBase< Rhs > &  b, MatrixBase< Dest > &  dest  ) const

inherited

_solve_impl() [3/3]

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

inlineinherited

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

analyzePattern()

PardisoLDLT< MatrixType, Options > & Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::cols ( ) const

inlineinherited

{ return m_size; }

compute()

PardisoLDLT< MatrixType, Options > & Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::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]

PardisoLDLT< MatrixType, Options > & Eigen::SparseSolverBase< PardisoLDLT< MatrixType, Options > >::derived ( )

inlineprotectedinherited

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

derived() [2/2]

const PardisoLDLT< MatrixType, Options > & Eigen::SparseSolverBase< PardisoLDLT< MatrixType, Options > >::derived ( ) const

inlineprotectedinherited

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

factorize()

PardisoLDLT< MatrixType, Options > & Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::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 Options>

void Eigen::PardisoLDLT< MatrixType, Options >::getMatrix ( const MatrixType &  matrix )

inline

    {
      // 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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >::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< PardisoLDLT< MatrixType, Options > >

template<typename MatrixType , int Options>

friend class PardisoImpl< PardisoLDLT< MatrixType, Options > >

friend

Member Data Documentation

m_analysisIsOk

bool Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_analysisIsOk

protectedinherited

m_factorizationIsOk

bool Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_factorizationIsOk

protectedinherited

m_info

ComputationInfo Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_info

mutableprotectedinherited

m_iparm

ParameterType Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_iparm

mutableprotectedinherited

m_isInitialized

bool Eigen::SparseSolverBase< PardisoLDLT< MatrixType, Options > >::m_isInitialized

mutableprotectedinherited

m_matrix

SparseMatrixType Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_matrix

mutableprotectedinherited

m_msglvl

StorageIndex Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_msglvl

protectedinherited

m_perm

IntColVectorType Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_perm

mutableprotectedinherited

m_pt

void* Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_pt[64]

mutableprotectedinherited

m_size

Index Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_size

protectedinherited

m_type

StorageIndex Eigen::PardisoImpl< PardisoLDLT< MatrixType, Options > >::m_type

protectedinherited

---

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

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