Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
IncompleteCholesky.h
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) 2015 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_INCOMPLETE_CHOlESKY_H
12 #define EIGEN_INCOMPLETE_CHOlESKY_H
13 
14 #include <vector>
15 #include <list>
16 
17 #include "./InternalHeaderCheck.h"
18 
19 namespace Eigen {
46 template <typename Scalar, int UpLo_ = Lower, typename OrderingType_ = AMDOrdering<int> >
47 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,UpLo_,OrderingType_> >
48 {
49  protected:
51  using Base::m_isInitialized;
52  public:
53  typedef typename NumTraits<Scalar>::Real RealScalar;
54  typedef OrderingType_ OrderingType;
55  typedef typename OrderingType::PermutationType PermutationType;
56  typedef typename PermutationType::StorageIndex StorageIndex;
61  typedef std::vector<std::list<StorageIndex> > VectorList;
62  enum { UpLo = UpLo_ };
63  enum {
64  ColsAtCompileTime = Dynamic,
65  MaxColsAtCompileTime = Dynamic
66  };
67  public:
68 
75  IncompleteCholesky() : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) {}
76 
79  template<typename MatrixType>
80  IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false)
81  {
82  compute(matrix);
83  }
84 
86  EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_L.rows(); }
87 
89  EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_L.cols(); }
90 
91 
101  {
102  eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
103  return m_info;
104  }
105 
108  void setInitialShift(RealScalar shift) { m_initialShift = shift; }
109 
112  template<typename MatrixType>
113  void analyzePattern(const MatrixType& mat)
114  {
115  OrderingType ord;
116  PermutationType pinv;
117  ord(mat.template selfadjointView<UpLo>(), pinv);
118  if(pinv.size()>0) m_perm = pinv.inverse();
119  else m_perm.resize(0);
120  m_L.resize(mat.rows(), mat.cols());
121  m_analysisIsOk = true;
122  m_isInitialized = true;
123  m_info = Success;
124  }
125 
133  template<typename MatrixType>
134  void factorize(const MatrixType& mat);
135 
142  template<typename MatrixType>
143  void compute(const MatrixType& mat)
144  {
145  analyzePattern(mat);
146  factorize(mat);
147  }
148 
149  // internal
150  template<typename Rhs, typename Dest>
151  void _solve_impl(const Rhs& b, Dest& x) const
152  {
153  eigen_assert(m_factorizationIsOk && "factorize() should be called first");
154  if (m_perm.rows() == b.rows()) x = m_perm * b;
155  else x = b;
156  x = m_scale.asDiagonal() * x;
157  x = m_L.template triangularView<Lower>().solve(x);
158  x = m_L.adjoint().template triangularView<Upper>().solve(x);
159  x = m_scale.asDiagonal() * x;
160  if (m_perm.rows() == b.rows())
161  x = m_perm.inverse() * x;
162  }
163 
165  const FactorType& matrixL() const { eigen_assert(m_factorizationIsOk && "factorize() should be called first"); return m_L; }
166 
168  const VectorRx& scalingS() const { eigen_assert(m_factorizationIsOk && "factorize() should be called first"); return m_scale; }
169 
171  const PermutationType& permutationP() const { eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); return m_perm; }
172 
173  protected:
174  FactorType m_L; // The lower part stored in CSC
175  VectorRx m_scale; // The vector for scaling the matrix
176  RealScalar m_initialShift; // The initial shift parameter
177  bool m_analysisIsOk;
178  bool m_factorizationIsOk;
179  ComputationInfo m_info;
180  PermutationType m_perm;
181 
182  private:
183  inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
184 };
185 
186 // Based on the following paper:
187 // C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
188 // Limited memory, SIAM J. Sci. Comput. 21(1), pp. 24-45, 1999
189 // http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
190 template<typename Scalar, int UpLo_, typename OrderingType>
191 template<typename MatrixType_>
193 {
194  using std::sqrt;
195  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
196 
197  // Dropping strategy : 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
198 
199  // Apply the fill-reducing permutation computed in analyzePattern()
200  if (m_perm.rows() == mat.rows() ) // To detect the null permutation
201  {
202  // The temporary is needed to make sure that the diagonal entry is properly sorted
203  FactorType tmp(mat.rows(), mat.cols());
204  tmp = mat.template selfadjointView<UpLo_>().twistedBy(m_perm);
205  m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
206  }
207  else
208  {
209  m_L.template selfadjointView<Lower>() = mat.template selfadjointView<UpLo_>();
210  }
211 
212  Index n = m_L.cols();
213  Index nnz = m_L.nonZeros();
214  Map<VectorSx> vals(m_L.valuePtr(), nnz); //values
215  Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
216  Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
217  VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
218  VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
219  VectorSx col_vals(n); // Store a nonzero values in each column
220  VectorIx col_irow(n); // Row indices of nonzero elements in each column
221  VectorIx col_pattern(n);
222  col_pattern.fill(-1);
223  StorageIndex col_nnz;
224 
225 
226  // Computes the scaling factors
227  m_scale.resize(n);
228  m_scale.setZero();
229  for (Index j = 0; j < n; j++)
230  for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
231  {
232  m_scale(j) += numext::abs2(vals(k));
233  if(rowIdx[k]!=j)
234  m_scale(rowIdx[k]) += numext::abs2(vals(k));
235  }
236 
237  m_scale = m_scale.cwiseSqrt().cwiseSqrt();
238 
239  for (Index j = 0; j < n; ++j)
240  if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
241  m_scale(j) = RealScalar(1)/m_scale(j);
242  else
243  m_scale(j) = 1;
244 
245  // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
246 
247  // Scale and compute the shift for the matrix
248  RealScalar mindiag = NumTraits<RealScalar>::highest();
249  for (Index j = 0; j < n; j++)
250  {
251  for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
252  vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
253  eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
254  mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
255  }
256 
257  FactorType L_save = m_L;
258 
259  RealScalar shift = 0;
260  if(mindiag <= RealScalar(0.))
261  shift = m_initialShift - mindiag;
262 
263  m_info = NumericalIssue;
264 
265  // Try to perform the incomplete factorization using the current shift
266  int iter = 0;
267  do
268  {
269  // Apply the shift to the diagonal elements of the matrix
270  for (Index j = 0; j < n; j++)
271  vals[colPtr[j]] += shift;
272 
273  // jki version of the Cholesky factorization
274  Index j=0;
275  for (; j < n; ++j)
276  {
277  // Left-looking factorization of the j-th column
278  // First, load the j-th column into col_vals
279  Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
280  col_nnz = 0;
281  for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
282  {
283  StorageIndex l = rowIdx[i];
284  col_vals(col_nnz) = vals[i];
285  col_irow(col_nnz) = l;
286  col_pattern(l) = col_nnz;
287  col_nnz++;
288  }
289  {
290  typename std::list<StorageIndex>::iterator k;
291  // Browse all previous columns that will update column j
292  for(k = listCol[j].begin(); k != listCol[j].end(); k++)
293  {
294  Index jk = firstElt(*k); // First element to use in the column
295  eigen_internal_assert(rowIdx[jk]==j);
296  Scalar v_j_jk = numext::conj(vals[jk]);
297 
298  jk += 1;
299  for (Index i = jk; i < colPtr[*k+1]; i++)
300  {
301  StorageIndex l = rowIdx[i];
302  if(col_pattern[l]<0)
303  {
304  col_vals(col_nnz) = vals[i] * v_j_jk;
305  col_irow[col_nnz] = l;
306  col_pattern(l) = col_nnz;
307  col_nnz++;
308  }
309  else
310  col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
311  }
312  updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
313  }
314  }
315 
316  // Scale the current column
317  if(numext::real(diag) <= 0)
318  {
319  if(++iter>=10)
320  return;
321 
322  // increase shift
323  shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
324  // restore m_L, col_pattern, and listCol
325  vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
326  rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
327  colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
328  col_pattern.fill(-1);
329  for(Index i=0; i<n; ++i)
330  listCol[i].clear();
331 
332  break;
333  }
334 
335  RealScalar rdiag = sqrt(numext::real(diag));
336  vals[colPtr[j]] = rdiag;
337  for (Index k = 0; k<col_nnz; ++k)
338  {
339  Index i = col_irow[k];
340  //Scale
341  col_vals(k) /= rdiag;
342  //Update the remaining diagonals with col_vals
343  vals[colPtr[i]] -= numext::abs2(col_vals(k));
344  }
345  // Select the largest p elements
346  // p is the original number of elements in the column (without the diagonal)
347  Index p = colPtr[j+1] - colPtr[j] - 1 ;
348  Ref<VectorSx> cvals = col_vals.head(col_nnz);
349  Ref<VectorIx> cirow = col_irow.head(col_nnz);
350  internal::QuickSplit(cvals,cirow, p);
351  // Insert the largest p elements in the matrix
352  Index cpt = 0;
353  for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
354  {
355  vals[i] = col_vals(cpt);
356  rowIdx[i] = col_irow(cpt);
357  // restore col_pattern:
358  col_pattern(col_irow(cpt)) = -1;
359  cpt++;
360  }
361  // Get the first smallest row index and put it after the diagonal element
362  Index jk = colPtr(j)+1;
363  updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
364  }
365 
366  if(j==n)
367  {
368  m_factorizationIsOk = true;
369  m_info = Success;
370  }
371  } while(m_info!=Success);
372 }
373 
374 template<typename Scalar, int UpLo_, typename OrderingType>
375 inline void IncompleteCholesky<Scalar,UpLo_, OrderingType>::updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol)
376 {
377  if (jk < colPtr(col+1) )
378  {
379  Index p = colPtr(col+1) - jk;
380  Index minpos;
381  rowIdx.segment(jk,p).minCoeff(&minpos);
382  minpos += jk;
383  if (rowIdx(minpos) != rowIdx(jk))
384  {
385  //Swap
386  std::swap(rowIdx(jk),rowIdx(minpos));
387  std::swap(vals(jk),vals(minpos));
388  }
389  firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
390  listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
391  }
392 }
393 
394 } // end namespace Eigen
395 
396 #endif
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:48
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:100
IncompleteCholesky(const MatrixType &matrix)
Definition: IncompleteCholesky.h:80
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:86
void factorize(const MatrixType &mat)
Performs the numerical factorization of the input matrix mat.
const VectorRx & scalingS() const
Definition: IncompleteCholesky.h:168
IncompleteCholesky()
Definition: IncompleteCholesky.h:75
const PermutationType & permutationP() const
Definition: IncompleteCholesky.h:171
void compute(const MatrixType &mat)
Definition: IncompleteCholesky.h:143
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: IncompleteCholesky.h:89
const FactorType & matrixL() const
Definition: IncompleteCholesky.h:165
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector using the sparsity pattern of mat.
Definition: IncompleteCholesky.h:113
void setInitialShift(RealScalar shift)
Set the initial shift parameter .
Definition: IncompleteCholesky.h:108
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:285
Index cols() const
Definition: SparseMatrix.h:142
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:626
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const int Dynamic
Definition: Constants.h:24
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231