EIC code displayed by LXR master/​include/​Rivet/​Math/​eigen3/​src/​SparseLU/​SparseLU.h	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-04-19 09:06:42

 // 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 RivetEigen {
 
 template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::StorageIndex> > class SparseLU;
 template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType;
 template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType;
 
 template <bool Conjugate,class SparseLUType>
 class SparseLUTransposeView : public SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> >
 {
 protected:
   typedef SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> > APIBase;
   using APIBase::m_isInitialized;
 public:
   typedef typename SparseLUType::Scalar Scalar;
   typedef typename SparseLUType::StorageIndex StorageIndex;
   typedef typename SparseLUType::MatrixType MatrixType;
   typedef typename SparseLUType::OrderingType OrderingType;
 
   enum {
     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
   };
 
   SparseLUTransposeView() : m_sparseLU(NULL) {}
   SparseLUTransposeView(const SparseLUTransposeView& view) {
     this->m_sparseLU = view.m_sparseLU;
   }
   void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;}
   void setSparseLU(SparseLUType* sparseLU) {m_sparseLU = sparseLU;}
   using APIBase::_solve_impl;
   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_sparseLU->info() == Success && "The matrix should be factorized first");
     EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
                         THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
 
 
     // this ugly const_cast_derived() helps to detect aliasing when applying the permutations
     for(Index j = 0; j < B.cols(); ++j){
       X.col(j) = m_sparseLU->colsPermutation() * B.const_cast_derived().col(j);
     }
     //Forward substitution with transposed or adjoint of U
     m_sparseLU->matrixU().template solveTransposedInPlace<Conjugate>(X);
 
     //Backward substitution with transposed or adjoint of L
     m_sparseLU->matrixL().template solveTransposedInPlace<Conjugate>(X);
 
     // Permute back the solution
     for (Index j = 0; j < B.cols(); ++j)
       X.col(j) = m_sparseLU->rowsPermutation().transpose() * X.col(j);
     return true;
   }
   inline Index rows() const { return m_sparseLU->rows(); }
   inline Index cols() const { return m_sparseLU->cols(); }
 
 private:
   SparseLUType *m_sparseLU;
   SparseLUTransposeView& operator=(const SparseLUTransposeView&);
 };
 
 
 /** \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 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  \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> 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); 
   * \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
   *
   * \implsparsesolverconcept
   * 
   * \sa \ref TutorialSparseSolverConcept
   * \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;
 
     enum {
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
     };
     
   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);
     } 
 
     /** \returns an expression of the transposed of the factored matrix.
       *
       * A typical usage is to solve for the transposed problem A^T x = b:
       * \code
       * solver.compute(A);
       * x = solver.transpose().solve(b);
       * \endcode
       *
       * \sa adjoint(), solve()
       */
     const SparseLUTransposeView<false,SparseLU<_MatrixType,_OrderingType> > transpose()
     {
       SparseLUTransposeView<false,  SparseLU<_MatrixType,_OrderingType> > transposeView;
       transposeView.setSparseLU(this);
       transposeView.setIsInitialized(this->m_isInitialized);
       return transposeView;
     }
 
 
     /** \returns an expression of the adjoint of the factored matrix
       *
       * A typical usage is to solve for the adjoint problem A' x = b:
       * \code
       * solver.compute(A);
       * x = solver.adjoint().solve(b);
       * \endcode
       *
       * For real scalar types, this function is equivalent to transpose().
       *
       * \sa transpose(), solve()
       */
     const SparseLUTransposeView<true, SparseLU<_MatrixType,_OrderingType> > adjoint()
     {
       SparseLUTransposeView<true,  SparseLU<_MatrixType,_OrderingType> > adjointView;
       adjointView.setSparseLU(this);
       adjointView.setIsInitialized(this->m_isInitialized);
       return adjointView;
     }
     
     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 successful,
       *          \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;
     }
 
     Index nnzL() const { return m_nnzL; };
     Index nnzU() const { return m_nnzU; };
 
   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");
   
   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 ... 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() const { return m_mapL.rows(); }
   Index cols() const { return m_mapL.cols(); }
   template<typename Dest>
   void solveInPlace( MatrixBase<Dest> &X) const
   {
     m_mapL.solveInPlace(X);
   }
   template<bool Conjugate, typename Dest>
   void solveTransposedInPlace( MatrixBase<Dest> &X) const
   {
     m_mapL.template solveTransposedInPlace<Conjugate>(X);
   }
 
   const MappedSupernodalType& m_mapL;
 };
 
 template<typename MatrixLType, typename MatrixUType>
 struct SparseLUMatrixUReturnType : internal::no_assignment_operator
 {
   typedef typename MatrixLType::Scalar Scalar;
   SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU)
   : m_mapL(mapL),m_mapU(mapU)
   { }
   Index rows() const { return m_mapL.rows(); }
   Index cols() const { 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
       {
         // FIXME: the following lines should use Block expressions and not Map!
         Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
         Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X.coeffRef(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
   }
 
   template<bool Conjugate, typename Dest>   void solveTransposedInPlace(MatrixBase<Dest> &X) const
   {
     using numext::conj;
     Index nrhs = X.cols();
     Index n    = X.rows();
     // Forward solve with U
     for (Index k = 0; k <=  m_mapL.nsuper(); 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];
 
       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(jcol, j) -= X(irow, j) * (Conjugate? conj(it.value()): it.value());
           }
         }
       }
       if (nsupc == 1)
       {
         for (Index j = 0; j < nrhs; j++)
         {
           X(fsupc, j) /= (Conjugate? conj(m_mapL.valuePtr()[luptr]) : 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,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
         if(Conjugate)
           U = A.adjoint().template triangularView<Lower>().solve(U);
         else
           U = A.transpose().template triangularView<Lower>().solve(U);
       }
     }// End For U-solve
   }
 
 
   const MatrixLType& m_mapL;
   const MatrixUType& m_mapU;
 };
 
 } // End namespace RivetEigen 
 
 #endif

This page was automatically generated by the 2.3.7 LXR engine. The LXR team