IncompleteCholesky.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 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_INCOMPLETE_CHOlESKY_H
11 #define EIGEN_INCOMPLETE_CHOlESKY_H
13 #include <Eigen/OrderingMethods>
14 #include <list>
15 
16 namespace Eigen {
29 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
31 {
32  public:
34  typedef _OrderingType OrderingType;
36  typedef typename MatrixType::Index Index;
40  typedef std::vector<std::list<Index> > VectorList;
41  enum { UpLo = _UpLo };
42  public:
44  IncompleteCholesky(const MatrixType& matrix) : m_shift(1),m_factorizationIsOk(false)
45  {
46  compute(matrix);
47  }
48 
49  Index rows() const { return m_L.rows(); }
50 
51  Index cols() const { return m_L.cols(); }
52 
53 
60  {
61  eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
62  return m_info;
63  }
64 
68  void setShift( Scalar shift) { m_shift = shift; }
69 
73  template<typename MatrixType>
74  void analyzePattern(const MatrixType& mat)
75  {
76  OrderingType ord;
77  ord(mat.template selfadjointView<UpLo>(), m_perm);
78  m_analysisIsOk = true;
79  }
80 
81  template<typename MatrixType>
82  void factorize(const MatrixType& amat);
83 
84  template<typename MatrixType>
85  void compute (const MatrixType& matrix)
86  {
87  analyzePattern(matrix);
88  factorize(matrix);
89  }
90 
91  template<typename Rhs, typename Dest>
92  void _solve(const Rhs& b, Dest& x) const
93  {
94  eigen_assert(m_factorizationIsOk && "factorize() should be called first");
95  if (m_perm.rows() == b.rows())
96  x = m_perm.inverse() * b;
97  else
98  x = b;
99  x = m_scal.asDiagonal() * x;
100  x = m_L.template triangularView<UnitLower>().solve(x);
101  x = m_L.adjoint().template triangularView<Upper>().solve(x);
102  if (m_perm.rows() == b.rows())
103  x = m_perm * x;
104  x = m_scal.asDiagonal() * x;
105  }
106  template<typename Rhs> inline const internal::solve_retval<IncompleteCholesky, Rhs>
107  solve(const MatrixBase<Rhs>& b) const
108  {
109  eigen_assert(m_factorizationIsOk && "IncompleteLLT did not succeed");
110  eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
111  eigen_assert(cols()==b.rows()
112  && "IncompleteLLT::solve(): invalid number of rows of the right hand side matrix b");
113  return internal::solve_retval<IncompleteCholesky, Rhs>(*this, b.derived());
114  }
115  protected:
116  SparseMatrix<Scalar,ColMajor> m_L; // The lower part stored in CSC
117  ScalarType m_scal; // The vector for scaling the matrix
118  Scalar m_shift; //The initial shift parameter
123  PermutationType m_perm;
124 
125  private:
126  template <typename IdxType, typename SclType>
127  inline void updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol);
128 };
129 
130 template<typename Scalar, int _UpLo, typename OrderingType>
131 template<typename _MatrixType>
133 {
134  using std::sqrt;
135  using std::min;
136  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
137 
138  // Dropping strategies : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added
139 
140  // Apply the fill-reducing permutation computed in analyzePattern()
141  if (m_perm.rows() == mat.rows() ) // To detect the null permutation
142  m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
143  else
144  m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
145 
146  Index n = m_L.cols();
147  Index nnz = m_L.nonZeros();
148  Map<ScalarType> vals(m_L.valuePtr(), nnz); //values
149  Map<IndexType> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
150  Map<IndexType> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
151  IndexType firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
152  VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
153  ScalarType curCol(n); // Store a nonzero values in each column
154  IndexType irow(n); // Row indices of nonzero elements in each column
155 
156 
157  // Computes the scaling factors
158  m_scal.resize(n);
159  for (int j = 0; j < n; j++)
160  {
161  m_scal(j) = m_L.col(j).norm();
162  m_scal(j) = sqrt(m_scal(j));
163  }
164  // Scale and compute the shift for the matrix
165  Scalar mindiag = vals[0];
166  for (int j = 0; j < n; j++){
167  for (int k = colPtr[j]; k < colPtr[j+1]; k++)
168  vals[k] /= (m_scal(j) * m_scal(rowIdx[k]));
169  mindiag = (min)(vals[colPtr[j]], mindiag);
170  }
171 
172  if(mindiag < Scalar(0.)) m_shift = m_shift - mindiag;
173  // Apply the shift to the diagonal elements of the matrix
174  for (int j = 0; j < n; j++)
175  vals[colPtr[j]] += m_shift;
176  // jki version of the Cholesky factorization
177  for (int j=0; j < n; ++j)
178  {
179  //Left-looking factorize the column j
180  // First, load the jth column into curCol
181  Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
182  curCol.setZero();
183  irow.setLinSpaced(n,0,n-1);
184  for (int i = colPtr[j] + 1; i < colPtr[j+1]; i++)
185  {
186  curCol(rowIdx[i]) = vals[i];
187  irow(rowIdx[i]) = rowIdx[i];
188  }
189  std::list<int>::iterator k;
190  // Browse all previous columns that will update column j
191  for(k = listCol[j].begin(); k != listCol[j].end(); k++)
192  {
193  int jk = firstElt(*k); // First element to use in the column
194  jk += 1;
195  for (int i = jk; i < colPtr[*k+1]; i++)
196  {
197  curCol(rowIdx[i]) -= vals[i] * vals[jk] ;
198  }
199  updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
200  }
201 
202  // Scale the current column
203  if(RealScalar(diag) <= 0)
204  {
205  std::cerr << "\nNegative diagonal during Incomplete factorization... "<< j << "\n";
207  return;
208  }
209  RealScalar rdiag = sqrt(RealScalar(diag));
210  vals[colPtr[j]] = rdiag;
211  for (int i = j+1; i < n; i++)
212  {
213  //Scale
214  curCol(i) /= rdiag;
215  //Update the remaining diagonals with curCol
216  vals[colPtr[i]] -= curCol(i) * curCol(i);
217  }
218  // Select the largest p elements
219  // p is the original number of elements in the column (without the diagonal)
220  int p = colPtr[j+1] - colPtr[j] - 1 ;
221  internal::QuickSplit(curCol, irow, p);
222  // Insert the largest p elements in the matrix
223  int cpt = 0;
224  for (int i = colPtr[j]+1; i < colPtr[j+1]; i++)
225  {
226  vals[i] = curCol(cpt);
227  rowIdx[i] = irow(cpt);
228  cpt ++;
229  }
230  // Get the first smallest row index and put it after the diagonal element
231  Index jk = colPtr(j)+1;
232  updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
233  }
234  m_factorizationIsOk = true;
235  m_isInitialized = true;
236  m_info = Success;
237 }
238 
239 template<typename Scalar, int _UpLo, typename OrderingType>
240 template <typename IdxType, typename SclType>
241 inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(const IdxType& colPtr, IdxType& rowIdx, SclType& vals, const Index& col, const Index& jk, IndexType& firstElt, VectorList& listCol)
242 {
243  if (jk < colPtr(col+1) )
244  {
245  Index p = colPtr(col+1) - jk;
246  Index minpos;
247  rowIdx.segment(jk,p).minCoeff(&minpos);
248  minpos += jk;
249  if (rowIdx(minpos) != rowIdx(jk))
250  {
251  //Swap
252  std::swap(rowIdx(jk),rowIdx(minpos));
253  std::swap(vals(jk),vals(minpos));
254  }
255  firstElt(col) = jk;
256  listCol[rowIdx(jk)].push_back(col);
257  }
258 }
259 namespace internal {
260 
261 template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs>
262 struct solve_retval<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
263  : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
264 {
266  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
267 
268  template<typename Dest> void evalTo(Dest& dst) const
269  {
270  dec()._solve(rhs(),dst);
271  }
272 };
273 
274 } // end namespace internal
275 
276 } // end namespace Eigen
277 
278 #endif
void _solve(const Rhs &b, Dest &x) const
Index cols() const
Definition: SparseMatrix.h:121
IntermediateState sqrt(const Expression &arg)
const AdjointReturnType adjoint() const
MatrixType::RealScalar RealScalar
SparseMatrix< Scalar, ColMajor > MatrixType
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
Modified Incomplete Cholesky with dual threshold.
void factorize(const MatrixType &amat)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
void setShift(Scalar shift)
Set the initial shift parameter.
Index nonZeros() const
Definition: SparseMatrix.h:246
Matrix< Scalar, Dynamic, 1 > ScalarType
const Index * outerIndexPtr() const
Definition: SparseMatrix.h:149
PermutationMatrix< Dynamic, Dynamic, Index > PermutationType
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
Derived & setZero(Index size)
Transpose< PermutationBase > inverse() const
Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
Definition: IncompleteLUT.h:28
void rhs(const real_t *x, real_t *f)
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector.
std::vector< std::list< Index > > VectorList
void compute(const MatrixType &matrix)
Matrix< Index, Dynamic, 1 > IndexType
IncompleteCholesky(const MatrixType &matrix)
const internal::solve_retval< IncompleteCholesky, Rhs > solve(const MatrixBase< Rhs > &b) const
const Index * innerIndexPtr() const
Definition: SparseMatrix.h:140
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
ComputationInfo info() const
Reports whether previous computation was successful.
ColXpr col(Index i)
Definition: BlockMethods.h:708
NumTraits< Scalar >::Real RealScalar
void updateList(const IdxType &colPtr, IdxType &rowIdx, SclType &vals, const Index &col, const Index &jk, IndexType &firstElt, VectorList &listCol)
#define eigen_assert(x)
SparseMatrix< Scalar, ColMajor > m_L
const Scalar * valuePtr() const
Definition: SparseMatrix.h:131
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:119
internal::traits< Derived >::Index Index


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:41