1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
5// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
12#define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
13
14namespace Eigen {
15namespace internal {
16
17/** \ingroup SparseLU_Module
18 * \brief a class to manipulate the L supernodal factor from the SparseLU factorization
19 *
20 * This class  contain the data to easily store
21 * and manipulate the supernodes during the factorization and solution phase of Sparse LU.
22 * Only the lower triangular matrix has supernodes.
23 *
24 * NOTE : This class corresponds to the SCformat structure in SuperLU
25 *
26 */
27/* TODO
28 * InnerIterator as for sparsematrix
29 * SuperInnerIterator to iterate through all supernodes
30 * Function for triangular solve
31 */
32template <typename _Scalar, typename _Index>
33class MappedSuperNodalMatrix
34{
35  public:
36    typedef _Scalar Scalar;
37    typedef _Index Index;
38    typedef Matrix<Index,Dynamic,1> IndexVector;
39    typedef Matrix<Scalar,Dynamic,1> ScalarVector;
40  public:
41    MappedSuperNodalMatrix()
42    {
43
44    }
45    MappedSuperNodalMatrix(Index m, Index n,  ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
46             IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
47    {
48      setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col);
49    }
50
51    ~MappedSuperNodalMatrix()
52    {
53
54    }
55    /**
56     * Set appropriate pointers for the lower triangular supernodal matrix
57     * These infos are available at the end of the numerical factorization
58     * FIXME This class will be modified such that it can be use in the course
59     * of the factorization.
60     */
61    void setInfos(Index m, Index n, ScalarVector& nzval, IndexVector& nzval_colptr, IndexVector& rowind,
62             IndexVector& rowind_colptr, IndexVector& col_to_sup, IndexVector& sup_to_col )
63    {
64      m_row = m;
65      m_col = n;
66      m_nzval = nzval.data();
67      m_nzval_colptr = nzval_colptr.data();
68      m_rowind = rowind.data();
69      m_rowind_colptr = rowind_colptr.data();
70      m_nsuper = col_to_sup(n);
71      m_col_to_sup = col_to_sup.data();
72      m_sup_to_col = sup_to_col.data();
73    }
74
75    /**
76     * Number of rows
77     */
78    Index rows() { return m_row; }
79
80    /**
81     * Number of columns
82     */
83    Index cols() { return m_col; }
84
85    /**
86     * Return the array of nonzero values packed by column
87     *
88     * The size is nnz
89     */
90    Scalar* valuePtr() {  return m_nzval; }
91
92    const Scalar* valuePtr() const
93    {
94      return m_nzval;
95    }
96    /**
97     * Return the pointers to the beginning of each column in \ref valuePtr()
98     */
99    Index* colIndexPtr()
100    {
101      return m_nzval_colptr;
102    }
103
104    const Index* colIndexPtr() const
105    {
106      return m_nzval_colptr;
107    }
108
109    /**
110     * Return the array of compressed row indices of all supernodes
111     */
112    Index* rowIndex()  { return m_rowind; }
113
114    const Index* rowIndex() const
115    {
116      return m_rowind;
117    }
118
119    /**
120     * Return the location in \em rowvaluePtr() which starts each column
121     */
122    Index* rowIndexPtr() { return m_rowind_colptr; }
123
124    const Index* rowIndexPtr() const
125    {
126      return m_rowind_colptr;
127    }
128
129    /**
130     * Return the array of column-to-supernode mapping
131     */
132    Index* colToSup()  { return m_col_to_sup; }
133
134    const Index* colToSup() const
135    {
136      return m_col_to_sup;
137    }
138    /**
139     * Return the array of supernode-to-column mapping
140     */
141    Index* supToCol() { return m_sup_to_col; }
142
143    const Index* supToCol() const
144    {
145      return m_sup_to_col;
146    }
147
148    /**
149     * Return the number of supernodes
150     */
151    Index nsuper() const
152    {
153      return m_nsuper;
154    }
155
156    class InnerIterator;
157    template<typename Dest>
158    void solveInPlace( MatrixBase<Dest>&X) const;
159
160
161
162
163  protected:
164    Index m_row; // Number of rows
165    Index m_col; // Number of columns
166    Index m_nsuper; // Number of supernodes
167    Scalar* m_nzval; //array of nonzero values packed by column
168    Index* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
169    Index* m_rowind; // Array of compressed row indices of rectangular supernodes
170    Index* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
171    Index* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
172    Index* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
173
174  private :
175};
176
177/**
178  * \brief InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L
179  *
180  */
181template<typename Scalar, typename Index>
182class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
183{
184  public:
185     InnerIterator(const MappedSuperNodalMatrix& mat, Index outer)
186      : m_matrix(mat),
187        m_outer(outer),
188        m_supno(mat.colToSup()[outer]),
189        m_idval(mat.colIndexPtr()[outer]),
190        m_startidval(m_idval),
191        m_endidval(mat.colIndexPtr()[outer+1]),
192        m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
193        m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
194    {}
195    inline InnerIterator& operator++()
196    {
197      m_idval++;
198      m_idrow++;
199      return *this;
200    }
201    inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
202
203    inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
204
205    inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
206    inline Index row() const { return index(); }
207    inline Index col() const { return m_outer; }
208
209    inline Index supIndex() const { return m_supno; }
210
211    inline operator bool() const
212    {
213      return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
214                && (m_idrow < m_endidrow) );
215    }
216
217  protected:
218    const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix
219    const Index m_outer;                    // Current column
220    const Index m_supno;                    // Current SuperNode number
221    Index m_idval;                          // Index to browse the values in the current column
222    const Index m_startidval;               // Start of the column value
223    const Index m_endidval;                 // End of the column value
224    Index m_idrow;                          // Index to browse the row indices
225    Index m_endidrow;                       // End index of row indices of the current column
226};
227
228/**
229 * \brief Solve with the supernode triangular matrix
230 *
231 */
232template<typename Scalar, typename Index>
233template<typename Dest>
234void MappedSuperNodalMatrix<Scalar,Index>::solveInPlace( MatrixBase<Dest>&X) const
235{
236    Index n = X.rows();
237    Index nrhs = X.cols();
238    const Scalar * Lval = valuePtr();                 // Nonzero values
239    Matrix<Scalar,Dynamic,Dynamic> work(n, nrhs);     // working vector
240    work.setZero();
241    for (Index k = 0; k <= nsuper(); k ++)
242    {
243      Index fsupc = supToCol()[k];                    // First column of the current supernode
244      Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
245      Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
246      Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
247      Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
248      Index irow;                                     //Current index row
249
250      if (nsupc == 1 )
251      {
252        for (Index j = 0; j < nrhs; j++)
253        {
254          InnerIterator it(*this, fsupc);
255          ++it; // Skip the diagonal element
256          for (; it; ++it)
257          {
258            irow = it.row();
259            X(irow, j) -= X(fsupc, j) * it.value();
260          }
261        }
262      }
263      else
264      {
265        // The supernode has more than one column
266        Index luptr = colIndexPtr()[fsupc];
267        Index lda = colIndexPtr()[fsupc+1] - luptr;
268
269        // Triangular solve
270        Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
271        Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
272        U = A.template triangularView<UnitLower>().solve(U);
273
274        // Matrix-vector product
275        new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
276        work.block(0, 0, nrow, nrhs) = A * U;
277
278        //Begin Scatter
279        for (Index j = 0; j < nrhs; j++)
280        {
281          Index iptr = istart + nsupc;
282          for (Index i = 0; i < nrow; i++)
283          {
284            irow = rowIndex()[iptr];
285            X(irow, j) -= work(i, j); // Scatter operation
286            work(i, j) = Scalar(0);
287            iptr++;
288          }
289        }
290      }
291    }
292}
293
294} // end namespace internal
295
296} // end namespace Eigen
297
298#endif // EIGEN_SPARSELU_MATRIX_H
299