#include <PardisoSupport.h>

Inheritance diagram for Eigen::PardisoImpl< Derived >:

Collaboration diagram for Eigen::PardisoImpl< Derived >:

Public Types
enum	{ ScalarIsComplex = NumTraits<Scalar>::IsComplex, ColsAtCompileTime = Dynamic, MaxColsAtCompileTime = Dynamic }
typedef Traits::MatrixType	MatrixType
typedef Traits::Scalar	Scalar
typedef Traits::RealScalar	RealScalar
typedef Traits::StorageIndex	StorageIndex
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
	PardisoImpl ()
	~PardisoImpl ()
Index	cols () const
Index	rows () const
ComputationInfo	info () const
	Reports whether previous computation was successful.
ParameterType &	pardisoParameterArray ()
Derived &	analyzePattern (const MatrixType &matrix)
Derived &	factorize (const MatrixType &matrix)
Derived &	compute (const MatrixType &matrix)
template<typename Rhs , typename Dest >
void	_solve_impl (const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
template<typename BDerived , typename XDerived >
void	_solve_impl (const MatrixBase< BDerived > &b, MatrixBase< XDerived > &x) const
Protected Types
typedef SparseSolverBase< Derived >	Base
typedef internal::pardiso_traits < Derived >	Traits
Protected Member Functions
void	pardisoRelease ()
void	pardisoInit (int type)
void	manageErrorCode (Index error) const
Protected Attributes
SparseMatrixType	m_matrix
ComputationInfo	m_info
bool	m_analysisIsOk
bool	m_factorizationIsOk
Index	m_type
Index	m_msglvl
void *	m_pt [64]
ParameterType	m_iparm
IntColVectorType	m_perm
Index	m_size

Detailed Description

template<class Derived>
class Eigen::PardisoImpl< Derived >

Definition at line 99 of file PardisoSupport.h.

Member Typedef Documentation

template<class Derived>

typedef SparseSolverBase<Derived> Eigen::PardisoImpl< Derived >::Base [protected]

Reimplemented in Eigen::PardisoLDLT< MatrixType, Options >, Eigen::PardisoLLT< MatrixType, _UpLo >, and Eigen::PardisoLU< MatrixType >.

Definition at line 102 of file PardisoSupport.h.

template<class Derived>

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

Definition at line 117 of file PardisoSupport.h.

template<class Derived>

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

Definition at line 116 of file PardisoSupport.h.

template<class Derived>

typedef Traits::MatrixType Eigen::PardisoImpl< Derived >::MatrixType

Definition at line 110 of file PardisoSupport.h.

template<class Derived>

typedef Array<StorageIndex,64,1,DontAlign> Eigen::PardisoImpl< Derived >::ParameterType

Definition at line 118 of file PardisoSupport.h.

template<class Derived>

typedef Traits::RealScalar Eigen::PardisoImpl< Derived >::RealScalar

Reimplemented in Eigen::PardisoLDLT< MatrixType, Options >, Eigen::PardisoLLT< MatrixType, _UpLo >, and Eigen::PardisoLU< MatrixType >.

Definition at line 112 of file PardisoSupport.h.

template<class Derived>

typedef Traits::Scalar Eigen::PardisoImpl< Derived >::Scalar

Reimplemented in Eigen::PardisoLDLT< MatrixType, Options >, Eigen::PardisoLLT< MatrixType, _UpLo >, and Eigen::PardisoLU< MatrixType >.

Definition at line 111 of file PardisoSupport.h.

template<class Derived>

typedef SparseMatrix<Scalar,RowMajor,StorageIndex> Eigen::PardisoImpl< Derived >::SparseMatrixType

Definition at line 114 of file PardisoSupport.h.

template<class Derived>

typedef Traits::StorageIndex Eigen::PardisoImpl< Derived >::StorageIndex

Reimplemented in Eigen::PardisoLDLT< MatrixType, Options >, and Eigen::PardisoLLT< MatrixType, _UpLo >.

Definition at line 113 of file PardisoSupport.h.

template<class Derived>

typedef internal::pardiso_traits<Derived> Eigen::PardisoImpl< Derived >::Traits [protected]

Definition at line 106 of file PardisoSupport.h.

template<class Derived>

typedef Matrix<Scalar,Dynamic,1> Eigen::PardisoImpl< Derived >::VectorType

Definition at line 115 of file PardisoSupport.h.

Member Enumeration Documentation

template<class Derived>

anonymous enum

Enumerator:

ScalarIsComplex
ColsAtCompileTime
MaxColsAtCompileTime

Definition at line 119 of file PardisoSupport.h.

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

Constructor & Destructor Documentation

template<class Derived>

Eigen::PardisoImpl< Derived >::PardisoImpl ( ) [inline]

Definition at line 125 of file PardisoSupport.h.

    {
      eigen_assert((sizeof(StorageIndex) >= sizeof(_INTEGER_t) && sizeof(StorageIndex) <= 8) && "Non-supported index type");
      m_iparm.setZero();
      m_msglvl = 0; // No output
      m_isInitialized = false;
    }

template<class Derived>

Eigen::PardisoImpl< Derived >::~PardisoImpl ( ) [inline]

Definition at line 133 of file PardisoSupport.h.

    {
      pardisoRelease();
    }

Member Function Documentation

template<class Derived>

template<typename Rhs , typename Dest >

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

template<class Derived>

template<typename BDerived , typename XDerived >

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

Definition at line 321 of file PardisoSupport.h.

{
  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, m_size,
                                                            m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
                                                            m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
                                                            rhs_ptr, x.derived().data());

  manageErrorCode(error);
}

template<class Derived >

Derived & Eigen::PardisoImpl< Derived >::analyzePattern ( const MatrixType & matrix )

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()

Definition at line 280 of file PardisoSupport.h.

{
  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, 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();
}

template<class Derived>

Index Eigen::PardisoImpl< Derived >::cols ( void ) const [inline]

Definition at line 138 of file PardisoSupport.h.

{ return m_size; }

template<class Derived >

Derived & Eigen::PardisoImpl< Derived >::compute ( const MatrixType & matrix )

Definition at line 258 of file PardisoSupport.h.

{
  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, 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();
}

template<class Derived >

Derived & Eigen::PardisoImpl< Derived >::factorize ( const MatrixType & matrix )

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()

Definition at line 302 of file PardisoSupport.h.

{
  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, 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();
}

template<class Derived>

ComputationInfo Eigen::PardisoImpl< Derived >::info ( ) const [inline]

Reports whether previous computation was successful.

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

Definition at line 146 of file PardisoSupport.h.

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

template<class Derived>

void Eigen::PardisoImpl< Derived >::manageErrorCode ( Index error ) const [inline, protected]

Definition at line 230 of file PardisoSupport.h.

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

template<class Derived>

void Eigen::PardisoImpl< Derived >::pardisoInit ( int type ) [inline, protected]

Definition at line 192 of file PardisoSupport.h.

    {
      m_type = type;
      bool symmetric = std::abs(m_type) < 10;
      m_iparm[0] = 1;   // No solver default
      m_iparm[1] = 3;   // use Metis for the ordering
      m_iparm[2] = 1;   // 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
      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[59] = 1;  // Automatic switch between In-Core and Out-of-Core modes
      
      memset(m_pt, 0, sizeof(m_pt));
    }

template<class Derived>

ParameterType& Eigen::PardisoImpl< Derived >::pardisoParameterArray ( ) [inline]

Warning:: for advanced usage only.

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

Definition at line 155 of file PardisoSupport.h.

    {
      return m_iparm;
    }

template<class Derived>

void Eigen::PardisoImpl< Derived >::pardisoRelease ( ) [inline, protected]

Definition at line 182 of file PardisoSupport.h.

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