11 #ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H 12 #define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H 32 template <
typename _Scalar,
typename _StorageIndex>
33 class MappedSuperNodalMatrix
36 typedef _Scalar Scalar;
37 typedef _StorageIndex StorageIndex;
38 typedef Matrix<StorageIndex,Dynamic,1> IndexVector;
39 typedef Matrix<Scalar,Dynamic,1> ScalarVector;
41 MappedSuperNodalMatrix()
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 )
48 setInfos(m, n, nzval, nzval_colptr, rowind, rowind_colptr, col_to_sup, sup_to_col);
51 ~MappedSuperNodalMatrix()
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 )
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();
78 Index rows() {
return m_row; }
83 Index cols() {
return m_col; }
90 Scalar* valuePtr() {
return m_nzval; }
92 const Scalar* valuePtr()
const 99 StorageIndex* colIndexPtr()
101 return m_nzval_colptr;
104 const StorageIndex* colIndexPtr()
const 106 return m_nzval_colptr;
112 StorageIndex* rowIndex() {
return m_rowind; }
114 const StorageIndex* rowIndex()
const 122 StorageIndex* rowIndexPtr() {
return m_rowind_colptr; }
124 const StorageIndex* rowIndexPtr()
const 126 return m_rowind_colptr;
132 StorageIndex* colToSup() {
return m_col_to_sup; }
134 const StorageIndex* colToSup()
const 141 StorageIndex* supToCol() {
return m_sup_to_col; }
143 const StorageIndex* supToCol()
const 157 template<
typename Dest>
158 void solveInPlace( MatrixBase<Dest>&X)
const;
168 StorageIndex* m_nzval_colptr;
169 StorageIndex* m_rowind;
170 StorageIndex* m_rowind_colptr;
171 StorageIndex* m_col_to_sup;
172 StorageIndex* m_sup_to_col;
181 template<
typename Scalar,
typename StorageIndex>
182 class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator
185 InnerIterator(
const MappedSuperNodalMatrix& mat,
Index 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])
195 inline InnerIterator& operator++()
201 inline Scalar value()
const {
return m_matrix.valuePtr()[m_idval]; }
203 inline Scalar& valueRef() {
return const_cast<Scalar&
>(m_matrix.valuePtr()[m_idval]); }
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; }
209 inline Index supIndex()
const {
return m_supno; }
211 inline operator bool()
const 213 return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
214 && (m_idrow < m_endidrow) );
218 const MappedSuperNodalMatrix& m_matrix;
222 const Index m_startidval;
223 const Index m_endidval;
232 template<
typename Scalar,
typename Index_>
233 template<
typename Dest>
234 void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X)
const 239 Index n = int(X.rows());
241 const Scalar * Lval = valuePtr();
242 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);
244 for (
Index k = 0; k <= nsuper(); k ++)
246 Index fsupc = supToCol()[k];
247 Index istart = rowIndexPtr()[fsupc];
248 Index nsupr = rowIndexPtr()[fsupc+1] - istart;
249 Index nsupc = supToCol()[k+1] - fsupc;
250 Index nrow = nsupr - nsupc;
255 for (
Index j = 0; j < nrhs; j++)
257 InnerIterator it(*
this, fsupc);
262 X(irow, j) -= X(fsupc, j) * it.value();
269 Index luptr = colIndexPtr()[fsupc];
270 Index lda = colIndexPtr()[fsupc+1] - luptr;
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);
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;
282 for (
Index j = 0; j < nrhs; j++)
284 Index iptr = istart + nsupc;
285 for (
Index i = 0; i < nrow; i++)
287 irow = rowIndex()[iptr];
288 X(irow, j) -= work(i, j);
289 work(i, j) = Scalar(0);
301 #endif // EIGEN_SPARSELU_MATRIX_H Namespace containing all symbols from the Eigen library.
Definition: Core:287
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Eigen_Colamd.h:50