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
16namespace Eigen {
17/**
18 * \brief Modified Incomplete Cholesky with dual threshold
19 *
20 * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
21 *              Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
22 *
23 * \tparam _MatrixType The type of the sparse matrix. It should be a symmetric
24 *                     matrix. It is advised to give  a row-oriented sparse matrix
25 * \tparam _UpLo The triangular part of the matrix to reference.
26 * \tparam _OrderingType
27 */
28
29template <typename Scalar, int _UpLo = Lower, typename _OrderingType = NaturalOrdering<int> >
30class IncompleteCholesky : internal::noncopyable
31{
32  public:
33    typedef SparseMatrix<Scalar,ColMajor> MatrixType;
34    typedef _OrderingType OrderingType;
35    typedef typename MatrixType::RealScalar RealScalar;
36    typedef typename MatrixType::Index Index;
37    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
38    typedef Matrix<Scalar,Dynamic,1> ScalarType;
39    typedef Matrix<Index,Dynamic, 1> IndexType;
40    typedef std::vector<std::list<Index> > VectorList;
41    enum { UpLo = _UpLo };
42  public:
43    IncompleteCholesky() : m_shift(1),m_factorizationIsOk(false) {}
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
54    /** \brief Reports whether previous computation was successful.
55      *
56      * \returns \c Success if computation was succesful,
57      *          \c NumericalIssue if the matrix appears to be negative.
58      */
59    ComputationInfo info() const
60    {
61      eigen_assert(m_isInitialized && "IncompleteLLT is not initialized.");
62      return m_info;
63    }
64
65    /**
66     * \brief Set the initial shift parameter
67     */
68    void setShift( Scalar shift) { m_shift = shift; }
69
70    /**
71    * \brief Computes the fill reducing permutation vector.
72    */
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
119    bool m_analysisIsOk;
120    bool m_factorizationIsOk;
121    bool m_isInitialized;
122    ComputationInfo m_info;
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
130template<typename Scalar, int _UpLo, typename OrderingType>
131template<typename _MatrixType>
132void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
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";
206      m_info = NumericalIssue;
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
239template<typename Scalar, int _UpLo, typename OrderingType>
240template <typename IdxType, typename SclType>
241inline 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}
259namespace internal {
260
261template<typename _Scalar, int _UpLo, typename OrderingType, typename Rhs>
262struct solve_retval<IncompleteCholesky<_Scalar,  _UpLo, OrderingType>, Rhs>
263  : solve_retval_base<IncompleteCholesky<_Scalar, _UpLo, OrderingType>, Rhs>
264{
265  typedef IncompleteCholesky<_Scalar, _UpLo, OrderingType> Dec;
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
279