Quelle SparseLU.h Sprache: C

// 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;

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 Eigen

#endif

quality88%

¤ Dauer der Verarbeitung: 0.15 Sekunden ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.