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
17namespace Eigen {
18/**
19  * \brief Modified Incomplete Cholesky with dual threshold
20  *
21  * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
22  *              Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
23  *
24  * \tparam Scalar the scalar type of the input matrices
25  * \tparam _UpLo The triangular part that will be used for the computations. It can be Lower
26    *               or Upper. Default is Lower.
27  * \tparam _OrderingType The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<int>,
28  *                       unless EIGEN_MPL2_ONLY is defined, in which case the default is NaturalOrdering<int>.
29  *
30  * \implsparsesolverconcept
31  *
32  * It performs the following incomplete factorization: \f$ S P A P' S \approx L L' \f$
33  * where L is a lower triangular factor, S is a diagonal scaling matrix, and P is a
34  * fill-in reducing permutation as computed by the ordering method.
35  *
36  * \b Shifting \b strategy: Let \f$ B = S P A P' S \f$  be the scaled matrix on which the factorization is carried out,
37  * and \f$ \beta \f$ be the minimum value of the diagonal. If \f$ \beta > 0 \f$ then, the factorization is directly performed
38  * on the matrix B. Otherwise, the factorization is performed on the shifted matrix \f$ B + (\sigma+|\beta| I \f$ where
39  * \f$ \sigma \f$ is the initial shift value as returned and set by setInitialShift() method. The default value is \f$ \sigma = 10^{-3} \f$.
40  * If the factorization fails, then the shift in doubled until it succeed or a maximum of ten attempts. If it still fails, as returned by
41  * the info() method, then you can either increase the initial shift, or better use another preconditioning technique.
42  *
43  */
44template <typename Scalar, int _UpLo = Lower, typename _OrderingType =
45#ifndef EIGEN_MPL2_ONLY
46AMDOrdering<int>
47#else
48NaturalOrdering<int>
49#endif
50>
51class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> >
52{
53  protected:
54    typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base;
55    using Base::m_isInitialized;
56  public:
57    typedef typename NumTraits<Scalar>::Real RealScalar;
58    typedef _OrderingType OrderingType;
59    typedef typename OrderingType::PermutationType PermutationType;
60    typedef typename PermutationType::StorageIndex StorageIndex;
61    typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType;
62    typedef Matrix<Scalar,Dynamic,1> VectorSx;
63    typedef Matrix<RealScalar,Dynamic,1> VectorRx;
64    typedef Matrix<StorageIndex,Dynamic, 1> VectorIx;
65    typedef std::vector<std::list<StorageIndex> > VectorList;
66    enum { UpLo = _UpLo };
67    enum {
68      ColsAtCompileTime = Dynamic,
69      MaxColsAtCompileTime = Dynamic
70    };
71  public:
72
73    /** Default constructor leaving the object in a partly non-initialized stage.
74      *
75      * You must call compute() or the pair analyzePattern()/factorize() to make it valid.
76      *
77      * \sa IncompleteCholesky(const MatrixType&)
78      */
79    IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {}
80
81    /** Constructor computing the incomplete factorization for the given matrix \a matrix.
82      */
83    template<typename MatrixType>
84    IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false)
85    {
86      compute(matrix);
87    }
88
89    /** \returns number of rows of the factored matrix */
90    Index rows() const { return m_L.rows(); }
91
92    /** \returns number of columns of the factored matrix */
93    Index cols() const { return m_L.cols(); }
94
95
96    /** \brief Reports whether previous computation was successful.
97      *
98      * It triggers an assertion if \c *this has not been initialized through the respective constructor,
99      * or a call to compute() or analyzePattern().
100      *
101      * \returns \c Success if computation was successful,
102      *          \c NumericalIssue if the matrix appears to be negative.
103      */
104    ComputationInfo info() const
105    {
106      eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized.");
107      return m_info;
108    }
109
110    /** \brief Set the initial shift parameter \f$ \sigma \f$.
111      */
112    void setInitialShift(RealScalar shift) { m_initialShift = shift; }
113
114    /** \brief Computes the fill reducing permutation vector using the sparsity pattern of \a mat
115      */
116    template<typename MatrixType>
117    void analyzePattern(const MatrixType& mat)
118    {
119      OrderingType ord;
120      PermutationType pinv;
121      ord(mat.template selfadjointView<UpLo>(), pinv);
122      if(pinv.size()>0) m_perm = pinv.inverse();
123      else              m_perm.resize(0);
124      m_L.resize(mat.rows(), mat.cols());
125      m_analysisIsOk = true;
126      m_isInitialized = true;
127      m_info = Success;
128    }
129
130    /** \brief Performs the numerical factorization of the input matrix \a mat
131      *
132      * The method analyzePattern() or compute() must have been called beforehand
133      * with a matrix having the same pattern.
134      *
135      * \sa compute(), analyzePattern()
136      */
137    template<typename MatrixType>
138    void factorize(const MatrixType& mat);
139
140    /** Computes or re-computes the incomplete Cholesky factorization of the input matrix \a mat
141      *
142      * It is a shortcut for a sequential call to the analyzePattern() and factorize() methods.
143      *
144      * \sa analyzePattern(), factorize()
145      */
146    template<typename MatrixType>
147    void compute(const MatrixType& mat)
148    {
149      analyzePattern(mat);
150      factorize(mat);
151    }
152
153    // internal
154    template<typename Rhs, typename Dest>
155    void _solve_impl(const Rhs& b, Dest& x) const
156    {
157      eigen_assert(m_factorizationIsOk && "factorize() should be called first");
158      if (m_perm.rows() == b.rows())  x = m_perm * b;
159      else                            x = b;
160      x = m_scale.asDiagonal() * x;
161      x = m_L.template triangularView<Lower>().solve(x);
162      x = m_L.adjoint().template triangularView<Upper>().solve(x);
163      x = m_scale.asDiagonal() * x;
164      if (m_perm.rows() == b.rows())
165        x = m_perm.inverse() * x;
166    }
167
168    /** \returns the sparse lower triangular factor L */
169    const FactorType& matrixL() const { eigen_assert("m_factorizationIsOk"); return m_L; }
170
171    /** \returns a vector representing the scaling factor S */
172    const VectorRx& scalingS() const { eigen_assert("m_factorizationIsOk"); return m_scale; }
173
174    /** \returns the fill-in reducing permutation P (can be empty for a natural ordering) */
175    const PermutationType& permutationP() const { eigen_assert("m_analysisIsOk"); return m_perm; }
176
177  protected:
178    FactorType m_L;              // The lower part stored in CSC
179    VectorRx m_scale;            // The vector for scaling the matrix
180    RealScalar m_initialShift;   // The initial shift parameter
181    bool m_analysisIsOk;
182    bool m_factorizationIsOk;
183    ComputationInfo m_info;
184    PermutationType m_perm;
185
186  private:
187    inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol);
188};
189
190// Based on the following paper:
191//   C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with
192//   Limited memory, SIAM J. Sci. Comput.  21(1), pp. 24-45, 1999
193//   http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
194template<typename Scalar, int _UpLo, typename OrderingType>
195template<typename _MatrixType>
196void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat)
197{
198  using std::sqrt;
199  eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
200
201  // 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
202
203  // Apply the fill-reducing permutation computed in analyzePattern()
204  if (m_perm.rows() == mat.rows() ) // To detect the null permutation
205  {
206    // The temporary is needed to make sure that the diagonal entry is properly sorted
207    FactorType tmp(mat.rows(), mat.cols());
208    tmp = mat.template selfadjointView<_UpLo>().twistedBy(m_perm);
209    m_L.template selfadjointView<Lower>() = tmp.template selfadjointView<Lower>();
210  }
211  else
212  {
213    m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>();
214  }
215
216  Index n = m_L.cols();
217  Index nnz = m_L.nonZeros();
218  Map<VectorSx> vals(m_L.valuePtr(), nnz);         //values
219  Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz);  //Row indices
220  Map<VectorIx> colPtr( m_L.outerIndexPtr(), n+1); // Pointer to the beginning of each row
221  VectorIx firstElt(n-1); // for each j, points to the next entry in vals that will be used in the factorization
222  VectorList listCol(n);  // listCol(j) is a linked list of columns to update column j
223  VectorSx col_vals(n);   // Store a  nonzero values in each column
224  VectorIx col_irow(n);   // Row indices of nonzero elements in each column
225  VectorIx col_pattern(n);
226  col_pattern.fill(-1);
227  StorageIndex col_nnz;
228
229
230  // Computes the scaling factors
231  m_scale.resize(n);
232  m_scale.setZero();
233  for (Index j = 0; j < n; j++)
234    for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
235    {
236      m_scale(j) += numext::abs2(vals(k));
237      if(rowIdx[k]!=j)
238        m_scale(rowIdx[k]) += numext::abs2(vals(k));
239    }
240
241  m_scale = m_scale.cwiseSqrt().cwiseSqrt();
242
243  for (Index j = 0; j < n; ++j)
244    if(m_scale(j)>(std::numeric_limits<RealScalar>::min)())
245      m_scale(j) = RealScalar(1)/m_scale(j);
246    else
247      m_scale(j) = 1;
248
249  // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster)
250
251  // Scale and compute the shift for the matrix
252  RealScalar mindiag = NumTraits<RealScalar>::highest();
253  for (Index j = 0; j < n; j++)
254  {
255    for (Index k = colPtr[j]; k < colPtr[j+1]; k++)
256      vals[k] *= (m_scale(j)*m_scale(rowIdx[k]));
257    eigen_internal_assert(rowIdx[colPtr[j]]==j && "IncompleteCholesky: only the lower triangular part must be stored");
258    mindiag = numext::mini(numext::real(vals[colPtr[j]]), mindiag);
259  }
260
261  FactorType L_save = m_L;
262
263  RealScalar shift = 0;
264  if(mindiag <= RealScalar(0.))
265    shift = m_initialShift - mindiag;
266
267  m_info = NumericalIssue;
268
269  // Try to perform the incomplete factorization using the current shift
270  int iter = 0;
271  do
272  {
273    // Apply the shift to the diagonal elements of the matrix
274    for (Index j = 0; j < n; j++)
275      vals[colPtr[j]] += shift;
276
277    // jki version of the Cholesky factorization
278    Index j=0;
279    for (; j < n; ++j)
280    {
281      // Left-looking factorization of the j-th column
282      // First, load the j-th column into col_vals
283      Scalar diag = vals[colPtr[j]];  // It is assumed that only the lower part is stored
284      col_nnz = 0;
285      for (Index i = colPtr[j] + 1; i < colPtr[j+1]; i++)
286      {
287        StorageIndex l = rowIdx[i];
288        col_vals(col_nnz) = vals[i];
289        col_irow(col_nnz) = l;
290        col_pattern(l) = col_nnz;
291        col_nnz++;
292      }
293      {
294        typename std::list<StorageIndex>::iterator k;
295        // Browse all previous columns that will update column j
296        for(k = listCol[j].begin(); k != listCol[j].end(); k++)
297        {
298          Index jk = firstElt(*k); // First element to use in the column
299          eigen_internal_assert(rowIdx[jk]==j);
300          Scalar v_j_jk = numext::conj(vals[jk]);
301
302          jk += 1;
303          for (Index i = jk; i < colPtr[*k+1]; i++)
304          {
305            StorageIndex l = rowIdx[i];
306            if(col_pattern[l]<0)
307            {
308              col_vals(col_nnz) = vals[i] * v_j_jk;
309              col_irow[col_nnz] = l;
310              col_pattern(l) = col_nnz;
311              col_nnz++;
312            }
313            else
314              col_vals(col_pattern[l]) -= vals[i] * v_j_jk;
315          }
316          updateList(colPtr,rowIdx,vals, *k, jk, firstElt, listCol);
317        }
318      }
319
320      // Scale the current column
321      if(numext::real(diag) <= 0)
322      {
323        if(++iter>=10)
324          return;
325
326        // increase shift
327        shift = numext::maxi(m_initialShift,RealScalar(2)*shift);
328        // restore m_L, col_pattern, and listCol
329        vals = Map<const VectorSx>(L_save.valuePtr(), nnz);
330        rowIdx = Map<const VectorIx>(L_save.innerIndexPtr(), nnz);
331        colPtr = Map<const VectorIx>(L_save.outerIndexPtr(), n+1);
332        col_pattern.fill(-1);
333        for(Index i=0; i<n; ++i)
334          listCol[i].clear();
335
336        break;
337      }
338
339      RealScalar rdiag = sqrt(numext::real(diag));
340      vals[colPtr[j]] = rdiag;
341      for (Index k = 0; k<col_nnz; ++k)
342      {
343        Index i = col_irow[k];
344        //Scale
345        col_vals(k) /= rdiag;
346        //Update the remaining diagonals with col_vals
347        vals[colPtr[i]] -= numext::abs2(col_vals(k));
348      }
349      // Select the largest p elements
350      // p is the original number of elements in the column (without the diagonal)
351      Index p = colPtr[j+1] - colPtr[j] - 1 ;
352      Ref<VectorSx> cvals = col_vals.head(col_nnz);
353      Ref<VectorIx> cirow = col_irow.head(col_nnz);
354      internal::QuickSplit(cvals,cirow, p);
355      // Insert the largest p elements in the matrix
356      Index cpt = 0;
357      for (Index i = colPtr[j]+1; i < colPtr[j+1]; i++)
358      {
359        vals[i] = col_vals(cpt);
360        rowIdx[i] = col_irow(cpt);
361        // restore col_pattern:
362        col_pattern(col_irow(cpt)) = -1;
363        cpt++;
364      }
365      // Get the first smallest row index and put it after the diagonal element
366      Index jk = colPtr(j)+1;
367      updateList(colPtr,rowIdx,vals,j,jk,firstElt,listCol);
368    }
369
370    if(j==n)
371    {
372      m_factorizationIsOk = true;
373      m_info = Success;
374    }
375  } while(m_info!=Success);
376}
377
378template<typename Scalar, int _UpLo, typename OrderingType>
379inline 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)
380{
381  if (jk < colPtr(col+1) )
382  {
383    Index p = colPtr(col+1) - jk;
384    Index minpos;
385    rowIdx.segment(jk,p).minCoeff(&minpos);
386    minpos += jk;
387    if (rowIdx(minpos) != rowIdx(jk))
388    {
389      //Swap
390      std::swap(rowIdx(jk),rowIdx(minpos));
391      std::swap(vals(jk),vals(minpos));
392    }
393    firstElt(col) = internal::convert_index<StorageIndex,Index>(jk);
394    listCol[rowIdx(jk)].push_back(internal::convert_index<StorageIndex,Index>(col));
395  }
396}
397
398} // end namespace Eigen
399
400#endif
401