11#ifndef EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
12#define EIGEN_SPARSELU_SUPERNODAL_MATRIX_H
32template <
typename _Scalar,
typename _StorageIndex>
33class 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()
const {
return m_row; }
83 Index cols()
const {
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;
159 template<
bool Conjugate,
typename Dest>
160 void solveTransposedInPlace( MatrixBase<Dest>&X)
const;
171 StorageIndex* m_nzval_colptr;
172 StorageIndex* m_rowind;
173 StorageIndex* m_rowind_colptr;
174 StorageIndex* m_col_to_sup;
175 StorageIndex* m_sup_to_col;
184template<
typename Scalar,
typename StorageIndex>
185class MappedSuperNodalMatrix<Scalar,StorageIndex>::InnerIterator
188 InnerIterator(
const MappedSuperNodalMatrix& mat,
Index 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])
198 inline InnerIterator& operator++()
204 inline Scalar value()
const {
return m_matrix.valuePtr()[m_idval]; }
206 inline Scalar& valueRef() {
return const_cast<Scalar&
>(m_matrix.valuePtr()[m_idval]); }
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; }
212 inline Index supIndex()
const {
return m_supno; }
214 inline operator bool()
const
216 return ( (m_idval < m_endidval) && (m_idval >= m_startidval)
217 && (m_idrow < m_endidrow) );
221 const MappedSuperNodalMatrix& m_matrix;
225 const Index m_startidval;
226 const Index m_endidval;
235template<
typename Scalar,
typename Index_>
236template<
typename Dest>
237void MappedSuperNodalMatrix<Scalar,Index_>::solveInPlace( MatrixBase<Dest>&X)
const
242 Index n = int(X.rows());
244 const Scalar * Lval = valuePtr();
245 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);
247 for (
Index k = 0; k <= nsuper(); k ++)
249 Index fsupc = supToCol()[k];
250 Index istart = rowIndexPtr()[fsupc];
251 Index nsupr = rowIndexPtr()[fsupc+1] - istart;
252 Index nsupc = supToCol()[k+1] - fsupc;
253 Index nrow = nsupr - nsupc;
258 for (
Index j = 0; j < nrhs; j++)
260 InnerIterator it(*
this, fsupc);
265 X(irow, j) -= X(fsupc, j) * it.value();
272 Index luptr = colIndexPtr()[fsupc];
273 Index lda = colIndexPtr()[fsupc+1] - luptr;
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);
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;
285 for (
Index j = 0; j < nrhs; j++)
287 Index iptr = istart + nsupc;
288 for (
Index i = 0; i < nrow; i++)
290 irow = rowIndex()[iptr];
291 X(irow, j) -= work(i, j);
292 work(i, j) = Scalar(0);
300template<
typename Scalar,
typename Index_>
301template<
bool Conjugate,
typename Dest>
302void MappedSuperNodalMatrix<Scalar,Index_>::solveTransposedInPlace( MatrixBase<Dest>&X)
const
305 Index n = int(X.rows());
307 const Scalar * Lval = valuePtr();
308 Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs);
310 for (
Index k = nsuper(); k >= 0; k--)
312 Index fsupc = supToCol()[k];
313 Index istart = rowIndexPtr()[fsupc];
314 Index nsupr = rowIndexPtr()[fsupc+1] - istart;
315 Index nsupc = supToCol()[k+1] - fsupc;
316 Index nrow = nsupr - nsupc;
321 for (
Index j = 0; j < nrhs; j++)
323 InnerIterator it(*
this, fsupc);
328 X(fsupc,j) -= X(irow, j) * (Conjugate?
conj(it.value()):it.value());
335 Index luptr = colIndexPtr()[fsupc];
336 Index lda = colIndexPtr()[fsupc+1] - luptr;
339 for (
Index j = 0; j < nrhs; j++)
341 Index iptr = istart + nsupc;
342 for (
Index i = 0; i < nrow; i++)
344 irow = rowIndex()[iptr];
345 work.topRows(nrow)(i,j)= X(irow,j);
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) );
354 U = U - A.adjoint() * work.topRows(nrow);
356 U = U - A.transpose() * work.topRows(nrow);
359 new (&A) Map<
const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > ( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
361 U = A.adjoint().template triangularView<UnitUpper>().solve(U);
363 U = A.transpose().template triangularView<UnitUpper>().solve(U);
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