Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
SparseLU_SupernodalMatrix.h
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
27/* TODO
28 * InnerIterator as for sparsematrix
29 * SuperInnerIterator to iterate through all supernodes
30 * Function for triangular solve
31 */
32template <typename _Scalar, typename _StorageIndex>
33class 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:
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 }
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
78 Index rows() const { return m_row; }
79
83 Index cols() const { return m_col; }
84
90 Scalar* valuePtr() { return m_nzval; }
91
92 const Scalar* valuePtr() const
93 {
94 return m_nzval;
95 }
99 StorageIndex* colIndexPtr()
100 {
101 return m_nzval_colptr;
102 }
103
104 const StorageIndex* colIndexPtr() const
105 {
106 return m_nzval_colptr;
107 }
108
112 StorageIndex* rowIndex() { return m_rowind; }
113
114 const StorageIndex* rowIndex() const
115 {
116 return m_rowind;
117 }
118
122 StorageIndex* rowIndexPtr() { return m_rowind_colptr; }
123
124 const StorageIndex* rowIndexPtr() const
125 {
126 return m_rowind_colptr;
127 }
128
132 StorageIndex* colToSup() { return m_col_to_sup; }
133
134 const StorageIndex* colToSup() const
135 {
136 return m_col_to_sup;
137 }
141 StorageIndex* supToCol() { return m_sup_to_col; }
142
143 const StorageIndex* supToCol() const
144 {
145 return m_sup_to_col;
146 }
147
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 template<bool Conjugate, typename Dest>
160 void solveTransposedInPlace( MatrixBase<Dest>&X) const;
161
162
163
164
165
166 protected:
167 Index m_row; // Number of rows
168 Index m_col; // Number of columns
169 Index m_nsuper; // Number of supernodes
170 Scalar* m_nzval; //array of nonzero values packed by column
171 StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
172 StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes
173 StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
174 StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
175 StorageIndex* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
176
177 private :
178};
179
184template<typename Scalar, typename StorageIndex>
185class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator
186{
187 public:
188 InnerIterator(const MappedSuperNodalMatrix& mat, Index outer)
189 : m_matrix(mat),
190 m_outer(outer),
191 m_supno(mat.colToSup()[outer]),
192 m_idval(mat.colIndexPtr()[outer]),
193 m_startidval(m_idval),
194 m_endidval(mat.colIndexPtr()[outer+1]),
195 m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
196 m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
197 {}
198 inline InnerIterator& operator++()
199 {
200 m_idval++;
201 m_idrow++;
202 return *this;
203 }
204 inline Scalar value() const { return m_matrix.valuePtr()[m_idval]; }
205
206 inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix.valuePtr()[m_idval]); }
207
208 inline Index index() const { return m_matrix.rowIndex()[m_idrow]; }
209 inline Index row() const { return index(); }
210 inline Index col() const { return m_outer; }
211
212 inline Index supIndex() const { return m_supno; }
213
214 inline operator bool() const
215 {
216 return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
217 && (m_idrow < m_endidrow) );
218 }
219
220 protected:
221 const MappedSuperNodalMatrix& m_matrix; // Supernodal lower triangular matrix
222 const Index m_outer; // Current column
223 const Index m_supno; // Current SuperNode number
224 Index m_idval; // Index to browse the values in the current column
225 const Index m_startidval; // Start of the column value
226 const Index m_endidval; // End of the column value
227 Index m_idrow; // Index to browse the row indices
228 Index m_endidrow; // End index of row indices of the current column
229};
230
235template<typename Scalar, typename Index_>
236template<typename Dest>
237void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X) const
238{
239 /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */
240// eigen_assert(X.rows() <= NumTraits<Index>::highest());
241// eigen_assert(X.cols() <= NumTraits<Index>::highest());
242 Index n = int(X.rows());
243 Index nrhs = Index(X.cols());
244 const Scalar * Lval = valuePtr(); // Nonzero values
245 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector
246 work.setZero();
247 for (Index k = 0; k <= nsuper(); k ++)
248 {
249 Index fsupc = supToCol()[k]; // First column of the current supernode
250 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column
251 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode
252 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode
253 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode
254 Index irow; //Current index row
255
256 if (nsupc == 1 )
257 {
258 for (Index j = 0; j < nrhs; j++)
259 {
260 InnerIterator it(*this, fsupc);
261 ++it; // Skip the diagonal element
262 for (; it; ++it)
263 {
264 irow = it.row();
265 X(irow, j) -= X(fsupc, j) * it.value();
266 }
267 }
268 }
269 else
270 {
271 // The supernode has more than one column
272 Index luptr = colIndexPtr()[fsupc];
273 Index lda = colIndexPtr()[fsupc+1] - luptr;
274
275 // Triangular solve
276 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
277 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
278 U = A.template triangularView<UnitLower>().solve(U);
279
280 // Matrix-vector product
281 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
282 work.topRows(nrow).noalias() = A * U;
283
284 //Begin Scatter
285 for (Index j = 0; j < nrhs; j++)
286 {
287 Index iptr = istart + nsupc;
288 for (Index i = 0; i < nrow; i++)
289 {
290 irow = rowIndex()[iptr];
291 X(irow, j) -= work(i, j); // Scatter operation
292 work(i, j) = Scalar(0);
293 iptr++;
294 }
295 }
296 }
297 }
298}
299
300template<typename Scalar, typename Index_>
301template<bool Conjugate, typename Dest>
302void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<Dest>&X) const
303{
304 using numext::conj;
305 Index n = int(X.rows());
306 Index nrhs = Index(X.cols());
307 const Scalar * Lval = valuePtr(); // Nonzero values
308 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector
309 work.setZero();
310 for (Index k = nsuper(); k >= 0; k--)
311 {
312 Index fsupc = supToCol()[k]; // First column of the current supernode
313 Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column
314 Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode
315 Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode
316 Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode
317 Index irow; //Current index row
318
319 if (nsupc == 1 )
320 {
321 for (Index j = 0; j < nrhs; j++)
322 {
323 InnerIterator it(*this, fsupc);
324 ++it; // Skip the diagonal element
325 for (; it; ++it)
326 {
327 irow = it.row();
328 X(fsupc,j) -= X(irow, j) * (Conjugate?conj(it.value()):it.value());
329 }
330 }
331 }
332 else
333 {
334 // The supernode has more than one column
335 Index luptr = colIndexPtr()[fsupc];
336 Index lda = colIndexPtr()[fsupc+1] - luptr;
337
338 //Begin Gather
339 for (Index j = 0; j < nrhs; j++)
340 {
341 Index iptr = istart + nsupc;
342 for (Index i = 0; i < nrow; i++)
343 {
344 irow = rowIndex()[iptr];
345 work.topRows(nrow)(i,j)= X(irow,j); // Gather operation
346 iptr++;
347 }
348 }
349
350 // Matrix-vector product with transposed submatrix
351 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr+nsupc]), nrow, nsupc, OuterStride<>(lda) );
352 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
353 if(Conjugate)
354 U = U - A.adjoint() * work.topRows(nrow);
355 else
356 U = U - A.transpose() * work.topRows(nrow);
357
358 // Triangular solve (of transposed diagonal block)
359 new (&A) Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
360 if(Conjugate)
361 U = A.adjoint().template triangularView<UnitUpper>().solve(U);
362 else
363 U = A.transpose().template triangularView<UnitUpper>().solve(U);
364
365 }
366
367 }
368}
369
370
371} // end namespace internal
372
373} // end namespace Eigen
374
375#endif // EIGEN_SPARSELU_MATRIX_H
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74