src/SparseLU/SparseLU_SupernodalMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
0012 #define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
0013
0014 namespace Eigen {
0015 namespace internal {
0016
0017 /** \ingroup SparseLU_Module
0018  * \brief a class to manipulate the L supernodal factor from the SparseLU factorization
0019  *
0020  * This class  contain the data to easily store
0021  * and manipulate the supernodes during the factorization and solution phase of Sparse LU.
0022  * Only the lower triangular matrix has supernodes.
0023  *
0024  * NOTE : This class corresponds to the SCformat structure in SuperLU
0025  *
0026  */
0027 /* TODO
0028  * InnerIterator as for sparsematrix
0029  * SuperInnerIterator to iterate through all supernodes
0030  * Function for triangular solve
0031  */
0032 template <typename _Scalar, typename _StorageIndex>
0033 class MappedSuperNodalMatrix
0034 {
0035   public:
0036     typedef _Scalar Scalar;
0037     typedef _StorageIndex StorageIndex;
0038     typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
0039     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
0040   public:
0041     MappedSuperNodalMatrix()
0042     {
0043
0044     }
0045     MappedSuperNodalMatrix(Index m, Index n,  ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
0046              IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
0047     {
0048       setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col);
0049     }
0050
0051     ~MappedSuperNodalMatrix()
0052     {
0053
0054     }
0055     /**
0056      * Set appropriate pointers for the lower triangular supernodal matrix
0057      * These infos are available at the end of the numerical factorization
0058      * FIXME This class will be modified such that it can be use in the course
0059      * of the factorization.
0060      */
0061     void setInfos(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
0062              IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
0063     {
0064       m_row = m;
0065       m_col = n;
0066       m_nzval = nzval.data();
0067       m_nzval_colptr = nzval_colptr.data();
0068       m_rowind = rowind.data();
0069       m_rowind_colptr = rowind_colptr.data();
0070       m_nsuper = col_to_sup(n);
0071       m_col_to_sup = col_to_sup.data();
0072       m_sup_to_col = sup_to_col.data();
0073     }
0074
0075     /**
0076      * Number of rows
0077      */
0078     Index rows() const { return m_row; }
0079
0080     /**
0081      * Number of columns
0082      */
0083     Index cols() const { return m_col; }
0084
0085     /**
0086      * Return the array of nonzero values packed by column
0087      *
0088      * The size is nnz
0089      */
0090     Scalar* valuePtr() {  return m_nzval; }
0091
0092     const Scalar* valuePtr() const
0093     {
0094       return m_nzval;
0095     }
0096     /**
0097      * Return the pointers to the beginning of each column in \ref valuePtr()
0098      */
0099     StorageIndex* colIndexPtr()
0100     {
0101       return m_nzval_colptr;
0102     }
0103
0104     const StorageIndex* colIndexPtr() const
0105     {
0106       return m_nzval_colptr;
0107     }
0108
0109     /**
0110      * Return the array of compressed row indices of all supernodes
0111      */
0112     StorageIndex* rowIndex()  { return m_rowind; }
0113
0114     const StorageIndex* rowIndex() const
0115     {
0116       return m_rowind;
0117     }
0118
0119     /**
0120      * Return the location in \em rowvaluePtr() which starts each column
0121      */
0122     StorageIndex* rowIndexPtr() { return m_rowind_colptr; }
0123
0124     const StorageIndex* rowIndexPtr() const
0125     {
0126       return m_rowind_colptr;
0127     }
0128
0129     /**
0130      * Return the array of column-to-supernode mapping
0131      */
0132     StorageIndex* colToSup()  { return m_col_to_sup; }
0133
0134     const StorageIndex* colToSup() const
0135     {
0136       return m_col_to_sup;
0137     }
0138     /**
0139      * Return the array of supernode-to-column mapping
0140      */
0141     StorageIndex* supToCol() { return m_sup_to_col; }
0142
0143     const StorageIndex* supToCol() const
0144     {
0145       return m_sup_to_col;
0146     }
0147
0148     /**
0149      * Return the number of supernodes
0150      */
0151     Index nsuper() const
0152     {
0153       return m_nsuper;
0154     }
0155
0156     class InnerIterator;
0157     template<typename Dest>
0158     void solveInPlace( MatrixBase<Dest>&X) const;
0159     template<bool Conjugate, typename Dest>
0160     void solveTransposedInPlace( MatrixBase<Dest>&X) const;
0161
0162
0163
0164
0165
0166   protected:
0167     Index m_row; // Number of rows
0168     Index m_col; // Number of columns
0169     Index m_nsuper; // Number of supernodes
0170     Scalar* m_nzval; //array of nonzero values packed by column
0171     StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
0172     StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes
0173     StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
0174     StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
0175     StorageIndex* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
0176
0177   private :
0178 };
0179
0180 /**
0181   * \brief InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L
0182   *
0183   */
0184 template<typename Scalar, typename StorageIndex>
0185 class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator
0186 {
0187   public:
0188      InnerIterator(const MappedSuperNodalMatrix& mat, Index outer)
0189       : m_matrix(mat),
0190         m_outer(outer),
0191         m_supno(mat.colToSup()[outer]),
0192         m_idval(mat.colIndexPtr()[outer]),
0193         m_startidval(m_idval),
0194         m_endidval(mat.colIndexPtr()[outer+1]),
0195         m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
0196         m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
0197     {}
0198     inline InnerIterator& operator++()
0199     {
0200       m_idval++;
0201       m_idrow++;
0202       return *this;
0203     }
0204     inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
0205
0206     inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
0207
0208     inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
0209     inline Index row() const { return index(); }
0210     inline Index col() const { return m_outer; }
0211
0212     inline Index supIndex() const { return m_supno; }
0213
0214     inline operator bool() const
0215     {
0216       return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
0217                 && (m_idrow < m_endidrow) );
0218     }
0219
0220   protected:
0221     const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix
0222     const Index m_outer;                    // Current column
0223     const Index m_supno;                    // Current SuperNode number
0224     Index m_idval;                          // Index to browse the values in the current column
0225     const Index m_startidval;               // Start of the column value
0226     const Index m_endidval;                 // End of the column value
0227     Index m_idrow;                          // Index to browse the row indices
0228     Index m_endidrow;                       // End index of row indices of the current column
0229 };
0230
0231 /**
0232  * \brief Solve with the supernode triangular matrix
0233  *
0234  */
0235 template<typename Scalar, typename Index_>
0236 template<typename Dest>
0237 void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) const
0238 {
0239     /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */
0240 //    eigen_assert(X.rows() <= NumTraits<Index>::highest());
0241 //    eigen_assert(X.cols() <= NumTraits<Index>::highest());
0242     Index n    = int(X.rows());
0243     Index nrhs = Index(X.cols());
0244     const Scalar * Lval = valuePtr();                 // Nonzero values
0245     Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);     // working vector
0246     work.setZero();
0247     for (Index k = 0; k <= nsuper(); k ++)
0248     {
0249       Index fsupc = supToCol()[k];                    // First column of the current supernode
0250       Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
0251       Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
0252       Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
0253       Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
0254       Index irow;                                     //Current index row
0255
0256       if (nsupc == 1 )
0257       {
0258         for (Index j = 0; j < nrhs; j++)
0259         {
0260           InnerIterator it(*this, fsupc);
0261           ++it; // Skip the diagonal element
0262           for (; it; ++it)
0263           {
0264             irow = it.row();
0265             X(irow, j) -= X(fsupc, j) * it.value();
0266           }
0267         }
0268       }
0269       else
0270       {
0271         // The supernode has more than one column
0272         Index luptr = colIndexPtr()[fsupc];
0273         Index lda = colIndexPtr()[fsupc+1] - luptr;
0274
0275         // Triangular solve
0276         Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
0277         Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
0278         U = A.template triangularView<UnitLower>().solve(U);
0279
0280         // Matrix-vector product
0281         new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
0282         work.topRows(nrow).noalias() = A * U;
0283
0284         //Begin Scatter
0285         for (Index j = 0; j < nrhs; j++)
0286         {
0287           Index iptr = istart + nsupc;
0288           for (Index i = 0; i < nrow; i++)
0289           {
0290             irow = rowIndex()[iptr];
0291             X(irow, j) -= work(i, j); // Scatter operation
0292             work(i, j) = Scalar(0);
0293             iptr++;
0294           }
0295         }
0296       }
0297     }
0298 }
0299
0300 template<typename Scalar, typename Index_>
0301 template<bool Conjugate, typename Dest>
0302 void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<Dest>&X) const
0303 {
0304     using numext::conj;
0305   Index n    = int(X.rows());
0306   Index nrhs = Index(X.cols());
0307   const Scalar * Lval = valuePtr();                 // Nonzero values
0308   Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);     // working vector
0309   work.setZero();
0310   for (Index k = nsuper(); k >= 0; k--)
0311   {
0312     Index fsupc = supToCol()[k];                    // First column of the current supernode
0313     Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
0314     Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
0315     Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
0316     Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
0317     Index irow;                                     //Current index row
0318
0319     if (nsupc == 1 )
0320     {
0321       for (Index j = 0; j < nrhs; j++)
0322       {
0323         InnerIterator it(*this, fsupc);
0324         ++it; // Skip the diagonal element
0325         for (; it; ++it)
0326         {
0327           irow = it.row();
0328           X(fsupc,j) -= X(irow, j) * (Conjugate?conj(it.value()):it.value());
0329         }
0330       }
0331     }
0332     else
0333     {
0334       // The supernode has more than one column
0335       Index luptr = colIndexPtr()[fsupc];
0336       Index lda = colIndexPtr()[fsupc+1] - luptr;
0337
0338       //Begin Gather
0339       for (Index j = 0; j < nrhs; j++)
0340       {
0341         Index iptr = istart + nsupc;
0342         for (Index i = 0; i < nrow; i++)
0343         {
0344           irow = rowIndex()[iptr];
0345           work.topRows(nrow)(i,j)= X(irow,j); // Gather operation
0346           iptr++;
0347         }
0348       }
0349
0350       // Matrix-vector product with transposed submatrix
0351       Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
0352       Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
0353       if(Conjugate)
0354         U = U - A.adjoint() * work.topRows(nrow);
0355       else
0356         U = U - A.transpose() * work.topRows(nrow);
0357
0358       // Triangular solve (of transposed diagonal block)
0359       new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
0360       if(Conjugate)
0361         U = A.adjoint().template triangularView<UnitUpper>().solve(U);
0362       else
0363         U = A.transpose().template triangularView<UnitUpper>().solve(U);
0364
0365     }
0366
0367   }
0368 }
0369
0370
0371 } // end namespace internal
0372
0373 } // end namespace Eigen
0374
0375 #endif // EIGEN_SPARSELU_MATRIX_H