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  }
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() 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  }
100  {
101  return m_nzval_colptr;
102  }
103 
104  const StorageIndex* colIndexPtr() const
105  {
106  return m_nzval_colptr;
107  }
108 
113 
114  const StorageIndex* rowIndex() const
115  {
116  return m_rowind;
117  }
118 
123 
124  const StorageIndex* rowIndexPtr() const
125  {
126  return m_rowind_colptr;
127  }
128 
133 
134  const StorageIndex* colToSup() const
135  {
136  return m_col_to_sup;
137  }
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>
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 
184 template<typename Scalar, typename StorageIndex>
186 {
187  public:
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  {}
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 
235 template<typename Scalar, typename Index_>
236 template<typename Dest>
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
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 
300 template<typename Scalar, typename Index_>
301 template<bool Conjugate, typename Dest>
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) );
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
gtsam.examples.DogLegOptimizerExample.int
int
Definition: DogLegOptimizerExample.py:111
Eigen::internal::MappedSuperNodalMatrix::m_nsuper
Index m_nsuper
Definition: SparseLU_SupernodalMatrix.h:169
Eigen::conj
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:574
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_idrow
Index m_idrow
Definition: SparseLU_SupernodalMatrix.h:227
Eigen::internal::MappedSuperNodalMatrix::nsuper
Index nsuper() const
Definition: SparseLU_SupernodalMatrix.h:151
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_supno
const Index m_supno
Definition: SparseLU_SupernodalMatrix.h:223
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::InnerIterator::row
EIGEN_STRONG_INLINE Index row() const
Definition: CoreIterators.h:59
Eigen::internal::MappedSuperNodalMatrix::~MappedSuperNodalMatrix
~MappedSuperNodalMatrix()
Definition: SparseLU_SupernodalMatrix.h:51
Eigen::internal::MappedSuperNodalMatrix::ScalarVector
Matrix< Scalar, Dynamic, 1 > ScalarVector
Definition: SparseLU_SupernodalMatrix.h:39
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_endidrow
Index m_endidrow
Definition: SparseLU_SupernodalMatrix.h:228
Eigen::internal::MappedSuperNodalMatrix::rows
Index rows() const
Definition: SparseLU_SupernodalMatrix.h:78
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_outer
const Index m_outer
Definition: SparseLU_SupernodalMatrix.h:222
Eigen::internal::MappedSuperNodalMatrix::colToSup
StorageIndex * colToSup()
Definition: SparseLU_SupernodalMatrix.h:132
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::row
Index row() const
Definition: SparseLU_SupernodalMatrix.h:209
X
#define X
Definition: icosphere.cpp:20
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::internal::MappedSuperNodalMatrix::m_nzval_colptr
StorageIndex * m_nzval_colptr
Definition: SparseLU_SupernodalMatrix.h:171
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::InnerIterator
InnerIterator(const MappedSuperNodalMatrix &mat, Index outer)
Definition: SparseLU_SupernodalMatrix.h:188
Eigen::internal::MappedSuperNodalMatrix::m_sup_to_col
StorageIndex * m_sup_to_col
Definition: SparseLU_SupernodalMatrix.h:175
Eigen::internal::MappedSuperNodalMatrix::solveInPlace
void solveInPlace(MatrixBase< Dest > &X) const
Solve with the supernode triangular matrix.
Definition: SparseLU_SupernodalMatrix.h:237
Eigen::internal::MappedSuperNodalMatrix::m_rowind_colptr
StorageIndex * m_rowind_colptr
Definition: SparseLU_SupernodalMatrix.h:173
n
int n
Definition: BiCGSTAB_simple.cpp:1
Eigen::internal::MappedSuperNodalMatrix::colIndexPtr
StorageIndex * colIndexPtr()
Definition: SparseLU_SupernodalMatrix.h:99
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_idval
Index m_idval
Definition: SparseLU_SupernodalMatrix.h:224
Eigen::internal::MappedSuperNodalMatrix::m_rowind
StorageIndex * m_rowind
Definition: SparseLU_SupernodalMatrix.h:172
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::value
Scalar value() const
Definition: SparseLU_SupernodalMatrix.h:204
A
Definition: test_numpy_dtypes.cpp:298
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_matrix
const MappedSuperNodalMatrix & m_matrix
Definition: SparseLU_SupernodalMatrix.h:221
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::internal::MappedSuperNodalMatrix::m_col
Index m_col
Definition: SparseLU_SupernodalMatrix.h:168
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::valueRef
Scalar & valueRef()
Definition: SparseLU_SupernodalMatrix.h:206
Eigen::internal::MappedSuperNodalMatrix::rowIndex
StorageIndex * rowIndex()
Definition: SparseLU_SupernodalMatrix.h:112
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_startidval
const Index m_startidval
Definition: SparseLU_SupernodalMatrix.h:225
Eigen::InnerIterator::value
EIGEN_STRONG_INLINE Scalar value() const
Definition: CoreIterators.h:46
Eigen::internal::MappedSuperNodalMatrix::IndexVector
Matrix< StorageIndex, Dynamic, 1 > IndexVector
Definition: SparseLU_SupernodalMatrix.h:38
Eigen::OuterStride
Convenience specialization of Stride to specify only an outer stride See class Map for some examples.
Definition: Stride.h:106
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::internal::MappedSuperNodalMatrix::rowIndex
const StorageIndex * rowIndex() const
Definition: SparseLU_SupernodalMatrix.h:114
Eigen::internal::MappedSuperNodalMatrix::MappedSuperNodalMatrix
MappedSuperNodalMatrix(Index m, Index n, ScalarVector &nzval, IndexVector &nzval_colptr, IndexVector &rowind, IndexVector &rowind_colptr, IndexVector &col_to_sup, IndexVector &sup_to_col)
Definition: SparseLU_SupernodalMatrix.h:45
conj
AnnoyingScalar conj(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:104
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::internal::MappedSuperNodalMatrix::m_nzval
Scalar * m_nzval
Definition: SparseLU_SupernodalMatrix.h:170
Eigen::InnerIterator
An InnerIterator allows to loop over the element of any matrix expression.
Definition: CoreIterators.h:33
lda
* lda
Definition: eigenvalues.cpp:59
Eigen::internal::MappedSuperNodalMatrix::solveTransposedInPlace
void solveTransposedInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU_SupernodalMatrix.h:302
Eigen::internal::MappedSuperNodalMatrix::rowIndexPtr
const StorageIndex * rowIndexPtr() const
Definition: SparseLU_SupernodalMatrix.h:124
Eigen::internal::MappedSuperNodalMatrix::setInfos
void setInfos(Index m, Index n, ScalarVector &nzval, IndexVector &nzval_colptr, IndexVector &rowind, IndexVector &rowind_colptr, IndexVector &col_to_sup, IndexVector &sup_to_col)
Definition: SparseLU_SupernodalMatrix.h:61
Eigen::Conjugate
Definition: ForwardDeclarations.h:87
Eigen::internal::MappedSuperNodalMatrix::valuePtr
Scalar * valuePtr()
Definition: SparseLU_SupernodalMatrix.h:90
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::operator++
InnerIterator & operator++()
Definition: SparseLU_SupernodalMatrix.h:198
Eigen::internal::MappedSuperNodalMatrix::valuePtr
const Scalar * valuePtr() const
Definition: SparseLU_SupernodalMatrix.h:92
Eigen::internal::MappedSuperNodalMatrix::supToCol
const StorageIndex * supToCol() const
Definition: SparseLU_SupernodalMatrix.h:143
Eigen::internal::MappedSuperNodalMatrix::colToSup
const StorageIndex * colToSup() const
Definition: SparseLU_SupernodalMatrix.h:134
Eigen::internal::MappedSuperNodalMatrix::MappedSuperNodalMatrix
MappedSuperNodalMatrix()
Definition: SparseLU_SupernodalMatrix.h:41
Eigen::internal::MappedSuperNodalMatrix
a class to manipulate the L supernodal factor from the SparseLU factorization
Definition: SparseLU_SupernodalMatrix.h:33
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar * data() const
Definition: PlainObjectBase.h:247
U
@ U
Definition: testDecisionTree.cpp:342
Eigen::internal::MappedSuperNodalMatrix::colIndexPtr
const StorageIndex * colIndexPtr() const
Definition: SparseLU_SupernodalMatrix.h:104
Eigen::Matrix< StorageIndex, Dynamic, 1 >
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::index
Index index() const
Definition: SparseLU_SupernodalMatrix.h:208
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::supIndex
Index supIndex() const
Definition: SparseLU_SupernodalMatrix.h:212
Eigen::internal::MappedSuperNodalMatrix::rowIndexPtr
StorageIndex * rowIndexPtr()
Definition: SparseLU_SupernodalMatrix.h:122
Eigen::internal::MappedSuperNodalMatrix::m_col_to_sup
StorageIndex * m_col_to_sup
Definition: SparseLU_SupernodalMatrix.h:174
internal
Definition: BandTriangularSolver.h:13
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setZero
EIGEN_DEVICE_FUNC Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setZero(Index size)
Definition: CwiseNullaryOp.h:562
Eigen::internal::MappedSuperNodalMatrix::InnerIterator
InnerIterator class to iterate over nonzero values of the current column in the supernodal matrix L.
Definition: SparseLU_SupernodalMatrix.h:185
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::m_endidval
const Index m_endidval
Definition: SparseLU_SupernodalMatrix.h:226
Eigen::internal::MappedSuperNodalMatrix::supToCol
StorageIndex * supToCol()
Definition: SparseLU_SupernodalMatrix.h:141
Eigen::internal::MappedSuperNodalMatrix::InnerIterator::col
Index col() const
Definition: SparseLU_SupernodalMatrix.h:210
Eigen::internal::MappedSuperNodalMatrix::cols
Index cols() const
Definition: SparseLU_SupernodalMatrix.h:83
Eigen::internal::MappedSuperNodalMatrix::StorageIndex
_StorageIndex StorageIndex
Definition: SparseLU_SupernodalMatrix.h:37
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::internal::MappedSuperNodalMatrix::m_row
Index m_row
Definition: SparseLU_SupernodalMatrix.h:167
Eigen::internal::MappedSuperNodalMatrix::Scalar
_Scalar Scalar
Definition: SparseLU_SupernodalMatrix.h:36
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Sun Dec 22 2024 04:13:38