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:
|
---|
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 | */
|
---|
181 | template<typename Scalar, typename Index>
|
---|
182 | class 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 | */
|
---|
232 | template<typename Scalar, typename Index>
|
---|
233 | template<typename Dest>
|
---|
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, ColMajor> 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, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
|
---|
271 | Map< Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 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, ColMajor>, 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
|
---|