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 //
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
12 #include "Eigen/src/IterativeLinearSolvers/IncompleteLUT.h"
13 #include <Eigen/OrderingMethods>
14 #include <list>
15 
16 namespace Eigen {
29 template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> >
30 class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
31 {
32  protected:
33  typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
34  using Base::m_isInitialized;
35  public:
36  typedef typename NumTraits<Scalar>::Real RealScalar;
37  typedef _OrderingType OrderingType;
38  typedef typename OrderingType::PermutationType PermutationType;
39  typedef typename PermutationType::StorageIndex StorageIndex;
40  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
41  typedef FactorType MatrixType;
42  typedef Matrix<Scalar,Dynamic,1> VectorSx;
43  typedef Matrix<RealScalar,Dynamic,1> VectorRx;
44  typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
45  typedef std::vector<std::list<StorageIndex> > VectorList;
46  enum { UpLo = _UpLo };
47  public:
48  IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
49 
50  template<typename MatrixType>
51  IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
52  {
53  compute(matrix);
54  }
55 
56  Index rows() const { return m_L.rows(); }
57 
58  Index cols() const { return m_L.cols(); }
59 
60 
66  ComputationInfo info() const
67  {
68  eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
69  return m_info;
70  }
71 
75  void setInitialShift(RealScalar shift) { m_initialShift = shift; }
76 
80  template<typename MatrixType>
81  void analyzePattern(const MatrixType& mat)
82  {
83  OrderingType ord;
84  PermutationType pinv;
85  ord(mat.template selfadjointView<UpLo>(), pinv);
86  if(pinv.size()>0) m_perm = pinv.inverse();
87  else m_perm.resize(0);
88  m_analysisIsOk = true;
89  }
90 
91  template<typename MatrixType>
92  void factorize(const MatrixType& amat);
93 
94  template<typename MatrixType>
95  void compute(const MatrixType& matrix)
96  {
97  analyzePattern(matrix);
98  factorize(matrix);
99  }
100 
101  template<typename Rhs, typename Dest>
102  void _solve_impl(const Rhs& b, Dest& x) const
103  {
104  eigen_assert(m_factorizationIsOk && "factorize() should be called first");
105  if (m_perm.rows() == b.rows()) x = m_perm * b;
106  else x = b;
107  x = m_scale.asDiagonal() * x;
108  x = m_L.template triangularView<Lower>().solve(x);
109  x = m_L.adjoint().template triangularView<Upper>().solve(x);
110  x = m_scale.asDiagonal() * x;
111  if (m_perm.rows() == b.rows())
112  x = m_perm.inverse() * x;
113 
114  }
115 
116  protected:
117  FactorType m_L; // The lower part stored in CSC
118  VectorRx m_scale; // The vector for scaling the matrix
119  RealScalar m_initialShift; // The initial shift parameter
120  bool m_analysisIsOk;
121  bool m_factorizationIsOk;
122  ComputationInfo m_info;
123  PermutationType m_perm;
124 
125  private:
126  inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
127 };
128 
129 template<typename Scalar, int _UpLo, typename OrderingType>
130 template<typename _MatrixType>
131 void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
132 {
133  using std::sqrt;
134  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
135 
136  // 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
137 
138  m_L.resize(mat.rows(), mat.cols());
139 
140  // Apply the fill-reducing permutation computed in analyzePattern()
141  if (m_perm.rows() == mat.rows() ) // To detect the null permutation
142  {
143  // The temporary is needed to make sure that the diagonal entry is properly sorted
144  FactorType tmp(mat.rows(), mat.cols());
145  tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
146  m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
147  }
148  else
149  {
150  m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
151  }
152 
153  Index n = m_L.cols();
154  Index nnz = m_L.nonZeros();
155  Map<VectorSx> vals(m_L.valuePtr(), nnz); //values
156  Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices
157  Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
158  VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
159  VectorList listCol(n); // listCol(j) is a linked list of columns to update column j
160  VectorSx col_vals(n); // Store a nonzero values in each column
161  VectorIx col_irow(n); // Row indices of nonzero elements in each column
162  VectorIx col_pattern(n);
163  col_pattern.fill(-1);
164  StorageIndex col_nnz;
165 
166 
167  // Computes the scaling factors
168  m_scale.resize(n);
169  m_scale.setZero();
170  for (Index j = 0; j < n; j++)
171  for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
172  {
173  m_scale(j) += numext::abs2(vals(k));
174  if(rowIdx[k]!=j)
175  m_scale(rowIdx[k]) += numext::abs2(vals(k));
176  }
177 
178  m_scale = m_scale.cwiseSqrt().cwiseSqrt();
179 
180  // Scale and compute the shift for the matrix
181  RealScalar mindiag = NumTraits<RealScalar>::highest();
182  for (Index j = 0; j < n; j++)
183  {
184  for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
185  vals[k] /= (m_scale(j)*m_scale(rowIdx[k]));
186  eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
187  mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
188  }
189 
190  RealScalar shift = 0;
191  if(mindiag <= RealScalar(0.))
192  shift = m_initialShift - mindiag;
193 
194  // Apply the shift to the diagonal elements of the matrix
195  for (Index j = 0; j < n; j++)
196  vals[colPtr[j]] += shift;
197 
198  // jki version of the Cholesky factorization
199  for (Index j=0; j < n; ++j)
200  {
201  // Left-looking factorization of the j-th column
202  // First, load the j-th column into col_vals
203  Scalar diag = vals[colPtr[j]]; // It is assumed that only the lower part is stored
204  col_nnz = 0;
205  for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
206  {
207  StorageIndex l = rowIdx[i];
208  col_vals(col_nnz) = vals[i];
209  col_irow(col_nnz) = l;
210  col_pattern(l) = col_nnz;
211  col_nnz++;
212  }
213  {
214  typename std::list<StorageIndex>::iterator k;
215  // Browse all previous columns that will update column j
216  for(k = listCol[j].begin(); k != listCol[j].end(); k++)
217  {
218  Index jk = firstElt(*k); // First element to use in the column
219  eigen_internal_assert(rowIdx[jk]==j);
220  Scalar v_j_jk = numext::conj(vals[jk]);
221 
222  jk += 1;
223  for (Index i = jk; i < colPtr[*k+1]; i++)
224  {
225  StorageIndex l = rowIdx[i];
226  if(col_pattern[l]<0)
227  {
228  col_vals(col_nnz) = vals[i] * v_j_jk;
229  col_irow[col_nnz] = l;
230  col_pattern(l) = col_nnz;
231  col_nnz++;
232  }
233  else
234  col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
235  }
236  updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
237  }
238  }
239 
240  // Scale the current column
241  if(numext::real(diag) <= 0)
242  {
243  std::cerr << "\nNegative diagonal during Incomplete factorization at position " << j << " (value = " << diag << ")\n";
244  m_info = NumericalIssue;
245  return;
246  }
247 
248  RealScalar rdiag = sqrt(numext::real(diag));
249  vals[colPtr[j]] = rdiag;
250  for (Index k = 0; k<col_nnz; ++k)
251  {
252  Index i = col_irow[k];
253  //Scale
254  col_vals(k) /= rdiag;
255  //Update the remaining diagonals with col_vals
256  vals[colPtr[i]] -= numext::abs2(col_vals(k));
257  }
258  // Select the largest p elements
259  // p is the original number of elements in the column (without the diagonal)
260  Index p = colPtr[j+1] - colPtr[j] - 1 ;
261  Ref<VectorSx> cvals = col_vals.head(col_nnz);
262  Ref<VectorIx> cirow = col_irow.head(col_nnz);
263  internal::QuickSplit(cvals,cirow, p);
264  // Insert the largest p elements in the matrix
265  Index cpt = 0;
266  for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
267  {
268  vals[i] = col_vals(cpt);
269  rowIdx[i] = col_irow(cpt);
270  // restore col_pattern:
271  col_pattern(col_irow(cpt)) = -1;
272  cpt++;
273  }
274  // Get the first smallest row index and put it after the diagonal element
275  Index jk = colPtr(j)+1;
276  updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
277  }
278  m_factorizationIsOk = true;
279  m_isInitialized = true;
280  m_info = Success;
281 }
282 
283 template<typename Scalar, int _UpLo, typename OrderingType>
284 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)
285 {
286  if (jk < colPtr(col+1) )
287  {
288  Index p = colPtr(col+1) - jk;
289  Index minpos;
290  rowIdx.segment(jk,p).minCoeff(&minpos);
291  minpos += jk;
292  if (rowIdx(minpos) != rowIdx(jk))
293  {
294  //Swap
295  std::swap(rowIdx(jk),rowIdx(minpos));
296  std::swap(vals(jk),vals(minpos));
297  }
298  firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
299  listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
300  }
301 }
302 
303 } // end namespace Eigen
304 
305 #endif
Modified Incomplete Cholesky with dual threshold.
Definition: IncompleteCholesky.h:30
Namespace containing all symbols from the Eigen library.
Definition: CXX11Meta.h:13
void setInitialShift(RealScalar shift)
Set the initial shift parameter.
Definition: IncompleteCholesky.h:75
void analyzePattern(const MatrixType &mat)
Computes the fill reducing permutation vector.
Definition: IncompleteCholesky.h:81
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteCholesky.h:66