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 
21 template <bool Conjugate,class SparseLUType>
22 class SparseLUTransposeView : public SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> >
23 {
24 protected:
27 public:
28  typedef typename SparseLUType::Scalar Scalar;
29  typedef typename SparseLUType::StorageIndex StorageIndex;
32 
33  enum {
34  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
35  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
36  };
37 
40  this->m_sparseLU = view.m_sparseLU;
41  }
42  void setIsInitialized(const bool isInitialized) {this->m_isInitialized = isInitialized;}
43  void setSparseLU(SparseLUType* sparseLU) {m_sparseLU = sparseLU;}
45  template<typename Rhs, typename Dest>
46  bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const
47  {
48  Dest& X(X_base.derived());
49  eigen_assert(m_sparseLU->info() == Success && "The matrix should be factorized first");
50  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
51  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
52 
53 
54  // this ugly const_cast_derived() helps to detect aliasing when applying the permutations
55  for(Index j = 0; j < B.cols(); ++j){
56  X.col(j) = m_sparseLU->colsPermutation() * B.const_cast_derived().col(j);
57  }
58  //Forward substitution with transposed or adjoint of U
59  m_sparseLU->matrixU().template solveTransposedInPlace<Conjugate>(X);
60 
61  //Backward substitution with transposed or adjoint of L
62  m_sparseLU->matrixL().template solveTransposedInPlace<Conjugate>(X);
63 
64  // Permute back the solution
65  for (Index j = 0; j < B.cols(); ++j)
66  X.col(j) = m_sparseLU->rowsPermutation().transpose() * X.col(j);
67  return true;
68  }
69  inline Index rows() const { return m_sparseLU->rows(); }
70  inline Index cols() const { return m_sparseLU->cols(); }
71 
72 private:
73  SparseLUType *m_sparseLU;
75 };
76 
77 
130 template <typename _MatrixType, typename _OrderingType>
131 class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex>
132 {
133  protected:
136  public:
137  using APIBase::_solve_impl;
138 
139  typedef _MatrixType MatrixType;
140  typedef _OrderingType OrderingType;
141  typedef typename MatrixType::Scalar Scalar;
143  typedef typename MatrixType::StorageIndex StorageIndex;
150 
151  enum {
152  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
153  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
154  };
155 
156  public:
157 
159  {
160  initperfvalues();
161  }
162  explicit SparseLU(const MatrixType& matrix)
163  : m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
164  {
165  initperfvalues();
166  compute(matrix);
167  }
168 
170  {
171  // Free all explicit dynamic pointers
172  }
173 
174  void analyzePattern (const MatrixType& matrix);
175  void factorize (const MatrixType& matrix);
177 
182  void compute (const MatrixType& matrix)
183  {
184  // Analyze
186  //Factorize
187  factorize(matrix);
188  }
189 
201  {
203  transposeView.setSparseLU(this);
204  transposeView.setIsInitialized(this->m_isInitialized);
205  return transposeView;
206  }
207 
208 
222  {
224  adjointView.setSparseLU(this);
225  adjointView.setIsInitialized(this->m_isInitialized);
226  return adjointView;
227  }
228 
229  inline Index rows() const { return m_mat.rows(); }
230  inline Index cols() const { return m_mat.cols(); }
232  void isSymmetric(bool sym)
233  {
234  m_symmetricmode = sym;
235  }
236 
244  {
246  }
254  {
256  }
257 
262  inline const PermutationType& rowsPermutation() const
263  {
264  return m_perm_r;
265  }
270  inline const PermutationType& colsPermutation() const
271  {
272  return m_perm_c;
273  }
275  void setPivotThreshold(const RealScalar& thresh)
276  {
277  m_diagpivotthresh = thresh;
278  }
279 
280 #ifdef EIGEN_PARSED_BY_DOXYGEN
281 
287  template<typename Rhs>
288  inline const Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const;
289 #endif // EIGEN_PARSED_BY_DOXYGEN
290 
300  {
301  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
302  return m_info;
303  }
304 
308  std::string lastErrorMessage() const
309  {
310  return m_lastError;
311  }
312 
313  template<typename Rhs, typename Dest>
314  bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const
315  {
316  Dest& X(X_base.derived());
317  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
318  EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
319  THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
320 
321  // Permute the right hand side to form X = Pr*B
322  // on return, X is overwritten by the computed solution
323  X.resize(B.rows(),B.cols());
324 
325  // this ugly const_cast_derived() helps to detect aliasing when applying the permutations
326  for(Index j = 0; j < B.cols(); ++j)
327  X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
328 
329  //Forward substitution with L
330  this->matrixL().solveInPlace(X);
331  this->matrixU().solveInPlace(X);
332 
333  // Permute back the solution
334  for (Index j = 0; j < B.cols(); ++j)
335  X.col(j) = colsPermutation().inverse() * X.col(j);
336 
337  return true;
338  }
339 
351  {
352  using std::abs;
353  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
354  // Initialize with the determinant of the row matrix
355  Scalar det = Scalar(1.);
356  // Note that the diagonal blocks of U are stored in supernodes,
357  // which are available in the L part :)
358  for (Index j = 0; j < this->cols(); ++j)
359  {
360  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
361  {
362  if(it.index() == j)
363  {
364  det *= abs(it.value());
365  break;
366  }
367  }
368  }
369  return det;
370  }
371 
381  {
382  using std::log;
383  using std::abs;
384 
385  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
386  Scalar det = Scalar(0.);
387  for (Index j = 0; j < this->cols(); ++j)
388  {
389  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
390  {
391  if(it.row() < j) continue;
392  if(it.row() == j)
393  {
394  det += log(abs(it.value()));
395  break;
396  }
397  }
398  }
399  return det;
400  }
401 
407  {
408  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
409  // Initialize with the determinant of the row matrix
410  Index det = 1;
411  // Note that the diagonal blocks of U are stored in supernodes,
412  // which are available in the L part :)
413  for (Index j = 0; j < this->cols(); ++j)
414  {
415  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
416  {
417  if(it.index() == j)
418  {
419  if(it.value()<0)
420  det = -det;
421  else if(it.value()==0)
422  return 0;
423  break;
424  }
425  }
426  }
427  return det * m_detPermR * m_detPermC;
428  }
429 
435  {
436  eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
437  // Initialize with the determinant of the row matrix
438  Scalar det = Scalar(1.);
439  // Note that the diagonal blocks of U are stored in supernodes,
440  // which are available in the L part :)
441  for (Index j = 0; j < this->cols(); ++j)
442  {
443  for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
444  {
445  if(it.index() == j)
446  {
447  det *= it.value();
448  break;
449  }
450  }
451  }
452  return (m_detPermR * m_detPermC) > 0 ? det : -det;
453  }
454 
455  Index nnzL() const { return m_nnzL; };
456  Index nnzU() const { return m_nnzU; };
457 
458  protected:
459  // Functions
461  {
462  m_perfv.panel_size = 16;
463  m_perfv.relax = 1;
464  m_perfv.maxsuper = 128;
465  m_perfv.rowblk = 16;
466  m_perfv.colblk = 8;
467  m_perfv.fillfactor = 20;
468  }
469 
470  // Variables
474  std::string m_lastError;
475  NCMatrix m_mat; // The input (permuted ) matrix
476  SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
478  PermutationType m_perm_c; // Column permutation
479  PermutationType m_perm_r ; // Row permutation
480  IndexVector m_etree; // Column elimination tree
481 
483 
484  // SparseLU options
486  // values for performance
488  RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
489  Index m_nnzL, m_nnzU; // Nonzeros in L and U factors
490  Index m_detPermR, m_detPermC; // Determinants of the permutation matrices
491  private:
492  // Disable copy constructor
493  SparseLU (const SparseLU& );
494 }; // End class SparseLU
495 
496 
497 
498 // Functions needed by the anaysis phase
509 template <typename MatrixType, typename OrderingType>
511 {
512 
513  //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
514 
515  // Firstly, copy the whole input matrix.
516  m_mat = mat;
517 
518  // Compute fill-in ordering
519  OrderingType ord;
520  ord(m_mat,m_perm_c);
521 
522  // Apply the permutation to the column of the input matrix
523  if (m_perm_c.size())
524  {
525  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.
526  // Then, permute only the column pointers
527  ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?const_cast<StorageIndex*>(mat.outerIndexPtr()):0);
528 
529  // 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.
530  if(!mat.isCompressed())
531  IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
532 
533  // Apply the permutation and compute the nnz per column.
534  for (Index i = 0; i < mat.cols(); i++)
535  {
536  m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
537  m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
538  }
539  }
540 
541  // Compute the column elimination tree of the permuted matrix
542  IndexVector firstRowElt;
543  internal::coletree(m_mat, m_etree,firstRowElt);
544 
545  // In symmetric mode, do not do postorder here
546  if (!m_symmetricmode) {
547  IndexVector post, iwork;
548  // Post order etree
549  internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post);
550 
551 
552  // Renumber etree in postorder
553  Index m = m_mat.cols();
554  iwork.resize(m+1);
555  for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
556  m_etree = iwork;
557 
558  // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
559  PermutationType post_perm(m);
560  for (Index i = 0; i < m; i++)
561  post_perm.indices()(i) = post(i);
562 
563  // Combine the two permutations : postorder the permutation for future use
564  if(m_perm_c.size()) {
565  m_perm_c = post_perm * m_perm_c;
566  }
567 
568  } // end postordering
569 
570  m_analysisIsOk = true;
571 }
572 
573 // Functions needed by the numerical factorization phase
574 
575 
594 template <typename MatrixType, typename OrderingType>
596 {
597  using internal::emptyIdxLU;
598  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
599  eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
600 
601  m_isInitialized = true;
602 
603  // Apply the column permutation computed in analyzepattern()
604  // m_mat = matrix * m_perm_c.inverse();
605  m_mat = matrix;
606  if (m_perm_c.size())
607  {
608  m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
609  //Then, permute only the column pointers
610  const StorageIndex * outerIndexPtr;
611  if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
612  else
613  {
614  StorageIndex* outerIndexPtr_t = new StorageIndex[matrix.cols()+1];
615  for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
616  outerIndexPtr = outerIndexPtr_t;
617  }
618  for (Index i = 0; i < matrix.cols(); i++)
619  {
620  m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
621  m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
622  }
623  if(!matrix.isCompressed()) delete[] outerIndexPtr;
624  }
625  else
626  { //FIXME This should not be needed if the empty permutation is handled transparently
627  m_perm_c.resize(matrix.cols());
628  for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
629  }
630 
631  Index m = m_mat.rows();
632  Index n = m_mat.cols();
633  Index nnz = m_mat.nonZeros();
634  Index maxpanel = m_perfv.panel_size * m;
635  // Allocate working storage common to the factor routines
636  Index lwork = 0;
637  Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
638  if (info)
639  {
640  m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
641  m_factorizationIsOk = false;
642  return ;
643  }
644 
645  // Set up pointers for integer working arrays
646  IndexVector segrep(m); segrep.setZero();
647  IndexVector parent(m); parent.setZero();
648  IndexVector xplore(m); xplore.setZero();
649  IndexVector repfnz(maxpanel);
650  IndexVector panel_lsub(maxpanel);
651  IndexVector xprune(n); xprune.setZero();
652  IndexVector marker(m*internal::LUNoMarker); marker.setZero();
653 
654  repfnz.setConstant(-1);
655  panel_lsub.setConstant(-1);
656 
657  // Set up pointers for scalar working arrays
658  ScalarVector dense;
659  dense.setZero(maxpanel);
660  ScalarVector tempv;
661  tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
662 
663  // Compute the inverse of perm_c
664  PermutationType iperm_c(m_perm_c.inverse());
665 
666  // Identify initial relaxed snodes
667  IndexVector relax_end(n);
668  if ( m_symmetricmode == true )
669  Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
670  else
671  Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
672 
673 
674  m_perm_r.resize(m);
675  m_perm_r.indices().setConstant(-1);
676  marker.setConstant(-1);
677  m_detPermR = 1; // Record the determinant of the row permutation
678 
679  m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
680  m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
681 
682  // Work on one 'panel' at a time. A panel is one of the following :
683  // (a) a relaxed supernode at the bottom of the etree, or
684  // (b) panel_size contiguous columns, <panel_size> defined by the user
685  Index jcol;
686  Index pivrow; // Pivotal row number in the original row matrix
687  Index nseg1; // Number of segments in U-column above panel row jcol
688  Index nseg; // Number of segments in each U-column
689  Index irep;
690  Index i, k, jj;
691  for (jcol = 0; jcol < n; )
692  {
693  // Adjust panel size so that a panel won't overlap with the next relaxed snode.
694  Index panel_size = m_perfv.panel_size; // upper bound on panel width
695  for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
696  {
697  if (relax_end(k) != emptyIdxLU)
698  {
699  panel_size = k - jcol;
700  break;
701  }
702  }
703  if (k == n)
704  panel_size = n - jcol;
705 
706  // Symbolic outer factorization on a panel of columns
707  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);
708 
709  // Numeric sup-panel updates in topological order
710  Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
711 
712  // Sparse LU within the panel, and below the panel diagonal
713  for ( jj = jcol; jj< jcol + panel_size; jj++)
714  {
715  k = (jj - jcol) * m; // Column index for w-wide arrays
716 
717  nseg = nseg1; // begin after all the panel segments
718  //Depth-first-search for the current column
719  VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
720  VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
721  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);
722  if ( info )
723  {
724  m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
725  m_info = NumericalIssue;
726  m_factorizationIsOk = false;
727  return;
728  }
729  // Numeric updates to this column
730  VectorBlock<ScalarVector> dense_k(dense, k, m);
731  VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1);
732  info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
733  if ( info )
734  {
735  m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
736  m_info = NumericalIssue;
737  m_factorizationIsOk = false;
738  return;
739  }
740 
741  // Copy the U-segments to ucol(*)
742  info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
743  if ( info )
744  {
745  m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
746  m_info = NumericalIssue;
747  m_factorizationIsOk = false;
748  return;
749  }
750 
751  // Form the L-segment
752  info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
753  if ( info )
754  {
755  m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
756  std::ostringstream returnInfo;
757  returnInfo << info;
758  m_lastError += returnInfo.str();
759  m_info = NumericalIssue;
760  m_factorizationIsOk = false;
761  return;
762  }
763 
764  // Update the determinant of the row permutation matrix
765  // FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot.
766  if (pivrow != jj) m_detPermR = -m_detPermR;
767 
768  // Prune columns (0:jj-1) using column jj
769  Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
770 
771  // Reset repfnz for this column
772  for (i = 0; i < nseg; i++)
773  {
774  irep = segrep(i);
775  repfnz_k(irep) = emptyIdxLU;
776  }
777  } // end SparseLU within the panel
778  jcol += panel_size; // Move to the next panel
779  } // end for -- end elimination
780 
781  m_detPermR = m_perm_r.determinant();
782  m_detPermC = m_perm_c.determinant();
783 
784  // Count the number of nonzeros in factors
785  Base::countnz(n, m_nnzL, m_nnzU, m_glu);
786  // Apply permutation to the L subscripts
787  Base::fixupL(n, m_perm_r.indices(), m_glu);
788 
789  // Create supernode matrix L
790  m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
791  // Create the column major upper sparse matrix U;
792  new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, StorageIndex> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() );
793 
794  m_info = Success;
795  m_factorizationIsOk = true;
796 }
797 
798 template<typename MappedSupernodalType>
800 {
802  explicit SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL)
803  { }
804  Index rows() const { return m_mapL.rows(); }
805  Index cols() const { return m_mapL.cols(); }
806  template<typename Dest>
808  {
809  m_mapL.solveInPlace(X);
810  }
811  template<bool Conjugate, typename Dest>
813  {
814  m_mapL.template solveTransposedInPlace<Conjugate>(X);
815  }
816 
817  const MappedSupernodalType& m_mapL;
818 };
819 
820 template<typename MatrixLType, typename MatrixUType>
822 {
823  typedef typename MatrixLType::Scalar Scalar;
824  SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU)
825  : m_mapL(mapL),m_mapU(mapU)
826  { }
827  Index rows() const { return m_mapL.rows(); }
828  Index cols() const { return m_mapL.cols(); }
829 
830  template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const
831  {
832  Index nrhs = X.cols();
833  Index n = X.rows();
834  // Backward solve with U
835  for (Index k = m_mapL.nsuper(); k >= 0; k--)
836  {
837  Index fsupc = m_mapL.supToCol()[k];
838  Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension
839  Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
840  Index luptr = m_mapL.colIndexPtr()[fsupc];
841 
842  if (nsupc == 1)
843  {
844  for (Index j = 0; j < nrhs; j++)
845  {
846  X(fsupc, j) /= m_mapL.valuePtr()[luptr];
847  }
848  }
849  else
850  {
851  // FIXME: the following lines should use Block expressions and not Map!
852  Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
854  U = A.template triangularView<Upper>().solve(U);
855  }
856 
857  for (Index j = 0; j < nrhs; ++j)
858  {
859  for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
860  {
861  typename MatrixUType::InnerIterator it(m_mapU, jcol);
862  for ( ; it; ++it)
863  {
864  Index irow = it.index();
865  X(irow, j) -= X(jcol, j) * it.value();
866  }
867  }
868  }
869  } // End For U-solve
870  }
871 
872  template<bool Conjugate, typename Dest> void solveTransposedInPlace(MatrixBase<Dest> &X) const
873  {
874  using numext::conj;
875  Index nrhs = X.cols();
876  Index n = X.rows();
877  // Forward solve with U
878  for (Index k = 0; k <= m_mapL.nsuper(); k++)
879  {
880  Index fsupc = m_mapL.supToCol()[k];
881  Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension
882  Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
883  Index luptr = m_mapL.colIndexPtr()[fsupc];
884 
885  for (Index j = 0; j < nrhs; ++j)
886  {
887  for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
888  {
889  typename MatrixUType::InnerIterator it(m_mapU, jcol);
890  for ( ; it; ++it)
891  {
892  Index irow = it.index();
893  X(jcol, j) -= X(irow, j) * (Conjugate? conj(it.value()): it.value());
894  }
895  }
896  }
897  if (nsupc == 1)
898  {
899  for (Index j = 0; j < nrhs; j++)
900  {
901  X(fsupc, j) /= (Conjugate? conj(m_mapL.valuePtr()[luptr]) : m_mapL.valuePtr()[luptr]);
902  }
903  }
904  else
905  {
906  Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
908  if(Conjugate)
909  U = A.adjoint().template triangularView<Lower>().solve(U);
910  else
911  U = A.transpose().template triangularView<Lower>().solve(U);
912  }
913  }// End For U-solve
914  }
915 
916 
917  const MatrixLType& m_mapL;
918  const MatrixUType& m_mapU;
919 };
920 
921 } // End namespace Eigen
922 
923 #endif
Eigen::SparseMatrix::cols
Index cols() const
Definition: SparseMatrix.h:140
Eigen::SparseLUMatrixLReturnType::Scalar
MappedSupernodalType::Scalar Scalar
Definition: SparseLU.h:801
Eigen::conj
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:574
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:444
Eigen::internal::perfvalues::fillfactor
Index fillfactor
Definition: SparseLU_Structs.h:104
Eigen::internal::treePostorder
void treePostorder(typename IndexVector::Scalar n, IndexVector &parent, IndexVector &post)
Post order a tree.
Definition: SparseColEtree.h:178
Eigen::SparseLUMatrixUReturnType::Scalar
MatrixLType::Scalar Scalar
Definition: SparseLU.h:823
Eigen::SparseLUMatrixUReturnType::m_mapL
const MatrixLType & m_mapL
Definition: SparseLU.h:917
Eigen::SparseLU::m_perm_c
PermutationType m_perm_c
Definition: SparseLU.h:478
Eigen::SparseLU::logAbsDeterminant
Scalar logAbsDeterminant() const
Definition: SparseLU.h:380
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::SparseMatrix< Scalar, ColMajor, StorageIndex >
Eigen::SparseLUMatrixLReturnType::cols
Index cols() const
Definition: SparseLU.h:805
Eigen::SparseLU::m_nnzL
Index m_nnzL
Definition: SparseLU.h:489
Eigen::SparseLU::m_Ustore
MappedSparseMatrix< Scalar, ColMajor, StorageIndex > m_Ustore
Definition: SparseLU.h:477
Eigen::SparseLUTransposeView::cols
Index cols() const
Definition: SparseLU.h:70
Eigen::SparseLUTransposeView::rows
Index rows() const
Definition: SparseLU.h:69
Eigen::SparseLU::m_perm_r
PermutationType m_perm_r
Definition: SparseLU.h:479
Eigen::SparseLU::transpose
const SparseLUTransposeView< false, SparseLU< _MatrixType, _OrderingType > > transpose()
Definition: SparseLU.h:200
Eigen::SparseLU::PermutationType
PermutationMatrix< Dynamic, Dynamic, StorageIndex > PermutationType
Definition: SparseLU.h:148
Eigen::SparseLU::lastErrorMessage
std::string lastErrorMessage() const
Definition: SparseLU.h:308
Eigen::SparseSolverBase< SparseLUTransposeView< Conjugate, SparseLUType > >::_solve_impl
void _solve_impl(const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
Definition: SparseSolverBase.h:111
Eigen::SparseLUMatrixUReturnType::SparseLUMatrixUReturnType
SparseLUMatrixUReturnType(const MatrixLType &mapL, const MatrixUType &mapU)
Definition: SparseLU.h:824
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::PermutationMatrix::indices
const IndicesType & indices() const
Definition: PermutationMatrix.h:360
Eigen::SparseLU::nnzL
Index nnzL() const
Definition: SparseLU.h:455
Eigen::SparseLU::SCMatrix
internal::MappedSuperNodalMatrix< Scalar, StorageIndex > SCMatrix
Definition: SparseLU.h:145
Eigen::SparseLU::m_glu
Base::GlobalLU_t m_glu
Definition: SparseLU.h:482
Eigen::SparseLU::m_lastError
std::string m_lastError
Definition: SparseLU.h:474
Eigen::SparseLU::m_analysisIsOk
bool m_analysisIsOk
Definition: SparseLU.h:473
Eigen::SparseLUMatrixUReturnType::solveInPlace
void solveInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:830
Eigen::SparseLUMatrixLReturnType::m_mapL
const MappedSupernodalType & m_mapL
Definition: SparseLU.h:817
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::SparseLU::adjoint
const SparseLUTransposeView< true, SparseLU< _MatrixType, _OrderingType > > adjoint()
Definition: SparseLU.h:221
Eigen::internal::perfvalues::relax
Index relax
Definition: SparseLU_Structs.h:98
Eigen::SparseLU::m_symmetricmode
bool m_symmetricmode
Definition: SparseLU.h:485
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:66
Eigen::SparseLU::m_diagpivotthresh
RealScalar m_diagpivotthresh
Definition: SparseLU.h:488
B
Definition: test_numpy_dtypes.cpp:299
Eigen::SparseLU::StorageIndex
MatrixType::StorageIndex StorageIndex
Definition: SparseLU.h:143
Eigen::SparseLU::m_Lstore
SCMatrix m_Lstore
Definition: SparseLU.h:476
Eigen::Success
@ Success
Definition: Constants.h:442
Eigen::SparseLU::m_nnzU
Index m_nnzU
Definition: SparseLU.h:489
Eigen::SparseLUMatrixUReturnType::rows
Index rows() const
Definition: SparseLU.h:827
Eigen::SparseLU::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:299
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
Eigen::SparseLU::compute
void compute(const MatrixType &matrix)
Definition: SparseLU.h:182
X
#define X
Definition: icosphere.cpp:20
Eigen::SparseLU::SparseLU
SparseLU()
Definition: SparseLU.h:158
Eigen::internal::perfvalues::maxsuper
Index maxsuper
Definition: SparseLU_Structs.h:101
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
Eigen::SparseLUTransposeView
Definition: SparseLU.h:22
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Eigen::SparseLU::setPivotThreshold
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:275
ei_declare_aligned_stack_constructed_variable
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
Definition: Memory.h:768
Eigen::SparseLU::SparseLU
SparseLU(const MatrixType &matrix)
Definition: SparseLU.h:162
Eigen::SparseLU::simplicialfactorize
void simplicialfactorize(const MatrixType &matrix)
Eigen::internal::perfvalues::rowblk
Index rowblk
Definition: SparseLU_Structs.h:102
Eigen::SparseLU::rows
Index rows() const
Definition: SparseLU.h:229
Eigen::SparseLU::IndexVector
Matrix< StorageIndex, Dynamic, 1 > IndexVector
Definition: SparseLU.h:147
Eigen::SparseLU::m_factorizationIsOk
bool m_factorizationIsOk
Definition: SparseLU.h:472
Eigen::SparseLUTransposeView::m_sparseLU
SparseLUType * m_sparseLU
Definition: SparseLU.h:73
Eigen::SparseLU::m_mat
NCMatrix m_mat
Definition: SparseLU.h:475
Eigen::SparseLU::APIBase
SparseSolverBase< SparseLU< _MatrixType, _OrderingType > > APIBase
Definition: SparseLU.h:134
Eigen::SparseLUMatrixUReturnType::cols
Index cols() const
Definition: SparseLU.h:828
n
int n
Definition: BiCGSTAB_simple.cpp:1
view
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set view
Definition: gnuplot_common_settings.hh:27
Eigen::PlainObjectBase::setConstant
EIGEN_DEVICE_FUNC Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
Eigen::SparseLUMatrixLReturnType::solveInPlace
void solveInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:807
Eigen::SparseLU::colsPermutation
const PermutationType & colsPermutation() const
Definition: SparseLU.h:270
Eigen::SparseLU::matrixU
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, StorageIndex > > matrixU() const
Definition: SparseLU.h:253
Eigen::SparseLUTransposeView::_solve_impl
bool _solve_impl(const MatrixBase< Rhs > &B, MatrixBase< Dest > &X_base) const
Definition: SparseLU.h:46
A
Definition: test_numpy_dtypes.cpp:298
Eigen::SparseLU::~SparseLU
~SparseLU()
Definition: SparseLU.h:169
Eigen::SparseLUTransposeView::SparseLUTransposeView
SparseLUTransposeView()
Definition: SparseLU.h:38
Eigen::SparseLU::m_etree
IndexVector m_etree
Definition: SparseLU.h:480
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::SparseLU::cols
Index cols() const
Definition: SparseLU.h:230
Eigen::SparseLUMatrixUReturnType
Definition: SparseLU.h:19
Eigen::SparseLU::_solve_impl
bool _solve_impl(const MatrixBase< Rhs > &B, MatrixBase< Dest > &X_base) const
Definition: SparseLU.h:314
Eigen::SparseLU::determinant
Scalar determinant()
Definition: SparseLU.h:434
Eigen::SparseLU::Scalar
MatrixType::Scalar Scalar
Definition: SparseLU.h:141
Eigen::SparseLUTransposeView::SparseLUTransposeView
SparseLUTransposeView(const SparseLUTransposeView &view)
Definition: SparseLU.h:39
Eigen::SparseLU::ColsAtCompileTime
@ ColsAtCompileTime
Definition: SparseLU.h:152
Eigen::SparseLU::RealScalar
MatrixType::RealScalar RealScalar
Definition: SparseLU.h:142
Eigen::SparseLU::m_perfv
internal::perfvalues m_perfv
Definition: SparseLU.h:487
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
Eigen::OuterStride
Convenience specialization of Stride to specify only an outer stride See class Map for some examples.
Definition: Stride.h:106
Eigen::internal::LUnumTempV
Index LUnumTempV(Index &m, Index &w, Index &t, Index &b)
Definition: SparseLU_Memory.h:39
Eigen::SparseLU::m_detPermR
Index m_detPermR
Definition: SparseLU.h:490
Eigen::PermutationBase::inverse
InverseReturnType inverse() const
Definition: PermutationMatrix.h:185
Eigen::SparseLU::OrderingType
_OrderingType OrderingType
Definition: SparseLU.h:140
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
Eigen::SparseLUTransposeView::setSparseLU
void setSparseLU(SparseLUType *sparseLU)
Definition: SparseLU.h:43
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
lda
* lda
Definition: eigenvalues.cpp:59
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::SparseLU::NCMatrix
SparseMatrix< Scalar, ColMajor, StorageIndex > NCMatrix
Definition: SparseLU.h:144
Eigen::internal::coletree
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::StorageIndex *perm=0)
Definition: SparseColEtree.h:61
Eigen::SparseLUTransposeView::MatrixType
SparseLUType::MatrixType MatrixType
Definition: SparseLU.h:30
Eigen::internal::perfvalues::panel_size
Index panel_size
Definition: SparseLU_Structs.h:97
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::SparseLU::nnzU
Index nnzU() const
Definition: SparseLU.h:456
Eigen::Conjugate
Definition: ForwardDeclarations.h:87
Eigen::SparseLU::MatrixType
_MatrixType MatrixType
Definition: SparseLU.h:139
Eigen::SparseSolverBase
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Eigen::SparseLUMatrixLReturnType::solveTransposedInPlace
void solveTransposedInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:812
Eigen::Solve
Pseudo expression representing a solving operation.
Definition: Solve.h:62
EIGEN_STATIC_ASSERT
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
Eigen::SparseLU::Base
internal::SparseLUImpl< Scalar, StorageIndex > Base
Definition: SparseLU.h:149
Eigen::internal::emptyIdxLU
@ emptyIdxLU
Definition: SparseLU_Memory.h:38
Eigen::SparseLUTransposeView::OrderingType
SparseLUType::OrderingType OrderingType
Definition: SparseLU.h:31
return
if n return
Definition: level1_cplx_impl.h:33
Eigen::PermutationMatrix< Dynamic, Dynamic, StorageIndex >
Eigen::SparseLU::absDeterminant
Scalar absDeterminant()
Definition: SparseLU.h:350
Eigen::internal::perfvalues::colblk
Index colblk
Definition: SparseLU_Structs.h:103
Eigen::SparseLUMatrixUReturnType::solveTransposedInPlace
void solveTransposedInPlace(MatrixBase< Dest > &X) const
Definition: SparseLU.h:872
Eigen::SparseLUTransposeView::setIsInitialized
void setIsInitialized(const bool isInitialized)
Definition: SparseLU.h:42
Eigen::SparseLU::m_detPermC
Index m_detPermC
Definition: SparseLU.h:490
Eigen::SparseLU::signDeterminant
Scalar signDeterminant()
Definition: SparseLU.h:406
Eigen::VectorBlock
Expression of a fixed-size or dynamic-size sub-vector.
Definition: ForwardDeclarations.h:85
Eigen::SparseLUTransposeView::Scalar
SparseLUType::Scalar Scalar
Definition: SparseLU.h:28
Eigen::SparseLUMatrixLReturnType::SparseLUMatrixLReturnType
SparseLUMatrixLReturnType(const MappedSupernodalType &mapL)
Definition: SparseLU.h:802
Eigen::SparseLU::ScalarVector
Matrix< Scalar, Dynamic, 1 > ScalarVector
Definition: SparseLU.h:146
Eigen::internal::MappedSuperNodalMatrix< Scalar, StorageIndex >
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::SparseLUTransposeView::ColsAtCompileTime
@ ColsAtCompileTime
Definition: SparseLU.h:34
U
@ U
Definition: testDecisionTree.cpp:349
Eigen::internal::SparseLUImpl
Definition: SparseLUImpl.h:20
test_DiscreteFactorGraph.OrderingType
OrderingType
Definition: test_DiscreteFactorGraph.py:20
Eigen::SparseLUTransposeView::StorageIndex
SparseLUType::StorageIndex StorageIndex
Definition: SparseLU.h:29
Eigen::SparseLUMatrixUReturnType::m_mapU
const MatrixUType & m_mapU
Definition: SparseLU.h:918
Eigen::internal::no_assignment_operator
Definition: XprHelper.h:109
Eigen::Matrix< Scalar, Dynamic, 1 >
abs
#define abs(x)
Definition: datatypes.h:17
Eigen::SparseLU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:510
Eigen::SparseLU::rowsPermutation
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:262
Eigen::SparseLUTransposeView::operator=
SparseLUTransposeView & operator=(const SparseLUTransposeView &)
Eigen::SparseLU
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
Eigen::SparseLU::factorize
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:595
Eigen::PlainObjectBase::setZero
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
NULL
#define NULL
Definition: ccolamd.c:609
Eigen::internal::LUNoMarker
@ LUNoMarker
Definition: SparseLU_Memory.h:37
Eigen::SparseLU::isSymmetric
void isSymmetric(bool sym)
Definition: SparseLU.h:232
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:440
Eigen::SparseLU::matrixL
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:243
Eigen::SparseLU::initperfvalues
void initperfvalues()
Definition: SparseLU.h:460
Eigen::SparseLUMatrixLReturnType::rows
Index rows() const
Definition: SparseLU.h:804
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
Eigen::SparseMatrix::rows
Index rows() const
Definition: SparseMatrix.h:138
Eigen::internal::LU_GlobalLU_t
Definition: SparseLU_Structs.h:77
Eigen::SparseLU::m_info
ComputationInfo m_info
Definition: SparseLU.h:471
Eigen::SparseLUMatrixLReturnType
Definition: SparseLU.h:18
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::SparseLUTransposeView::APIBase
SparseSolverBase< SparseLUTransposeView< Conjugate, SparseLUType > > APIBase
Definition: SparseLU.h:25
Eigen::MappedSparseMatrix< Scalar, ColMajor, StorageIndex >
Eigen::SparseLUTransposeView::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: SparseLU.h:35
Eigen::internal::perfvalues
Definition: SparseLU_Structs.h:96
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
Eigen::SparseSolverBase< SparseLUTransposeView< Conjugate, SparseLUType > >::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119
Eigen::SparseLU::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: SparseLU.h:153


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:35:44