• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
14 namespace Eigen {
15 namespace 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  */
32 template <typename _Scalar, typename _StorageIndex>
33 class MappedSuperNodalMatrix
34 {
35   public:
36     typedef _Scalar Scalar;
37     typedef _StorageIndex StorageIndex;
38     typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
39     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
40   public:
MappedSuperNodalMatrix()41     MappedSuperNodalMatrix()
42     {
43 
44     }
MappedSuperNodalMatrix(Index m,Index n,ScalarVector & nzval,IndexVector & nzval_colptr,IndexVector & rowind,IndexVector & rowind_colptr,IndexVector & col_to_sup,IndexVector & sup_to_col)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 
~MappedSuperNodalMatrix()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      */
setInfos(Index m,Index n,ScalarVector & nzval,IndexVector & nzval_colptr,IndexVector & rowind,IndexVector & rowind_colptr,IndexVector & col_to_sup,IndexVector & sup_to_col)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      */
rows()78     Index rows() { return m_row; }
79 
80     /**
81      * Number of columns
82      */
cols()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      */
valuePtr()90     Scalar* valuePtr() {  return m_nzval; }
91 
valuePtr()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      */
colIndexPtr()99     StorageIndex* colIndexPtr()
100     {
101       return m_nzval_colptr;
102     }
103 
colIndexPtr()104     const StorageIndex* colIndexPtr() const
105     {
106       return m_nzval_colptr;
107     }
108 
109     /**
110      * Return the array of compressed row indices of all supernodes
111      */
rowIndex()112     StorageIndex* rowIndex()  { return m_rowind; }
113 
rowIndex()114     const StorageIndex* rowIndex() const
115     {
116       return m_rowind;
117     }
118 
119     /**
120      * Return the location in \em rowvaluePtr() which starts each column
121      */
rowIndexPtr()122     StorageIndex* rowIndexPtr() { return m_rowind_colptr; }
123 
rowIndexPtr()124     const StorageIndex* rowIndexPtr() const
125     {
126       return m_rowind_colptr;
127     }
128 
129     /**
130      * Return the array of column-to-supernode mapping
131      */
colToSup()132     StorageIndex* colToSup()  { return m_col_to_sup; }
133 
colToSup()134     const StorageIndex* colToSup() const
135     {
136       return m_col_to_sup;
137     }
138     /**
139      * Return the array of supernode-to-column mapping
140      */
supToCol()141     StorageIndex* supToCol() { return m_sup_to_col; }
142 
supToCol()143     const StorageIndex* supToCol() const
144     {
145       return m_sup_to_col;
146     }
147 
148     /**
149      * Return the number of supernodes
150      */
nsuper()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     StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
169     StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes
170     StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
171     StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
172     StorageIndex* 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   */
181 template<typename Scalar, typename StorageIndex>
182 class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator
183 {
184   public:
InnerIterator(const MappedSuperNodalMatrix & mat,Index outer)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     }
value()201     inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
202 
valueRef()203     inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
204 
index()205     inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
row()206     inline Index row() const { return index(); }
col()207     inline Index col() const { return m_outer; }
208 
supIndex()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  */
232 template<typename Scalar, typename Index_>
233 template<typename Dest>
solveInPlace(MatrixBase<Dest> & X)234 void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) const
235 {
236     /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */
237 //    eigen_assert(X.rows() <= NumTraits<Index>::highest());
238 //    eigen_assert(X.cols() <= NumTraits<Index>::highest());
239     Index n    = int(X.rows());
240     Index nrhs = Index(X.cols());
241     const Scalar * Lval = valuePtr();                 // Nonzero values
242     Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);     // working vector
243     work.setZero();
244     for (Index k = 0; k <= nsuper(); k ++)
245     {
246       Index fsupc = supToCol()[k];                    // First column of the current supernode
247       Index istart = rowIndexPtr()[fsupc];            // Pointer index to the subscript of the current column
248       Index nsupr = rowIndexPtr()[fsupc+1] - istart;  // Number of rows in the current supernode
249       Index nsupc = supToCol()[k+1] - fsupc;          // Number of columns in the current supernode
250       Index nrow = nsupr - nsupc;                     // Number of rows in the non-diagonal part of the supernode
251       Index irow;                                     //Current index row
252 
253       if (nsupc == 1 )
254       {
255         for (Index j = 0; j < nrhs; j++)
256         {
257           InnerIterator it(*this, fsupc);
258           ++it; // Skip the diagonal element
259           for (; it; ++it)
260           {
261             irow = it.row();
262             X(irow, j) -= X(fsupc, j) * it.value();
263           }
264         }
265       }
266       else
267       {
268         // The supernode has more than one column
269         Index luptr = colIndexPtr()[fsupc];
270         Index lda = colIndexPtr()[fsupc+1] - luptr;
271 
272         // Triangular solve
273         Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
274         Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
275         U = A.template triangularView<UnitLower>().solve(U);
276 
277         // Matrix-vector product
278         new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
279         work.topRows(nrow).noalias() = A * U;
280 
281         //Begin Scatter
282         for (Index j = 0; j < nrhs; j++)
283         {
284           Index iptr = istart + nsupc;
285           for (Index i = 0; i < nrow; i++)
286           {
287             irow = rowIndex()[iptr];
288             X(irow, j) -= work(i, j); // Scatter operation
289             work(i, j) = Scalar(0);
290             iptr++;
291           }
292         }
293       }
294     }
295 }
296 
297 } // end namespace internal
298 
299 } // end namespace Eigen
300 
301 #endif // EIGEN_SPARSELU_MATRIX_H
302