Sparse supernodal LU factorization for general matrices. More...

Inheritance diagram for Eigen::SparseLU< MatrixType_, OrderingType_ >:

Public Types
enum	{ ColsAtCompileTime , MaxColsAtCompileTime }

typedef internal::SparseLUImpl< Scalar, StorageIndex >	Base

typedef Matrix< StorageIndex, Dynamic, 1 >	IndexVector

typedef MatrixType_	MatrixType

typedef SparseMatrix< Scalar, ColMajor, StorageIndex >	NCMatrix

typedef OrderingType_	OrderingType

typedef PermutationMatrix< Dynamic, Dynamic, StorageIndex >	PermutationType

typedef MatrixType::RealScalar	RealScalar

typedef MatrixType::Scalar	Scalar

typedef Matrix< Scalar, Dynamic, 1 >	ScalarVector

typedef internal::MappedSuperNodalMatrix< Scalar, StorageIndex >	SCMatrix

typedef MatrixType::StorageIndex	StorageIndex

Public Member Functions
template<typename Rhs , typename Dest >
bool	_solve_impl (const MatrixBase< Rhs > &B, MatrixBase< Dest > &X_base) const

Scalar	absDeterminant ()

const SparseLUTransposeView< true, SparseLU< MatrixType_, OrderingType_ > >	adjoint ()

void	analyzePattern (const MatrixType &matrix)

Index	cols () const

const PermutationType &	colsPermutation () const

void	compute (const MatrixType &matrix)

Scalar	determinant ()

void	factorize (const MatrixType &matrix)

ComputationInfo	info () const
	Reports whether previous computation was successful. More...

void	isSymmetric (bool sym)

std::string	lastErrorMessage () const

Scalar	logAbsDeterminant () const

SparseLUMatrixLReturnType< SCMatrix >	matrixL () const

SparseLUMatrixUReturnType< SCMatrix, Map< SparseMatrix< Scalar, ColMajor, StorageIndex > > >	matrixU () const

Index	nnzL () const

Index	nnzU () const

Index	rows () const

const PermutationType &	rowsPermutation () const

void	setPivotThreshold (const RealScalar &thresh)

Scalar	signDeterminant ()

void	simplicialfactorize (const MatrixType &matrix)

template<typename Rhs >
const Solve< SparseLU, Rhs >	solve (const MatrixBase< Rhs > &B) const

	SparseLU ()

	SparseLU (const MatrixType &matrix)

const SparseLUTransposeView< false, SparseLU< MatrixType_, OrderingType_ > >	transpose ()

	~SparseLU ()

Public Member Functions inherited from Eigen::SparseSolverBase< SparseLU< MatrixType_, OrderingType_ > >
SparseLU< MatrixType_, OrderingType_ > &	derived ()

const SparseLU< MatrixType_, OrderingType_ > &	derived () const

const Solve< SparseLU< MatrixType_, OrderingType_ >, Rhs >	solve (const MatrixBase< Rhs > &b) const

const Solve< SparseLU< MatrixType_, OrderingType_ >, Rhs >	solve (const SparseMatrixBase< Rhs > &b) const

	SparseSolverBase ()

	SparseSolverBase (SparseSolverBase &&other)

	~SparseSolverBase ()

Protected Types
typedef SparseSolverBase< SparseLU< MatrixType_, OrderingType_ > >	APIBase

Protected Member Functions
void	initperfvalues ()

Protected Attributes
bool	m_analysisIsOk

Index	m_detPermC

Index	m_detPermR

RealScalar	m_diagpivotthresh

IndexVector	m_etree

bool	m_factorizationIsOk

Base::GlobalLU_t	m_glu

ComputationInfo	m_info

std::string	m_lastError

SCMatrix	m_Lstore

NCMatrix	m_mat

Index	m_nnzL

Index	m_nnzU

internal::perfvalues	m_perfv

PermutationType	m_perm_c

PermutationType	m_perm_r

bool	m_symmetricmode

Map< SparseMatrix< Scalar, ColMajor, StorageIndex > >	m_Ustore

Protected Attributes inherited from Eigen::SparseSolverBase< SparseLU< MatrixType_, OrderingType_ > >
bool	m_isInitialized

Private Member Functions
	SparseLU (const SparseLU &)

Detailed Description

template<typename MatrixType_, typename OrderingType_>
class Eigen::SparseLU< MatrixType_, OrderingType_ >

Sparse supernodal LU factorization for general matrices.

This class implements the supernodal LU factorization for general matrices. It uses the main techniques from the sequential SuperLU package (http://crd-legacy.lbl.gov/~xiaoye/SuperLU/). It handles transparently real and complex arithmetic with single and double precision, depending on the scalar type of your input matrix. The code has been optimized to provide BLAS-3 operations during supernode-panel updates. It benefits directly from the built-in high-performant Eigen BLAS routines. Moreover, when the size of a supernode is very small, the BLAS calls are avoided to enable a better optimization from the compiler. For best performance, you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.

An important parameter of this class is the ordering method. It is used to reorder the columns (and eventually the rows) of the matrix to reduce the number of new elements that are created during numerical factorization. The cheapest method available is COLAMD. See the OrderingMethods module for the list of built-in and external ordering methods.

Simple example with key steps

VectorXd x(n), b(n);
SparseMatrix<double> A;
SparseLU<SparseMatrix<double>, COLAMDOrdering<int> >   solver;
// fill A and b;
// Compute the ordering permutation vector from the structural pattern of A
solver.analyzePattern(A); 
// Compute the numerical factorization 
solver.factorize(A); 
//Use the factors to solve the linear system 
x = solver.solve(b); 

Warning: The input matrix A should be in a compressed and column-major form. Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.

Note: Unlike the initial SuperLU implementation, there is no step to equilibrate the matrix. For badly scaled matrices, this step can be useful to reduce the pivoting during factorization. If this is the case for your matrices, you can try the basic scaling method at "unsupported/Eigen/src/IterativeSolvers/Scaling.h"

Template Parameters

MatrixType_	The type of the sparse matrix. It must be a column-major SparseMatrix<>
OrderingType_	The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD

This class follows the sparse solver concept .

See also: Sparse solver concept; OrderingMethods module

Definition at line 134 of file SparseLU.h.

Member Typedef Documentation

◆ APIBase

template<typename MatrixType_ , typename OrderingType_ >

typedef SparseSolverBase<SparseLU<MatrixType_,OrderingType_> > Eigen::SparseLU< MatrixType_, OrderingType_ >::APIBase

protected

Definition at line 137 of file SparseLU.h.

◆ Base

template<typename MatrixType_ , typename OrderingType_ >

typedef internal::SparseLUImpl<Scalar, StorageIndex> Eigen::SparseLU< MatrixType_, OrderingType_ >::Base

Definition at line 152 of file SparseLU.h.

◆ IndexVector

template<typename MatrixType_ , typename OrderingType_ >

typedef Matrix<StorageIndex,Dynamic,1> Eigen::SparseLU< MatrixType_, OrderingType_ >::IndexVector

Definition at line 150 of file SparseLU.h.

◆ MatrixType

template<typename MatrixType_ , typename OrderingType_ >

typedef MatrixType_ Eigen::SparseLU< MatrixType_, OrderingType_ >::MatrixType

Definition at line 142 of file SparseLU.h.

◆ NCMatrix

template<typename MatrixType_ , typename OrderingType_ >

typedef SparseMatrix<Scalar,ColMajor,StorageIndex> Eigen::SparseLU< MatrixType_, OrderingType_ >::NCMatrix

Definition at line 147 of file SparseLU.h.

◆ OrderingType

template<typename MatrixType_ , typename OrderingType_ >

typedef OrderingType_ Eigen::SparseLU< MatrixType_, OrderingType_ >::OrderingType

Definition at line 143 of file SparseLU.h.

◆ PermutationType

template<typename MatrixType_ , typename OrderingType_ >

typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> Eigen::SparseLU< MatrixType_, OrderingType_ >::PermutationType

Definition at line 151 of file SparseLU.h.

◆ RealScalar

template<typename MatrixType_ , typename OrderingType_ >

typedef MatrixType::RealScalar Eigen::SparseLU< MatrixType_, OrderingType_ >::RealScalar

Definition at line 145 of file SparseLU.h.

◆ Scalar

template<typename MatrixType_ , typename OrderingType_ >

typedef MatrixType::Scalar Eigen::SparseLU< MatrixType_, OrderingType_ >::Scalar

Definition at line 144 of file SparseLU.h.

◆ ScalarVector

template<typename MatrixType_ , typename OrderingType_ >

typedef Matrix<Scalar,Dynamic,1> Eigen::SparseLU< MatrixType_, OrderingType_ >::ScalarVector

Definition at line 149 of file SparseLU.h.

◆ SCMatrix

template<typename MatrixType_ , typename OrderingType_ >

typedef internal::MappedSuperNodalMatrix<Scalar, StorageIndex> Eigen::SparseLU< MatrixType_, OrderingType_ >::SCMatrix

Definition at line 148 of file SparseLU.h.

◆ StorageIndex

template<typename MatrixType_ , typename OrderingType_ >

typedef MatrixType::StorageIndex Eigen::SparseLU< MatrixType_, OrderingType_ >::StorageIndex

Definition at line 146 of file SparseLU.h.

Member Enumeration Documentation

◆ anonymous enum

template<typename MatrixType_ , typename OrderingType_ >

anonymous enum

Enumerator
ColsAtCompileTime
MaxColsAtCompileTime

Definition at line 154 of file SparseLU.h.

          {
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };

Constructor & Destructor Documentation

◆ SparseLU() [1/3]

template<typename MatrixType_ , typename OrderingType_ >

Eigen::SparseLU< MatrixType_, OrderingType_ >::SparseLU ( )

inline

Definition at line 161 of file SparseLU.h.

               :m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
     {
       initperfvalues(); 
     }

◆ SparseLU() [2/3]

template<typename MatrixType_ , typename OrderingType_ >

Eigen::SparseLU< MatrixType_, OrderingType_ >::SparseLU ( const MatrixType & matrix )

inlineexplicit

Definition at line 165 of file SparseLU.h.

       : m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
     {
       initperfvalues(); 
       compute(matrix);
     }

◆ ~SparseLU()

template<typename MatrixType_ , typename OrderingType_ >

Eigen::SparseLU< MatrixType_, OrderingType_ >::~SparseLU ( )

inline

Definition at line 172 of file SparseLU.h.

     {
       // Free all explicit dynamic pointers 
     }

◆ SparseLU() [3/3]

template<typename MatrixType_ , typename OrderingType_ >

Eigen::SparseLU< MatrixType_, OrderingType_ >::SparseLU ( const SparseLU< MatrixType_, OrderingType_ > & )

private

Member Function Documentation

◆ _solve_impl()

template<typename MatrixType_ , typename OrderingType_ >

template<typename Rhs , typename Dest >

bool Eigen::SparseLU< MatrixType_, OrderingType_ >::_solve_impl	(	const MatrixBase< Rhs > &	B,
		MatrixBase< Dest > &	X_base
	)		const

inline

Definition at line 317 of file SparseLU.h.

     {
       Dest& X(X_base.derived());
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
       EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
                         THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
       
       // Permute the right hand side to form X = Pr*B
       // on return, X is overwritten by the computed solution
       X.resize(B.rows(),B.cols());
  
       // this ugly const_cast_derived() helps to detect aliasing when applying the permutations
       for(Index j = 0; j < B.cols(); ++j)
         X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
       
       //Forward substitution with L
       this->matrixL().solveInPlace(X);
       this->matrixU().solveInPlace(X);
       
       // Permute back the solution 
       for (Index j = 0; j < B.cols(); ++j)
         X.col(j) = colsPermutation().inverse() * X.col(j);
       
       return true; 
     }

◆ absDeterminant()

template<typename MatrixType_ , typename OrderingType_ >

Scalar Eigen::SparseLU< MatrixType_, OrderingType_ >::absDeterminant ( )

inline

Returns: the absolute value of the determinant of the matrix of which *this is the QR decomposition.

Warning: a determinant can be very big or small, so for matrices of large enough dimension, there is a risk of overflow/underflow. One way to work around that is to use logAbsDeterminant() instead.

See also: logAbsDeterminant(), signDeterminant()

Definition at line 353 of file SparseLU.h.

     {
       using std::abs;
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       // Initialize with the determinant of the row matrix
       Scalar det = Scalar(1.);
       // Note that the diagonal blocks of U are stored in supernodes,
       // which are available in the  L part :)
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
           if(it.index() == j)
           {
             det *= abs(it.value());
             break;
           }
         }
       }
       return det;
     }

◆ adjoint()

template<typename MatrixType_ , typename OrderingType_ >

const SparseLUTransposeView<true, SparseLU<MatrixType_,OrderingType_> > Eigen::SparseLU< MatrixType_, OrderingType_ >::adjoint ( )

inline

Returns: an expression of the adjoint of the factored matrix

A typical usage is to solve for the adjoint problem A' x = b:

solver.compute(A);

x = solver.adjoint().solve(b);

For real scalar types, this function is equivalent to transpose().

See also: transpose(), solve()

Definition at line 224 of file SparseLU.h.

     {
       SparseLUTransposeView<true,  SparseLU<MatrixType_,OrderingType_> > adjointView;
       adjointView.setSparseLU(this);
       adjointView.setIsInitialized(this->m_isInitialized);
       return adjointView;
     }

◆ analyzePattern()

template<typename MatrixType , typename OrderingType >

void Eigen::SparseLU< MatrixType, OrderingType >::analyzePattern ( const MatrixType & mat )

Compute the column permutation to minimize the fill-in

Apply this permutation to the input matrix -
Compute the column elimination tree on the permuted matrix
Postorder the elimination tree and the column permutation

Definition at line 513 of file SparseLU.h.

 {
   
   //TODO  It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
   
   // Firstly, copy the whole input matrix. 
   m_mat = mat;
   
   // Compute fill-in ordering
   OrderingType ord; 
   ord(m_mat,m_perm_c);
   
   // Apply the permutation to the column of the input  matrix
   if (m_perm_c.size())
   {
     m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used.  
     // Then, permute only the column pointers
     ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?const_cast<StorageIndex*>(mat.outerIndexPtr()):0);
     
     // If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is thus needed.
     if(!mat.isCompressed()) 
       IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
     
     // Apply the permutation and compute the nnz per column.
     for (Index i = 0; i < mat.cols(); i++)
     {
       m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
       m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
     }
   }
   
   // Compute the column elimination tree of the permuted matrix 
   IndexVector firstRowElt;
   internal::coletree(m_mat, m_etree,firstRowElt); 
      
   // In symmetric mode, do not do postorder here
   if (!m_symmetricmode) {
     IndexVector post, iwork; 
     // Post order etree
     internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post); 
       
    
     // Renumber etree in postorder 
     Index m = m_mat.cols(); 
     iwork.resize(m+1);
     for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
     m_etree = iwork;
     
     // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
     PermutationType post_perm(m); 
     for (Index i = 0; i < m; i++) 
       post_perm.indices()(i) = post(i); 
         
     // Combine the two permutations : postorder the permutation for future use
     if(m_perm_c.size()) {
       m_perm_c = post_perm * m_perm_c;
     }
     
   } // end postordering 
   
   m_analysisIsOk = true; 
 }

◆ cols()

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::cols ( void ) const

inline

Definition at line 233 of file SparseLU.h.

233 { return m_mat.cols(); }

◆ colsPermutation()

template<typename MatrixType_ , typename OrderingType_ >

const PermutationType& Eigen::SparseLU< MatrixType_, OrderingType_ >::colsPermutation ( ) const

inline

Returns: a reference to the column matrix permutation \( P_c^T \) such that \(P_r A P_c^T = L U\)

See also: rowsPermutation()

Definition at line 273 of file SparseLU.h.

     {
       return m_perm_c;
     }

◆ compute()

template<typename MatrixType_ , typename OrderingType_ >

void Eigen::SparseLU< MatrixType_, OrderingType_ >::compute ( const MatrixType & matrix )

inline

Compute the symbolic and numeric factorization of the input sparse matrix. The input matrix should be in column-major storage.

Definition at line 185 of file SparseLU.h.

     {
       // Analyze 
       analyzePattern(matrix); 
       //Factorize
       factorize(matrix);
     } 

◆ determinant()

template<typename MatrixType_ , typename OrderingType_ >

Scalar Eigen::SparseLU< MatrixType_, OrderingType_ >::determinant ( )

inline

Returns: The determinant of the matrix.

See also: absDeterminant(), logAbsDeterminant()

Definition at line 437 of file SparseLU.h.

     {
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       // Initialize with the determinant of the row matrix
       Scalar det = Scalar(1.);
       // Note that the diagonal blocks of U are stored in supernodes,
       // which are available in the  L part :)
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
           if(it.index() == j)
           {
             det *= it.value();
             break;
           }
         }
       }
       return (m_detPermR * m_detPermC) > 0 ? det : -det;
     }

◆ factorize()

template<typename MatrixType , typename OrderingType >

void Eigen::SparseLU< MatrixType, OrderingType >::factorize ( const MatrixType & matrix )

Numerical factorization
Interleaved with the symbolic factorization On exit, info is

= 0: successful factorization

0: if info = i, and i is

<= A->ncol: U(i,i) is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.

> A->ncol: number of bytes allocated when memory allocation failure occurred, plus A->ncol. If lwork = -1, it is the estimated amount of space needed, plus A->ncol.

Definition at line 598 of file SparseLU.h.

 {
   using internal::emptyIdxLU;
   eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 
   eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
   
   m_isInitialized = true;
   
   // Apply the column permutation computed in analyzepattern()
   //   m_mat = matrix * m_perm_c.inverse(); 
   m_mat = matrix;
   if (m_perm_c.size()) 
   {
     m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
     //Then, permute only the column pointers
     const StorageIndex * outerIndexPtr;
     if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
     else
     {
       StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols()+1];
       for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
       outerIndexPtr = outerIndexPtr_t;
     }
     for (Index i = 0; i < matrix.cols(); i++)
     {
       m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
       m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
     }
     if(!matrix.isCompressed()) delete[] outerIndexPtr;
   } 
   else 
   { //FIXME This should not be needed if the empty permutation is handled transparently
     m_perm_c.resize(matrix.cols());
     for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
   }
   
   Index m = m_mat.rows();
   Index n = m_mat.cols();
   Index nnz = m_mat.nonZeros();
   Index maxpanel = m_perfv.panel_size * m;
   // Allocate working storage common to the factor routines
   Index lwork = 0;
   Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 
   if (info) 
   {
     m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
     m_factorizationIsOk = false;
     return ; 
   }
   
   // Set up pointers for integer working arrays 
   IndexVector segrep(m); segrep.setZero();
   IndexVector parent(m); parent.setZero();
   IndexVector xplore(m); xplore.setZero();
   IndexVector repfnz(maxpanel);
   IndexVector panel_lsub(maxpanel);
   IndexVector xprune(n); xprune.setZero();
   IndexVector marker(m*internal::LUNoMarker); marker.setZero();
   
   repfnz.setConstant(-1); 
   panel_lsub.setConstant(-1);
   
   // Set up pointers for scalar working arrays 
   ScalarVector dense; 
   dense.setZero(maxpanel);
   ScalarVector tempv; 
   tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
   
   // Compute the inverse of perm_c
   PermutationType iperm_c(m_perm_c.inverse()); 
   
   // Identify initial relaxed snodes
   IndexVector relax_end(n);
   if ( m_symmetricmode == true ) 
     Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
   else
     Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
   
   
   m_perm_r.resize(m); 
   m_perm_r.indices().setConstant(-1);
   marker.setConstant(-1);
   m_detPermR = 1; // Record the determinant of the row permutation
   
   m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
   m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
   
   // Work on one 'panel' at a time. A panel is one of the following :
   //  (a) a relaxed supernode at the bottom of the etree, or
   //  (b) panel_size contiguous columns, <panel_size> defined by the user
   Index jcol; 
   Index pivrow; // Pivotal row number in the original row matrix
   Index nseg1; // Number of segments in U-column above panel row jcol
   Index nseg; // Number of segments in each U-column 
   Index irep; 
   Index i, k, jj; 
   for (jcol = 0; jcol < n; )
   {
     // Adjust panel size so that a panel won't overlap with the next relaxed snode. 
     Index panel_size = m_perfv.panel_size; // upper bound on panel width
     for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
     {
       if (relax_end(k) != emptyIdxLU) 
       {
         panel_size = k - jcol; 
         break; 
       }
     }
     if (k == n) 
       panel_size = n - jcol; 
       
     // Symbolic outer factorization on a panel of columns 
     Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); 
     
     // Numeric sup-panel updates in topological order 
     Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 
     
     // Sparse LU within the panel, and below the panel diagonal 
     for ( jj = jcol; jj< jcol + panel_size; jj++) 
     {
       k = (jj - jcol) * m; // Column index for w-wide arrays 
       
       nseg = nseg1; // begin after all the panel segments
       //Depth-first-search for the current column
       VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
       VectorBlock<IndexVector> repfnz_k(repfnz, k, m); 
       info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 
       if ( info ) 
       {
         m_lastError =  "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       // Numeric updates to this column 
       VectorBlock<ScalarVector> dense_k(dense, k, m); 
       VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1); 
       info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 
       if ( info ) 
       {
         m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Copy the U-segments to ucol(*)
       info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 
       if ( info ) 
       {
         m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Form the L-segment 
       info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
       if ( info ) 
       {
         m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR";
 #ifndef EIGEN_NO_IO
         std::ostringstream returnInfo;
         returnInfo << " ... ZERO COLUMN AT ";
         returnInfo << info;
         m_lastError += returnInfo.str();
 #endif
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Update the determinant of the row permutation matrix
       // FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot.
       if (pivrow != jj) m_detPermR = -m_detPermR;
  
       // Prune columns (0:jj-1) using column jj
       Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
       
       // Reset repfnz for this column 
       for (i = 0; i < nseg; i++)
       {
         irep = segrep(i); 
         repfnz_k(irep) = emptyIdxLU; 
       }
     } // end SparseLU within the panel  
     jcol += panel_size;  // Move to the next panel
   } // end for -- end elimination 
   
   m_detPermR = m_perm_r.determinant();
   m_detPermC = m_perm_c.determinant();
   
   // Count the number of nonzeros in factors 
   Base::countnz(n, m_nnzL, m_nnzU, m_glu); 
   // Apply permutation  to the L subscripts 
   Base::fixupL(n, m_perm_r.indices(), m_glu);
   
   // Create supernode matrix L 
   m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); 
   // Create the column major upper sparse matrix  U; 
   new (&m_Ustore) Map<SparseMatrix<Scalar, ColMajor, StorageIndex>> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() );
   
   m_info = Success;
   m_factorizationIsOk = true;
 }

◆ info()

template<typename MatrixType_ , typename OrderingType_ >

ComputationInfo Eigen::SparseLU< MatrixType_, OrderingType_ >::info ( ) const

inline

Reports whether previous computation was successful.

Returns: Success if computation was successful, NumericalIssue if the LU factorization reports a problem, zero diagonal for instance InvalidInput if the input matrix is invalid

See also: iparm()

Definition at line 302 of file SparseLU.h.

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

◆ initperfvalues()

template<typename MatrixType_ , typename OrderingType_ >

void Eigen::SparseLU< MatrixType_, OrderingType_ >::initperfvalues ( )

inlineprotected

Definition at line 463 of file SparseLU.h.

     {
       m_perfv.panel_size = 16;
       m_perfv.relax = 1; 
       m_perfv.maxsuper = 128; 
       m_perfv.rowblk = 16; 
       m_perfv.colblk = 8; 
       m_perfv.fillfactor = 20;  
     }

◆ isSymmetric()

template<typename MatrixType_ , typename OrderingType_ >

void Eigen::SparseLU< MatrixType_, OrderingType_ >::isSymmetric ( bool sym )

inline

Indicate that the pattern of the input matrix is symmetric

Definition at line 235 of file SparseLU.h.

     {
       m_symmetricmode = sym;
     }

◆ lastErrorMessage()

template<typename MatrixType_ , typename OrderingType_ >

std::string Eigen::SparseLU< MatrixType_, OrderingType_ >::lastErrorMessage ( ) const

inline

Returns: A string describing the type of error

Definition at line 311 of file SparseLU.h.

     {
       return m_lastError; 
     }

◆ logAbsDeterminant()

template<typename MatrixType_ , typename OrderingType_ >

Scalar Eigen::SparseLU< MatrixType_, OrderingType_ >::logAbsDeterminant ( ) const

inline

Returns: the natural log of the absolute value of the determinant of the matrix of which **this is the QR decomposition

Note: This method is useful to work around the risk of overflow/underflow that's inherent to the determinant computation.

See also: absDeterminant(), signDeterminant()

Definition at line 383 of file SparseLU.h.

     {
       using std::log;
       using std::abs;
  
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       Scalar det = Scalar(0.);
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
           if(it.row() < j) continue;
           if(it.row() == j)
           {
             det += log(abs(it.value()));
             break;
           }
         }
       }
       return det;
     }

◆ matrixL()

template<typename MatrixType_ , typename OrderingType_ >

SparseLUMatrixLReturnType<SCMatrix> Eigen::SparseLU< MatrixType_, OrderingType_ >::matrixL ( ) const

inline

Returns: an expression of the matrix L, internally stored as supernodes The only operation available with this expression is the triangular solve
y = b; matrixL().solveInPlace(y);

Eigen::internal::y
const Scalar & y
Definition: MathFunctions.h:705

Definition at line 246 of file SparseLU.h.

     {
       return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
     }

◆ matrixU()

template<typename MatrixType_ , typename OrderingType_ >

SparseLUMatrixUReturnType<SCMatrix,Map<SparseMatrix<Scalar,ColMajor,StorageIndex> > > Eigen::SparseLU< MatrixType_, OrderingType_ >::matrixU ( ) const

inline

Returns: an expression of the matrix U, The only operation available with this expression is the triangular solve
y = b; matrixU().solveInPlace(y);

Definition at line 256 of file SparseLU.h.

     {
       return SparseLUMatrixUReturnType<SCMatrix, Map<SparseMatrix<Scalar,ColMajor,StorageIndex> > >(m_Lstore, m_Ustore);
     }

◆ nnzL()

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::nnzL ( ) const

inline

Definition at line 458 of file SparseLU.h.

458 { return m_nnzL; }

◆ nnzU()

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::nnzU ( ) const

inline

Definition at line 459 of file SparseLU.h.

459 { return m_nnzU; }

◆ rows()

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::rows ( void ) const

inline

Definition at line 232 of file SparseLU.h.

232 { return m_mat.rows(); }

◆ rowsPermutation()

template<typename MatrixType_ , typename OrderingType_ >

const PermutationType& Eigen::SparseLU< MatrixType_, OrderingType_ >::rowsPermutation ( ) const

inline

Returns: a reference to the row matrix permutation \( P_r \) such that \(P_r A P_c^T = L U\)

See also: colsPermutation()

Definition at line 265 of file SparseLU.h.

     {
       return m_perm_r;
     }

◆ setPivotThreshold()

template<typename MatrixType_ , typename OrderingType_ >

void Eigen::SparseLU< MatrixType_, OrderingType_ >::setPivotThreshold ( const RealScalar & thresh )

inline

Set the threshold used for a diagonal entry to be an acceptable pivot.

Definition at line 278 of file SparseLU.h.

     {
       m_diagpivotthresh = thresh; 
     }

◆ signDeterminant()

template<typename MatrixType_ , typename OrderingType_ >

Scalar Eigen::SparseLU< MatrixType_, OrderingType_ >::signDeterminant ( )

inline

Returns: A number representing the sign of the determinant

See also: absDeterminant(), logAbsDeterminant()

Definition at line 409 of file SparseLU.h.

     {
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       // Initialize with the determinant of the row matrix
       Index det = 1;
       // Note that the diagonal blocks of U are stored in supernodes,
       // which are available in the  L part :)
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
           if(it.index() == j)
           {
             if(it.value()<0)
               det = -det;
             else if(it.value()==0)
               return 0;
             break;
           }
         }
       }
       return det * m_detPermR * m_detPermC;
     }

◆ simplicialfactorize()

template<typename MatrixType_ , typename OrderingType_ >

void Eigen::SparseLU< MatrixType_, OrderingType_ >::simplicialfactorize ( const MatrixType & matrix )

◆ solve()

template<typename MatrixType_ , typename OrderingType_ >

template<typename Rhs >

const Solve<SparseLU, Rhs> Eigen::SparseLU< MatrixType_, OrderingType_ >::solve ( const MatrixBase< Rhs > & B ) const

inline

Returns: the solution X of \( A X = B \) using the current decomposition of A.

Warning: the destination matrix X in X = this->solve(B) must be colmun-major.

See also: compute()

◆ transpose()

template<typename MatrixType_ , typename OrderingType_ >

const SparseLUTransposeView<false,SparseLU<MatrixType_,OrderingType_> > Eigen::SparseLU< MatrixType_, OrderingType_ >::transpose ( )

inline

Returns: an expression of the transposed of the factored matrix.

A typical usage is to solve for the transposed problem A^T x = b:

solver.compute(A);

x = solver.transpose().solve(b);

See also: adjoint(), solve()

Definition at line 203 of file SparseLU.h.

     {
       SparseLUTransposeView<false,  SparseLU<MatrixType_,OrderingType_> > transposeView;
       transposeView.setSparseLU(this);
       transposeView.setIsInitialized(this->m_isInitialized);
       return transposeView;
     }

Member Data Documentation

◆ m_analysisIsOk

template<typename MatrixType_ , typename OrderingType_ >

bool Eigen::SparseLU< MatrixType_, OrderingType_ >::m_analysisIsOk

protected

Definition at line 476 of file SparseLU.h.

◆ m_detPermC

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::m_detPermC

protected

Definition at line 493 of file SparseLU.h.

◆ m_detPermR

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::m_detPermR

protected

Definition at line 493 of file SparseLU.h.

◆ m_diagpivotthresh

template<typename MatrixType_ , typename OrderingType_ >

RealScalar Eigen::SparseLU< MatrixType_, OrderingType_ >::m_diagpivotthresh

protected

Definition at line 491 of file SparseLU.h.

◆ m_etree

template<typename MatrixType_ , typename OrderingType_ >

IndexVector Eigen::SparseLU< MatrixType_, OrderingType_ >::m_etree

protected

Definition at line 483 of file SparseLU.h.

◆ m_factorizationIsOk

template<typename MatrixType_ , typename OrderingType_ >

bool Eigen::SparseLU< MatrixType_, OrderingType_ >::m_factorizationIsOk

protected

Definition at line 475 of file SparseLU.h.

◆ m_glu

template<typename MatrixType_ , typename OrderingType_ >

Base::GlobalLU_t Eigen::SparseLU< MatrixType_, OrderingType_ >::m_glu

protected

Definition at line 485 of file SparseLU.h.

◆ m_info

template<typename MatrixType_ , typename OrderingType_ >

ComputationInfo Eigen::SparseLU< MatrixType_, OrderingType_ >::m_info

mutableprotected

Definition at line 474 of file SparseLU.h.

◆ m_lastError

template<typename MatrixType_ , typename OrderingType_ >

std::string Eigen::SparseLU< MatrixType_, OrderingType_ >::m_lastError

protected

Definition at line 477 of file SparseLU.h.

◆ m_Lstore

template<typename MatrixType_ , typename OrderingType_ >

SCMatrix Eigen::SparseLU< MatrixType_, OrderingType_ >::m_Lstore

protected

Definition at line 479 of file SparseLU.h.

◆ m_mat

template<typename MatrixType_ , typename OrderingType_ >

NCMatrix Eigen::SparseLU< MatrixType_, OrderingType_ >::m_mat

protected

Definition at line 478 of file SparseLU.h.

◆ m_nnzL

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::m_nnzL

protected

Definition at line 492 of file SparseLU.h.

◆ m_nnzU

template<typename MatrixType_ , typename OrderingType_ >

Index Eigen::SparseLU< MatrixType_, OrderingType_ >::m_nnzU

protected

Definition at line 492 of file SparseLU.h.

◆ m_perfv

template<typename MatrixType_ , typename OrderingType_ >

internal::perfvalues Eigen::SparseLU< MatrixType_, OrderingType_ >::m_perfv

protected

Definition at line 490 of file SparseLU.h.

◆ m_perm_c

template<typename MatrixType_ , typename OrderingType_ >

PermutationType Eigen::SparseLU< MatrixType_, OrderingType_ >::m_perm_c

protected

Definition at line 481 of file SparseLU.h.

◆ m_perm_r

template<typename MatrixType_ , typename OrderingType_ >

PermutationType Eigen::SparseLU< MatrixType_, OrderingType_ >::m_perm_r

protected

Definition at line 482 of file SparseLU.h.

◆ m_symmetricmode

template<typename MatrixType_ , typename OrderingType_ >

bool Eigen::SparseLU< MatrixType_, OrderingType_ >::m_symmetricmode

protected

Definition at line 488 of file SparseLU.h.

◆ m_Ustore

template<typename MatrixType_ , typename OrderingType_ >

Map<SparseMatrix<Scalar,ColMajor,StorageIndex> > Eigen::SparseLU< MatrixType_, OrderingType_ >::m_Ustore

protected

Definition at line 480 of file SparseLU.h.

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

SparseLU.h

Public Types

Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Private Member Functions

Detailed Description

template<typename MatrixType_, typename OrderingType_> class Eigen::SparseLU< MatrixType_, OrderingType_ >

Member Typedef Documentation

◆ APIBase

◆ Base

◆ IndexVector

◆ MatrixType

◆ NCMatrix

◆ OrderingType

◆ PermutationType

◆ RealScalar

◆ Scalar

◆ ScalarVector

◆ SCMatrix

◆ StorageIndex

Member Enumeration Documentation

◆ anonymous enum

Constructor & Destructor Documentation

◆ SparseLU() [1/3]

◆ SparseLU() [2/3]

◆ ~SparseLU()

◆ SparseLU() [3/3]

Member Function Documentation

◆ _solve_impl()

◆ absDeterminant()

◆ adjoint()

◆ analyzePattern()

◆ cols()

◆ colsPermutation()

◆ compute()

◆ determinant()

◆ factorize()

◆ info()

◆ initperfvalues()

◆ isSymmetric()

◆ lastErrorMessage()

◆ logAbsDeterminant()

◆ matrixL()

◆ matrixU()

◆ nnzL()

◆ nnzU()

◆ rows()

◆ rowsPermutation()

◆ setPivotThreshold()

◆ signDeterminant()

◆ simplicialfactorize()

◆ solve()

◆ transpose()

Member Data Documentation

◆ m_analysisIsOk

◆ m_detPermC

◆ m_detPermR

◆ m_diagpivotthresh

◆ m_etree

◆ m_factorizationIsOk

◆ m_glu

◆ m_info

◆ m_lastError

◆ m_Lstore

◆ m_mat

◆ m_nnzL

◆ m_nnzU

◆ m_perfv

◆ m_perm_c

◆ m_perm_r

◆ m_symmetricmode

◆ m_Ustore

template<typename MatrixType_, typename OrderingType_>
class Eigen::SparseLU< MatrixType_, OrderingType_ >