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 _Index>
33 class 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:
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 Index* colIndexPtr()
100 {
101 return m_nzval_colptr;
102 }
103
colIndexPtr()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 */
rowIndex()112 Index* rowIndex() { return m_rowind; }
113
rowIndex()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 */
rowIndexPtr()122 Index* rowIndexPtr() { return m_rowind_colptr; }
123
rowIndexPtr()124 const Index* rowIndexPtr() const
125 {
126 return m_rowind_colptr;
127 }
128
129 /**
130 * Return the array of column-to-supernode mapping
131 */
colToSup()132 Index* colToSup() { return m_col_to_sup; }
133
colToSup()134 const Index* colToSup() const
135 {
136 return m_col_to_sup;
137 }
138 /**
139 * Return the array of supernode-to-column mapping
140 */
supToCol()141 Index* supToCol() { return m_sup_to_col; }
142
supToCol()143 const Index* 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 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 */
181 template<typename Scalar, typename Index>
182 class MappedSuperNodalMatrix<Scalar,Index>::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 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