SparseLU.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-2014 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 
12 #ifndef EIGEN_SPARSE_LU_H
13 #define EIGEN_SPARSE_LU_H
14 
15 namespace Eigen {
16 
17 template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::StorageIndex> > class SparseLU;
18 template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType;
19 template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType;
20 
73 template <typename _MatrixType, typename _OrderingType>
74 class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex>
75 {
76  protected:
79  public:
81 
82  typedef _MatrixType MatrixType;
83  typedef _OrderingType OrderingType;
84  typedef typename MatrixType::Scalar Scalar;
93 
94  enum {
95  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
96  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
97  };
98 
99  public:
101  {
102  initperfvalues();
103  }
104  explicit SparseLU(const MatrixType& matrix)
105  : m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
106  {
107  initperfvalues();
108  compute(matrix);
109  }
110 
112  {
113  // Free all explicit dynamic pointers
114  }
115 
116  void analyzePattern (const MatrixType& matrix);
117  void factorize (const MatrixType& matrix);
118  void simplicialfactorize(const MatrixType& matrix);
119 
124  void compute (const MatrixType& matrix)
125  {
126  // Analyze
127  analyzePattern(matrix);
128  //Factorize
129  factorize(matrix);
130  }
131 
132  inline Index rows() const { return m_mat.rows(); }
133  inline Index cols() const { return m_mat.cols(); }
135  void isSymmetric(bool sym)
136  {
137  m_symmetricmode = sym;
138  }
139 
147  {
149  }
157  {
159  }
160 
165  inline const PermutationType& rowsPermutation() const
166  {
167  return m_perm_r;
168  }
173  inline const PermutationType& colsPermutation() const
174  {
175  return m_perm_c;
176  }
178  void setPivotThreshold(const RealScalar& thresh)
179  {
180  m_diagpivotthresh = thresh;
181  }
182 
183 #ifdef EIGEN_PARSED_BY_DOXYGEN
184 
190  template<typename Rhs>
191  inline const Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const;
192 #endif // EIGEN_PARSED_BY_DOXYGEN
193 
203  {
204  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
205  return m_info;
206  }
207 
211  std::string lastErrorMessage() const
212  {
213  return m_lastError;
214  }
215 
216  template<typename Rhs, typename Dest>
217  bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const
218  {
219  Dest& X(X_base.derived());
220  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
221  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
222  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
223 
224  // Permute the right hand side to form X = Pr*B
225  // on return, X is overwritten by the computed solution
226  X.resize(B.rows(),B.cols());
227 
228  // this ugly const_cast_derived() helps to detect aliasing when applying the permutations
229  for(Index j = 0; j < B.cols(); ++j)
230  X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
231 
232  //Forward substitution with L
233  this->matrixL().solveInPlace(X);
234  this->matrixU().solveInPlace(X);
235 
236  // Permute back the solution
237  for (Index j = 0; j < B.cols(); ++j)
238  X.col(j) = colsPermutation().inverse() * X.col(j);
239 
240  return true;
241  }
242 
253  Scalar absDeterminant()
254  {
255  using std::abs;
256  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
257  // Initialize with the determinant of the row matrix
258  Scalar det = Scalar(1.);
259  // Note that the diagonal blocks of U are stored in supernodes,
260  // which are available in the L part :)
261  for (Index j = 0; j < this->cols(); ++j)
262  {
263  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
264  {
265  if(it.index() == j)
266  {
267  det *= abs(it.value());
268  break;
269  }
270  }
271  }
272  return det;
273  }
274 
283  Scalar logAbsDeterminant() const
284  {
285  using std::log;
286  using std::abs;
287 
288  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
289  Scalar det = Scalar(0.);
290  for (Index j = 0; j < this->cols(); ++j)
291  {
292  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
293  {
294  if(it.row() < j) continue;
295  if(it.row() == j)
296  {
297  det += log(abs(it.value()));
298  break;
299  }
300  }
301  }
302  return det;
303  }
304 
310  {
311  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
312  // Initialize with the determinant of the row matrix
313  Index det = 1;
314  // Note that the diagonal blocks of U are stored in supernodes,
315  // which are available in the L part :)
316  for (Index j = 0; j < this->cols(); ++j)
317  {
318  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
319  {
320  if(it.index() == j)
321  {
322  if(it.value()<0)
323  det = -det;
324  else if(it.value()==0)
325  return 0;
326  break;
327  }
328  }
329  }
330  return det * m_detPermR * m_detPermC;
331  }
332 
337  Scalar determinant()
338  {
339  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
340  // Initialize with the determinant of the row matrix
341  Scalar det = Scalar(1.);
342  // Note that the diagonal blocks of U are stored in supernodes,
343  // which are available in the L part :)
344  for (Index j = 0; j < this->cols(); ++j)
345  {
346  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
347  {
348  if(it.index() == j)
349  {
350  det *= it.value();
351  break;
352  }
353  }
354  }
355  return (m_detPermR * m_detPermC) > 0 ? det : -det;
356  }
357 
358  protected:
359  // Functions
361  {
362  m_perfv.panel_size = 16;
363  m_perfv.relax = 1;
364  m_perfv.maxsuper = 128;
365  m_perfv.rowblk = 16;
366  m_perfv.colblk = 8;
367  m_perfv.fillfactor = 20;
368  }
369 
370  // Variables
374  std::string m_lastError;
375  NCMatrix m_mat; // The input (permuted ) matrix
376  SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
378  PermutationType m_perm_c; // Column permutation
379  PermutationType m_perm_r ; // Row permutation
380  IndexVector m_etree; // Column elimination tree
381 
383 
384  // SparseLU options
386  // values for performance
388  RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
389  Index m_nnzL, m_nnzU; // Nonzeros in L and U factors
390  Index m_detPermR, m_detPermC; // Determinants of the permutation matrices
391  private:
392  // Disable copy constructor
393  SparseLU (const SparseLU& );
394 
395 }; // End class SparseLU
396 
397 
398 
399 // Functions needed by the anaysis phase
410 template <typename MatrixType, typename OrderingType>
412 {
413 
414  //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
415 
416  // Firstly, copy the whole input matrix.
417  m_mat = mat;
418 
419  // Compute fill-in ordering
420  OrderingType ord;
421  ord(m_mat,m_perm_c);
422 
423  // Apply the permutation to the column of the input matrix
424  if (m_perm_c.size())
425  {
426  m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used.
427  // Then, permute only the column pointers
428  ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?const_cast<StorageIndex*>(mat.outerIndexPtr()):0);
429 
430  // If the input matrix 'mat' is uncompressed, then the outer-indices do not match the ones of m_mat, and a copy is thus needed.
431  if(!mat.isCompressed())
432  IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
433 
434  // Apply the permutation and compute the nnz per column.
435  for (Index i = 0; i < mat.cols(); i++)
436  {
437  m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
438  m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
439  }
440  }
441 
442  // Compute the column elimination tree of the permuted matrix
443  IndexVector firstRowElt;
444  internal::coletree(m_mat, m_etree,firstRowElt);
445 
446  // In symmetric mode, do not do postorder here
447  if (!m_symmetricmode) {
448  IndexVector post, iwork;
449  // Post order etree
451 
452 
453  // Renumber etree in postorder
454  Index m = m_mat.cols();
455  iwork.resize(m+1);
456  for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
457  m_etree = iwork;
458 
459  // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
460  PermutationType post_perm(m);
461  for (Index i = 0; i < m; i++)
462  post_perm.indices()(i) = post(i);
463 
464  // Combine the two permutations : postorder the permutation for future use
465  if(m_perm_c.size()) {
466  m_perm_c = post_perm * m_perm_c;
467  }
468 
469  } // end postordering
470 
471  m_analysisIsOk = true;
472 }
473 
474 // Functions needed by the numerical factorization phase
475 
476 
495 template <typename MatrixType, typename OrderingType>
497 {
498  using internal::emptyIdxLU;
499  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
500  eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
501 
502  m_isInitialized = true;
503 
504 
505  // Apply the column permutation computed in analyzepattern()
506  // m_mat = matrix * m_perm_c.inverse();
507  m_mat = matrix;
508  if (m_perm_c.size())
509  {
510  m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
511  //Then, permute only the column pointers
512  const StorageIndex * outerIndexPtr;
513  if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
514  else
515  {
516  StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols()+1];
517  for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
518  outerIndexPtr = outerIndexPtr_t;
519  }
520  for (Index i = 0; i < matrix.cols(); i++)
521  {
522  m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
523  m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
524  }
525  if(!matrix.isCompressed()) delete[] outerIndexPtr;
526  }
527  else
528  { //FIXME This should not be needed if the empty permutation is handled transparently
529  m_perm_c.resize(matrix.cols());
530  for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
531  }
532 
533  Index m = m_mat.rows();
534  Index n = m_mat.cols();
535  Index nnz = m_mat.nonZeros();
536  Index maxpanel = m_perfv.panel_size * m;
537  // Allocate working storage common to the factor routines
538  Index lwork = 0;
540  if (info)
541  {
542  m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
543  m_factorizationIsOk = false;
544  return ;
545  }
546 
547  // Set up pointers for integer working arrays
548  IndexVector segrep(m); segrep.setZero();
549  IndexVector parent(m); parent.setZero();
550  IndexVector xplore(m); xplore.setZero();
551  IndexVector repfnz(maxpanel);
552  IndexVector panel_lsub(maxpanel);
553  IndexVector xprune(n); xprune.setZero();
554  IndexVector marker(m*internal::LUNoMarker); marker.setZero();
555 
556  repfnz.setConstant(-1);
557  panel_lsub.setConstant(-1);
558 
559  // Set up pointers for scalar working arrays
560  ScalarVector dense;
561  dense.setZero(maxpanel);
562  ScalarVector tempv;
563  tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
564 
565  // Compute the inverse of perm_c
566  PermutationType iperm_c(m_perm_c.inverse());
567 
568  // Identify initial relaxed snodes
569  IndexVector relax_end(n);
570  if ( m_symmetricmode == true )
571  Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
572  else
573  Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
574 
575 
576  m_perm_r.resize(m);
577  m_perm_r.indices().setConstant(-1);
578  marker.setConstant(-1);
579  m_detPermR = 1; // Record the determinant of the row permutation
580 
582  m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
583 
584  // Work on one 'panel' at a time. A panel is one of the following :
585  // (a) a relaxed supernode at the bottom of the etree, or
586  // (b) panel_size contiguous columns, <panel_size> defined by the user
587  Index jcol;
588  IndexVector panel_histo(n);
589  Index pivrow; // Pivotal row number in the original row matrix
590  Index nseg1; // Number of segments in U-column above panel row jcol
591  Index nseg; // Number of segments in each U-column
592  Index irep;
593  Index i, k, jj;
594  for (jcol = 0; jcol < n; )
595  {
596  // Adjust panel size so that a panel won't overlap with the next relaxed snode.
597  Index panel_size = m_perfv.panel_size; // upper bound on panel width
598  for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
599  {
600  if (relax_end(k) != emptyIdxLU)
601  {
602  panel_size = k - jcol;
603  break;
604  }
605  }
606  if (k == n)
607  panel_size = n - jcol;
608 
609  // Symbolic outer factorization on a panel of columns
610  Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu);
611 
612  // Numeric sup-panel updates in topological order
613  Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
614 
615  // Sparse LU within the panel, and below the panel diagonal
616  for ( jj = jcol; jj< jcol + panel_size; jj++)
617  {
618  k = (jj - jcol) * m; // Column index for w-wide arrays
619 
620  nseg = nseg1; // begin after all the panel segments
621  //Depth-first-search for the current column
622  VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
623  VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
624  info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu);
625  if ( info )
626  {
627  m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
629  m_factorizationIsOk = false;
630  return;
631  }
632  // Numeric updates to this column
633  VectorBlock<ScalarVector> dense_k(dense, k, m);
634  VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1);
635  info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
636  if ( info )
637  {
638  m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
640  m_factorizationIsOk = false;
641  return;
642  }
643 
644  // Copy the U-segments to ucol(*)
645  info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
646  if ( info )
647  {
648  m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
650  m_factorizationIsOk = false;
651  return;
652  }
653 
654  // Form the L-segment
655  info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
656  if ( info )
657  {
658  m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
659  std::ostringstream returnInfo;
660  returnInfo << info;
661  m_lastError += returnInfo.str();
663  m_factorizationIsOk = false;
664  return;
665  }
666 
667  // Update the determinant of the row permutation matrix
668  // FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot.
669  if (pivrow != jj) m_detPermR = -m_detPermR;
670 
671  // Prune columns (0:jj-1) using column jj
672  Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
673 
674  // Reset repfnz for this column
675  for (i = 0; i < nseg; i++)
676  {
677  irep = segrep(i);
678  repfnz_k(irep) = emptyIdxLU;
679  }
680  } // end SparseLU within the panel
681  jcol += panel_size; // Move to the next panel
682  } // end for -- end elimination
683 
686 
687  // Count the number of nonzeros in factors
689  // Apply permutation to the L subscripts
691 
692  // Create supernode matrix L
694  // Create the column major upper sparse matrix U;
696 
697  m_info = Success;
698  m_factorizationIsOk = true;
699 }
700 
701 template<typename MappedSupernodalType>
703 {
705  explicit SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL)
706  { }
707  Index rows() { return m_mapL.rows(); }
708  Index cols() { return m_mapL.cols(); }
709  template<typename Dest>
711  {
712  m_mapL.solveInPlace(X);
713  }
714  const MappedSupernodalType& m_mapL;
715 };
716 
717 template<typename MatrixLType, typename MatrixUType>
719 {
720  typedef typename MatrixLType::Scalar Scalar;
721  SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU)
722  : m_mapL(mapL),m_mapU(mapU)
723  { }
724  Index rows() { return m_mapL.rows(); }
725  Index cols() { return m_mapL.cols(); }
726 
727  template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const
728  {
729  Index nrhs = X.cols();
730  Index n = X.rows();
731  // Backward solve with U
732  for (Index k = m_mapL.nsuper(); k >= 0; k--)
733  {
734  Index fsupc = m_mapL.supToCol()[k];
735  Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension
736  Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
737  Index luptr = m_mapL.colIndexPtr()[fsupc];
738 
739  if (nsupc == 1)
740  {
741  for (Index j = 0; j < nrhs; j++)
742  {
743  X(fsupc, j) /= m_mapL.valuePtr()[luptr];
744  }
745  }
746  else
747  {
748  Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
750  U = A.template triangularView<Upper>().solve(U);
751  }
752 
753  for (Index j = 0; j < nrhs; ++j)
754  {
755  for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
756  {
757  typename MatrixUType::InnerIterator it(m_mapU, jcol);
758  for ( ; it; ++it)
759  {
760  Index irow = it.index();
761  X(irow, j) -= X(jcol, j) * it.value();
762  }
763  }
764  }
765  } // End For U-solve
766  }
767  const MatrixLType& m_mapL;
768  const MatrixUType& m_mapU;
769 };
770 
771 } // End namespace Eigen
772 
773 #endif
Matrix3f m
EIGEN_DEVICE_FUNC ColXpr col(Index i)
Definition: DenseBase.h:839
void _solve_impl(const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
SCMatrix m_Lstore
Definition: SparseLU.h:376
_MatrixType MatrixType
Definition: SparseLU.h:82
bool m_factorizationIsOk
Definition: SparseLU.h:372
SCALAR Scalar
Definition: bench_gemm.cpp:33
NCMatrix m_mat
Definition: SparseLU.h:375
SparseLUMatrixUReturnType(const MatrixLType &mapL, const MatrixUType &mapU)
Definition: SparseLU.h:721
Index column_bmod(const Index jcol, const Index nseg, BlockScalarVector dense, ScalarVector &tempv, BlockIndexVector segrep, BlockIndexVector repfnz, Index fpanelc, GlobalLU_t &glu)
Performs numeric block updates (sup-col) in topological order.
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
IndexVector m_etree
Definition: SparseLU.h:380
MatrixLType::Scalar Scalar
Definition: SparseLU.h:720
Index copy_to_ucol(const Index jcol, const Index nseg, IndexVector &segrep, BlockIndexVector repfnz, IndexVector &perm_r, BlockScalarVector dense, GlobalLU_t &glu)
Performs numeric block updates (sup-col) in topological order.
MatrixType::StorageIndex StorageIndex
Definition: SparseLU.h:86
const MatrixLType & m_mapL
Definition: SparseLU.h:767
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::StorageIndex *perm=0)
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:165
const StorageIndex * innerNonZeroPtr() const
Definition: SparseMatrix.h:175
#define min(a, b)
Definition: datatypes.h:19
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
ComputationInfo m_info
Definition: SparseLU.h:371
bool m_analysisIsOk
Definition: SparseLU.h:373
Matrix< Scalar, Dynamic, 1 > ScalarVector
Definition: SparseLU.h:89
A base class for sparse solvers.
int n
bool _solve_impl(const MatrixBase< Rhs > &B, MatrixBase< Dest > &X_base) const
Definition: SparseLU.h:217
void relax_snode(const Index n, IndexVector &et, const Index relax_columns, IndexVector &descendants, IndexVector &relax_end)
Identify the initial relaxed supernodes.
internal::perfvalues m_perfv
Definition: SparseLU.h:387
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
EIGEN_DEVICE_FUNC const LogReturnType log() const
void pruneL(const Index jcol, const IndexVector &perm_r, const Index pivrow, const Index nseg, const IndexVector &segrep, BlockIndexVector repfnz, IndexVector &xprune, GlobalLU_t &glu)
Prunes the L-structure.
PermutationMatrix< Dynamic, Dynamic, StorageIndex > PermutationType
Definition: SparseLU.h:91
const Solve< SparseLU< _MatrixType, _OrderingType >, Rhs > solve(const MatrixBase< Rhs > &b) const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
void heap_relax_snode(const Index n, IndexVector &et, const Index relax_columns, IndexVector &descendants, IndexVector &relax_end)
Identify the initial relaxed supernodes.
const MappedSupernodalType & m_mapL
Definition: SparseLU.h:714
void fixupL(const Index n, const IndexVector &perm_r, GlobalLU_t &glu)
Fix up the data storage lsub for L-subscripts.
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:124
void countnz(const Index n, Index &nnzL, Index &nnzU, GlobalLU_t &glu)
Count Nonzero elements in the factors.
std::string lastErrorMessage() const
Definition: SparseLU.h:211
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:496
const unsigned int RowMajorBit
Definition: Constants.h:61
void panel_bmod(const Index m, const Index w, const Index jcol, const Index nseg, ScalarVector &dense, ScalarVector &tempv, IndexVector &segrep, IndexVector &repfnz, GlobalLU_t &glu)
Performs numeric block updates (sup-panel) in topological order.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:202
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
void panel_dfs(const Index m, const Index w, const Index jcol, MatrixType &A, IndexVector &perm_r, Index &nseg, ScalarVector &dense, IndexVector &panel_lsub, IndexVector &segrep, IndexVector &repfnz, IndexVector &xprune, IndexVector &marker, IndexVector &parent, IndexVector &xplore, GlobalLU_t &glu)
Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
Index rows() const
Definition: SparseMatrix.h:136
void simplicialfactorize(const MatrixType &matrix)
void initperfvalues()
Definition: SparseLU.h:360
SparseLU(const MatrixType &matrix)
Definition: SparseLU.h:104
bool m_symmetricmode
Definition: SparseLU.h:385
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:178
Scalar absDeterminant()
Definition: SparseLU.h:253
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, StorageIndex > > matrixU() const
Definition: SparseLU.h:156
internal::traits< SparseMatrix< _Scalar, _Options, _StorageIndex > >::StorageIndex StorageIndex
Expression of a fixed-size or dynamic-size sub-vector.
Matrix< StorageIndex, Dynamic, 1 > IndexVector
Definition: SparseLU.h:90
PermutationType m_perm_c
Definition: SparseLU.h:378
void solveInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:727
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
RealScalar m_diagpivotthresh
Definition: SparseLU.h:388
Index m_detPermR
Definition: SparseLU.h:390
void compute(const MatrixType &matrix)
Definition: SparseLU.h:124
const PermutationType & colsPermutation() const
Definition: SparseLU.h:173
internal::traits< SparseMatrix< _Scalar, _Options, _StorageIndex > >::Scalar Scalar
return(x<=y?ADS(x):ADS(y))
Index cols() const
Definition: SparseMatrix.h:138
void isSymmetric(bool sym)
Definition: SparseLU.h:135
MatrixType::Scalar Scalar
Definition: SparseLU.h:84
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Index pivotL(const Index jcol, const RealScalar &diagpivotthresh, IndexVector &perm_r, IndexVector &iperm_c, Index &pivrow, GlobalLU_t &glu)
Performs the numerical pivotin on the current column of L, and the CDIV operation.
#define eigen_assert(x)
Definition: Macros.h:579
void setInfos(Index m, Index n, ScalarVector &nzval, IndexVector &nzval_colptr, IndexVector &rowind, IndexVector &rowind_colptr, IndexVector &col_to_sup, IndexVector &sup_to_col)
Scalar determinant()
Definition: SparseLU.h:337
Index LUnumTempV(Index &m, Index &w, Index &t, Index &b)
* lda
Definition: eigenvalues.cpp:59
EIGEN_DEVICE_FUNC Derived & setConstant(Index size, const Scalar &val)
Index column_dfs(const Index m, const Index jcol, IndexVector &perm_r, Index maxsuper, Index &nseg, BlockIndexVector lsub_col, IndexVector &segrep, BlockIndexVector repfnz, IndexVector &xprune, IndexVector &marker, IndexVector &parent, IndexVector &xplore, GlobalLU_t &glu)
Performs a symbolic factorization on column jcol and decide the supernode boundary.
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:644
Base::GlobalLU_t m_glu
Definition: SparseLU.h:382
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:166
MappedSparseMatrix< Scalar, ColMajor, StorageIndex > m_Ustore
Definition: SparseLU.h:377
internal::MappedSuperNodalMatrix< Scalar, StorageIndex > SCMatrix
Definition: SparseLU.h:88
InverseReturnType inverse() const
SparseMatrix< Scalar, ColMajor, StorageIndex > NCMatrix
Definition: SparseLU.h:87
_OrderingType OrderingType
Definition: SparseLU.h:83
PermutationType m_perm_r
Definition: SparseLU.h:379
internal::SparseLUImpl< Scalar, StorageIndex > Base
Definition: SparseLU.h:92
Index memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t &glu)
Allocate various working space for the numerical factorization phase.
const IndicesType & indices() const
static ConstMapType Map(const Scalar *data)
void treePostorder(typename IndexVector::Scalar n, IndexVector &parent, IndexVector &post)
Post order a tree.
Scalar logAbsDeterminant() const
Definition: SparseLU.h:283
SparseSolverBase< SparseLU< _MatrixType, _OrderingType > > APIBase
Definition: SparseLU.h:77
SparseLUMatrixLReturnType(const MappedSupernodalType &mapL)
Definition: SparseLU.h:705
Index rows() const
Definition: SparseLU.h:132
Index cols() const
Definition: SparseLU.h:133
Pseudo expression representing a solving operation.
Definition: Solve.h:62
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:411
const MatrixUType & m_mapU
Definition: SparseLU.h:768
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
MatrixType::RealScalar RealScalar
Definition: SparseLU.h:85
#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
void solveInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:710
void resize(Index newSize)
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:430
Scalar signDeterminant()
Definition: SparseLU.h:309
std::string m_lastError
Definition: SparseLU.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:146
std::ptrdiff_t j
Index m_detPermC
Definition: SparseLU.h:390
MappedSupernodalType::Scalar Scalar
Definition: SparseLU.h:704


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