SparseLU_SupernodalMatrix.h
Go to the documentation of this file.
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 
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 _StorageIndex>
34 {
35  public:
36  typedef _Scalar Scalar;
37  typedef _StorageIndex StorageIndex;
40  public:
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 
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() { return m_row; }
79 
83  Index cols() { 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 
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  StorageIndex* m_nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j
169  StorageIndex* m_rowind; // Array of compressed row indices of rectangular supernodes
170  StorageIndex* m_rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
171  StorageIndex* m_col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
172  StorageIndex* m_sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
173 
174  private :
175 };
176 
181 template<typename Scalar, typename StorageIndex>
183 {
184  public:
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  {}
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 
232 template<typename Scalar, typename Index_>
233 template<typename Dest>
235 {
236  /* Explicit type conversion as the Index type of MatrixBase<Dest> may be wider than Index */
237 // eigen_assert(X.rows() <= NumTraits<Index>::highest());
238 // eigen_assert(X.cols() <= NumTraits<Index>::highest());
239  Index n = int(X.rows());
240  Index nrhs = Index(X.cols());
241  const Scalar * Lval = valuePtr(); // Nonzero values
242  Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor> work(n, nrhs); // working vector
243  work.setZero();
244  for (Index k = 0; k <= nsuper(); k ++)
245  {
246  Index fsupc = supToCol()[k]; // First column of the current supernode
247  Index istart = rowIndexPtr()[fsupc]; // Pointer index to the subscript of the current column
248  Index nsupr = rowIndexPtr()[fsupc+1] - istart; // Number of rows in the current supernode
249  Index nsupc = supToCol()[k+1] - fsupc; // Number of columns in the current supernode
250  Index nrow = nsupr - nsupc; // Number of rows in the non-diagonal part of the supernode
251  Index irow; //Current index row
252 
253  if (nsupc == 1 )
254  {
255  for (Index j = 0; j < nrhs; j++)
256  {
257  InnerIterator it(*this, fsupc);
258  ++it; // Skip the diagonal element
259  for (; it; ++it)
260  {
261  irow = it.row();
262  X(irow, j) -= X(fsupc, j) * it.value();
263  }
264  }
265  }
266  else
267  {
268  // The supernode has more than one column
269  Index luptr = colIndexPtr()[fsupc];
270  Index lda = colIndexPtr()[fsupc+1] - luptr;
271 
272  // Triangular solve
273  Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(lda) );
275  U = A.template triangularView<UnitLower>().solve(U);
276 
277  // Matrix-vector product
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;
280 
281  //Begin Scatter
282  for (Index j = 0; j < nrhs; j++)
283  {
284  Index iptr = istart + nsupc;
285  for (Index i = 0; i < nrow; i++)
286  {
287  irow = rowIndex()[iptr];
288  X(irow, j) -= work(i, j); // Scatter operation
289  work(i, j) = Scalar(0);
290  iptr++;
291  }
292  }
293  }
294  }
295 }
296 
297 } // end namespace internal
298 
299 } // end namespace Eigen
300 
301 #endif // EIGEN_SPARSELU_MATRIX_H
Matrix3f m
EIGEN_DEVICE_FUNC Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setZero(Index size)
MappedSuperNodalMatrix(Index m, Index n, ScalarVector &nzval, IndexVector &nzval_colptr, IndexVector &rowind, IndexVector &rowind_colptr, IndexVector &col_to_sup, IndexVector &sup_to_col)
return int(ret)+1
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
void solveInPlace(MatrixBase< Dest > &X) const
Solve with the supernode triangular matrix.
int n
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L...
InnerIterator(const MappedSuperNodalMatrix &mat, Index outer)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
void setInfos(Index m, Index n, ScalarVector &nzval, IndexVector &nzval_colptr, IndexVector &rowind, IndexVector &rowind_colptr, IndexVector &col_to_sup, IndexVector &sup_to_col)
* lda
Definition: eigenvalues.cpp:59
Matrix< StorageIndex, Dynamic, 1 > IndexVector
#define X
Definition: icosphere.cpp:20
Convenience specialization of Stride to specify only an outer stride See class Map for some examples...
Definition: Stride.h:101
a class to manipulate the L supernodal factor from the SparseLU factorization
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
std::ptrdiff_t j
An InnerIterator allows to loop over the element of any matrix expression.
Definition: CoreIterators.h:33


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:25