| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> |
| // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| |
| #ifndef EIGEN_SPARSE_LU_H |
| #define EIGEN_SPARSE_LU_H |
| |
| namespace Eigen { |
| |
| template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::StorageIndex> > class SparseLU; |
| template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType; |
| template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType; |
| |
| /** \ingroup SparseLU_Module |
| * \class SparseLU |
| * |
| * \brief 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 arithmetics 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 \link OrderingMethods_Module the OrderingMethods module \endlink for the list of |
| * built-in and external ordering methods. |
| * |
| * Simple example with key steps |
| * \code |
| * VectorXd x(n), b(n); |
| * SparseMatrix<double, ColMajor> A; |
| * SparseLU<SparseMatrix<scalar, ColMajor>, COLAMDOrdering<Index> > 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); |
| * \endcode |
| * |
| * \warning The input matrix A should be in a \b compressed and \b 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" |
| * |
| * \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<> |
| * \tparam _OrderingType The ordering method to use, either AMD, COLAMD or METIS. Default is COLMAD |
| * |
| * |
| * \sa \ref TutorialSparseDirectSolvers |
| * \sa \ref OrderingMethods_Module |
| */ |
| template <typename _MatrixType, typename _OrderingType> |
| class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex> |
| { |
| protected: |
| typedef SparseSolverBase<SparseLU<_MatrixType,_OrderingType> > APIBase; |
| using APIBase::m_isInitialized; |
| public: |
| using APIBase::_solve_impl; |
| |
| typedef _MatrixType MatrixType; |
| typedef _OrderingType OrderingType; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef typename MatrixType::StorageIndex StorageIndex; |
| typedef SparseMatrix<Scalar,ColMajor,StorageIndex> NCMatrix; |
| typedef internal::MappedSuperNodalMatrix<Scalar, StorageIndex> SCMatrix; |
| typedef Matrix<Scalar,Dynamic,1> ScalarVector; |
| typedef Matrix<StorageIndex,Dynamic,1> IndexVector; |
| typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; |
| typedef internal::SparseLUImpl<Scalar, StorageIndex> Base; |
| |
| public: |
| SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) |
| { |
| initperfvalues(); |
| } |
| explicit SparseLU(const MatrixType& matrix):m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) |
| { |
| initperfvalues(); |
| compute(matrix); |
| } |
| |
| ~SparseLU() |
| { |
| // Free all explicit dynamic pointers |
| } |
| |
| void analyzePattern (const MatrixType& matrix); |
| void factorize (const MatrixType& matrix); |
| void simplicialfactorize(const MatrixType& matrix); |
| |
| /** |
| * Compute the symbolic and numeric factorization of the input sparse matrix. |
| * The input matrix should be in column-major storage. |
| */ |
| void compute (const MatrixType& matrix) |
| { |
| // Analyze |
| analyzePattern(matrix); |
| //Factorize |
| factorize(matrix); |
| } |
| |
| inline Index rows() const { return m_mat.rows(); } |
| inline Index cols() const { return m_mat.cols(); } |
| /** Indicate that the pattern of the input matrix is symmetric */ |
| void isSymmetric(bool sym) |
| { |
| m_symmetricmode = sym; |
| } |
| |
| /** \returns an expression of the matrix L, internally stored as supernodes |
| * The only operation available with this expression is the triangular solve |
| * \code |
| * y = b; matrixL().solveInPlace(y); |
| * \endcode |
| */ |
| SparseLUMatrixLReturnType<SCMatrix> matrixL() const |
| { |
| return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore); |
| } |
| /** \returns an expression of the matrix U, |
| * The only operation available with this expression is the triangular solve |
| * \code |
| * y = b; matrixU().solveInPlace(y); |
| * \endcode |
| */ |
| SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,StorageIndex> > matrixU() const |
| { |
| return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >(m_Lstore, m_Ustore); |
| } |
| |
| /** |
| * \returns a reference to the row matrix permutation \f$ P_r \f$ such that \f$P_r A P_c^T = L U\f$ |
| * \sa colsPermutation() |
| */ |
| inline const PermutationType& rowsPermutation() const |
| { |
| return m_perm_r; |
| } |
| /** |
| * \returns a reference to the column matrix permutation\f$ P_c^T \f$ such that \f$P_r A P_c^T = L U\f$ |
| * \sa rowsPermutation() |
| */ |
| inline const PermutationType& colsPermutation() const |
| { |
| return m_perm_c; |
| } |
| /** Set the threshold used for a diagonal entry to be an acceptable pivot. */ |
| void setPivotThreshold(const RealScalar& thresh) |
| { |
| m_diagpivotthresh = thresh; |
| } |
| |
| #ifdef EIGEN_PARSED_BY_DOXYGEN |
| /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. |
| * |
| * \warning the destination matrix X in X = this->solve(B) must be colmun-major. |
| * |
| * \sa compute() |
| */ |
| template<typename Rhs> |
| inline const Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const; |
| #endif // EIGEN_PARSED_BY_DOXYGEN |
| |
| /** \brief Reports whether previous computation was successful. |
| * |
| * \returns \c Success if computation was succesful, |
| * \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance |
| * \c InvalidInput if the input matrix is invalid |
| * |
| * \sa iparm() |
| */ |
| ComputationInfo info() const |
| { |
| eigen_assert(m_isInitialized && "Decomposition is not initialized."); |
| return m_info; |
| } |
| |
| /** |
| * \returns A string describing the type of error |
| */ |
| std::string lastErrorMessage() const |
| { |
| return m_lastError; |
| } |
| |
| template<typename Rhs, typename Dest> |
| bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const |
| { |
| 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; |
| } |
| |
| /** |
| * \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. |
| * |
| * \sa logAbsDeterminant(), signDeterminant() |
| */ |
| Scalar absDeterminant() |
| { |
| 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; |
| } |
| |
| /** \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. |
| * |
| * \sa absDeterminant(), signDeterminant() |
| */ |
| Scalar logAbsDeterminant() const |
| { |
| 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; |
| } |
| |
| /** \returns A number representing the sign of the determinant |
| * |
| * \sa absDeterminant(), logAbsDeterminant() |
| */ |
| Scalar signDeterminant() |
| { |
| 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; |
| } |
| |
| /** \returns The determinant of the matrix. |
| * |
| * \sa absDeterminant(), logAbsDeterminant() |
| */ |
| Scalar determinant() |
| { |
| 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; |
| } |
| |
| protected: |
| // Functions |
| void initperfvalues() |
| { |
| 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; |
| } |
| |
| // Variables |
| mutable ComputationInfo m_info; |
| bool m_factorizationIsOk; |
| bool m_analysisIsOk; |
| std::string m_lastError; |
| NCMatrix m_mat; // The input (permuted ) matrix |
| SCMatrix m_Lstore; // The lower triangular matrix (supernodal) |
| MappedSparseMatrix<Scalar,ColMajor,StorageIndex> m_Ustore; // The upper triangular matrix |
| PermutationType m_perm_c; // Column permutation |
| PermutationType m_perm_r ; // Row permutation |
| IndexVector m_etree; // Column elimination tree |
| |
| typename Base::GlobalLU_t m_glu; |
| |
| // SparseLU options |
| bool m_symmetricmode; |
| // values for performance |
| internal::perfvalues m_perfv; |
| RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot |
| Index m_nnzL, m_nnzU; // Nonzeros in L and U factors |
| Index m_detPermR, m_detPermC; // Determinants of the permutation matrices |
| private: |
| // Disable copy constructor |
| SparseLU (const SparseLU& ); |
| |
| }; // End class SparseLU |
| |
| |
| |
| // Functions needed by the anaysis phase |
| /** |
| * 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 |
| * |
| */ |
| template <typename MatrixType, typename OrderingType> |
| void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) |
| { |
| |
| //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; |
| } |
| |
| // Functions needed by the numerical factorization phase |
| |
| |
| /** |
| * - 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. |
| */ |
| template <typename MatrixType, typename OrderingType> |
| void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) |
| { |
| using internal::emptyIdxLU; |
| eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); |
| eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); |
| |
| typedef typename IndexVector::Scalar StorageIndex; |
| |
| 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; |
| IndexVector panel_histo(n); |
| 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 ... ZERO COLUMN AT "; |
| std::ostringstream returnInfo; |
| returnInfo << info; |
| m_lastError += returnInfo.str(); |
| 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) MappedSparseMatrix<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; |
| } |
| |
| template<typename MappedSupernodalType> |
| struct SparseLUMatrixLReturnType : internal::no_assignment_operator |
| { |
| typedef typename MappedSupernodalType::Scalar Scalar; |
| explicit SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL) |
| { } |
| Index rows() { return m_mapL.rows(); } |
| Index cols() { return m_mapL.cols(); } |
| template<typename Dest> |
| void solveInPlace( MatrixBase<Dest> &X) const |
| { |
| m_mapL.solveInPlace(X); |
| } |
| const MappedSupernodalType& m_mapL; |
| }; |
| |
| template<typename MatrixLType, typename MatrixUType> |
| struct SparseLUMatrixUReturnType : internal::no_assignment_operator |
| { |
| typedef typename MatrixLType::Scalar Scalar; |
| explicit SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) |
| : m_mapL(mapL),m_mapU(mapU) |
| { } |
| Index rows() { return m_mapL.rows(); } |
| Index cols() { return m_mapL.cols(); } |
| |
| template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const |
| { |
| Index nrhs = X.cols(); |
| Index n = X.rows(); |
| // Backward solve with U |
| for (Index k = m_mapL.nsuper(); k >= 0; k--) |
| { |
| Index fsupc = m_mapL.supToCol()[k]; |
| Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension |
| Index nsupc = m_mapL.supToCol()[k+1] - fsupc; |
| Index luptr = m_mapL.colIndexPtr()[fsupc]; |
| |
| if (nsupc == 1) |
| { |
| for (Index j = 0; j < nrhs; j++) |
| { |
| X(fsupc, j) /= m_mapL.valuePtr()[luptr]; |
| } |
| } |
| else |
| { |
| Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); |
| Map< Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); |
| U = A.template triangularView<Upper>().solve(U); |
| } |
| |
| for (Index j = 0; j < nrhs; ++j) |
| { |
| for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) |
| { |
| typename MatrixUType::InnerIterator it(m_mapU, jcol); |
| for ( ; it; ++it) |
| { |
| Index irow = it.index(); |
| X(irow, j) -= X(jcol, j) * it.value(); |
| } |
| } |
| } |
| } // End For U-solve |
| } |
| const MatrixLType& m_mapL; |
| const MatrixUType& m_mapU; |
| }; |
| |
| } // End namespace Eigen |
| |
| #endif |
