Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseLDLT.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24
25 /*
26
27 NOTE: the _symbolic, and _numeric functions has been adapted from
28       the LDL library:
29
30 LDL Copyright (c) 2005 by Timothy A. Davis.  All Rights Reserved.
31
32 LDL License:
33
34     Your use or distribution of LDL or any modified version of
35     LDL implies that you agree to this License.
36
37     This library is free software; you can redistribute it and/or
38     modify it under the terms of the GNU Lesser General Public
39     License as published by the Free Software Foundation; either
40     version 2.1 of the License, or (at your option) any later version.
41
42     This library is distributed in the hope that it will be useful,
43     but WITHOUT ANY WARRANTY; without even the implied warranty of
44     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
45     Lesser General Public License for more details.
46
47     You should have received a copy of the GNU Lesser General Public
48     License along with this library; if not, write to the Free Software
49     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
50     USA
51
52     Permission is hereby granted to use or copy this program under the
53     terms of the GNU LGPL, provided that the Copyright, this License,
54     and the Availability of the original version is retained on all copies.
55     User documentation of any code that uses this code or any modified
56     version of this code must cite the Copyright, this License, the
57     Availability note, and "Used by permission." Permission to modify
58     the code and to distribute modified code is granted, provided the
59     Copyright, this License, and the Availability note are retained,
60     and a notice that the code was modified is included.
61  */
62
63 #ifndef EIGEN_SPARSELDLT_H
64 #define EIGEN_SPARSELDLT_H
65
66 /** \ingroup Sparse_Module
67   *
68   * \class SparseLDLT
69   *
70   * \brief LDLT Cholesky decomposition of a sparse matrix and associated features
71   *
72   * \param MatrixType the type of the matrix of which we are computing the LDLT Cholesky decomposition
73   *
74   * \sa class LDLT, class LDLT
75   */
76 template<typename MatrixType, int Backend = DefaultBackend>
77 class SparseLDLT
78 {
79   protected:
80     typedef typename MatrixType::Scalar Scalar;
81     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
82     typedef SparseMatrix<Scalar,LowerTriangular|UnitDiagBit> CholMatrixType;
83     typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType;
84
85     enum {
86       SupernodalFactorIsDirty      = 0x10000,
87       MatrixLIsDirty               = 0x20000
88     };
89
90   public:
91
92     /** Creates a dummy LDLT factorization object with flags \a flags. */
93     SparseLDLT(int flags = 0)
94       : m_flags(flags), m_status(0)
95     {
96       ei_assert((MatrixType::Flags&RowMajorBit)==0);
97       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
98     }
99
100     /** Creates a LDLT object and compute the respective factorization of \a matrix using
101       * flags \a flags. */
102     SparseLDLT(const MatrixType& matrix, int flags = 0)
103       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
104     {
105       ei_assert((MatrixType::Flags&RowMajorBit)==0);
106       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
107       compute(matrix);
108     }
109
110     /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
111       *
112       * Setting a value greater than zero speeds up computation, and yields to an imcomplete
113       * factorization with fewer non zero coefficients. Such approximate factors are especially
114       * useful to initialize an iterative solver.
115       *
116       * \warning if precision is greater that zero, the LDLT factorization is not guaranteed to succeed
117       * even if the matrix is positive definite.
118       *
119       * Note that the exact meaning of this parameter might depends on the actual
120       * backend. Moreover, not all backends support this feature.
121       *
122       * \sa precision() */
123     void setPrecision(RealScalar v) { m_precision = v; }
124
125     /** \returns the current precision.
126       *
127       * \sa setPrecision() */
128     RealScalar precision() const { return m_precision; }
129
130     /** Sets the flags. Possible values are:
131       *  - CompleteFactorization
132       *  - IncompleteFactorization
133       *  - MemoryEfficient          (hint to use the memory most efficient method offered by the backend)
134       *  - SupernodalMultifrontal   (implies a complete factorization if supported by the backend,
135       *                              overloads the MemoryEfficient flags)
136       *  - SupernodalLeftLooking    (implies a complete factorization  if supported by the backend,
137       *                              overloads the MemoryEfficient flags)
138       *
139       * \sa flags() */
140     void settagss(int f) { m_flags = f; }
141     /** \returns the current flags */
142     int flags() const { return m_flags; }
143
144     /** Computes/re-computes the LDLT factorization */
145     void compute(const MatrixType& matrix);
146
147     /** Perform a symbolic factorization */
148     void _symbolic(const MatrixType& matrix);
149     /** Perform the actual factorization using the previously
150       * computed symbolic factorization */
151     bool _numeric(const MatrixType& matrix);
152
153     /** \returns the lower triangular matrix L */
154     inline const CholMatrixType& matrixL(void) const { return m_matrix; }
155
156     /** \returns the coefficients of the diagonal matrix D */
157     inline VectorType vectorD(void) const { return m_diag; }
158
159     template<typename Derived>
160     bool solveInPlace(MatrixBase<Derived> &b) const;
161
162     /** \returns true if the factorization succeeded */
163     inline bool succeeded(void) const { return m_succeeded; }
164
165   protected:
166     CholMatrixType m_matrix;
167     VectorType m_diag;
168     VectorXi m_parent; // elimination tree
169     VectorXi m_nonZerosPerCol;
170 //     VectorXi m_w; // workspace
171     RealScalar m_precision;
172     int m_flags;
173     mutable int m_status;
174     bool m_succeeded;
175 };
176
177 /** Computes / recomputes the LDLT decomposition of matrix \a a
178   * using the default algorithm.
179   */
180 template<typename MatrixType, int Backend>
181 void SparseLDLT<MatrixType,Backend>::compute(const MatrixType& a)
182 {
183   _symbolic(a);
184   m_succeeded = _numeric(a);
185 }
186
187 template<typename MatrixType, int Backend>
188 void SparseLDLT<MatrixType,Backend>::_symbolic(const MatrixType& a)
189 {
190   assert(a.rows()==a.cols());
191   const int size = a.rows();
192   m_matrix.resize(size, size);
193   m_parent.resize(size);
194   m_nonZerosPerCol.resize(size);
195   int * tags = ei_aligned_stack_new(int, size);
196
197   const int* Ap = a._outerIndexPtr();
198   const int* Ai = a._innerIndexPtr();
199   int* Lp = m_matrix._outerIndexPtr();
200   const int* P = 0;
201   int* Pinv = 0;
202
203   if (P)
204   {
205     /* If P is present then compute Pinv, the inverse of P */
206     for (int k = 0; k < size; ++k)
207       Pinv[P[k]] = k;
208   }
209   for (int k = 0; k < size; ++k)
210   {
211     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
212     m_parent[k] = -1;             /* parent of k is not yet known */
213     tags[k] = k;                  /* mark node k as visited */
214     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
215     int kk = P ? P[k] : k;  /* kth original, or permuted, column */
216     int p2 = Ap[kk+1];
217     for (int p = Ap[kk]; p < p2; ++p)
218     {
219       /* A (i,k) is nonzero (original or permuted A) */
220       int i = Pinv ? Pinv[Ai[p]] : Ai[p];
221       if (i < k)
222       {
223         /* follow path from i to root of etree, stop at flagged node */
224         for (; tags[i] != k; i = m_parent[i])
225         {
226           /* find parent of i if not yet determined */
227           if (m_parent[i] == -1)
228             m_parent[i] = k;
229           ++m_nonZerosPerCol[i];        /* L (k,i) is nonzero */
230           tags[i] = k;                  /* mark i as visited */
231         }
232       }
233     }
234   }
235   /* construct Lp index array from m_nonZerosPerCol column counts */
236   Lp[0] = 0;
237   for (int k = 0; k < size; ++k)
238     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
239
240   m_matrix.resizeNonZeros(Lp[size]);
241   ei_aligned_stack_delete(int, tags, size);
242 }
243
244 template<typename MatrixType, int Backend>
245 bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
246 {
247   assert(a.rows()==a.cols());
248   const int size = a.rows();
249   assert(m_parent.size()==size);
250   assert(m_nonZerosPerCol.size()==size);
251
252   const int* Ap = a._outerIndexPtr();
253   const int* Ai = a._innerIndexPtr();
254   const Scalar* Ax = a._valuePtr();
255   const int* Lp = m_matrix._outerIndexPtr();
256   int* Li = m_matrix._innerIndexPtr();
257   Scalar* Lx = m_matrix._valuePtr();
258   m_diag.resize(size);
259
260   Scalar * y = ei_aligned_stack_new(Scalar, size);
261   int * pattern = ei_aligned_stack_new(int, size);
262   int * tags = ei_aligned_stack_new(int, size);
263
264   const int* P = 0;
265   const int* Pinv = 0;
266   bool ok = true;
267
268   for (int k = 0; k < size; ++k)
269   {
270     /* compute nonzero pattern of kth row of L, in topological order */
271     y[k] = 0.0;                   /* Y(0:k) is now all zero */
272     int top = size;               /* stack for pattern is empty */
273     tags[k] = k;                  /* mark node k as visited */
274     m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
275     int kk = (P) ? (P[k]) : (k);  /* kth original, or permuted, column */
276     int p2 = Ap[kk+1];
277     for (int p = Ap[kk]; p < p2; ++p)
278     {
279       int i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
280       if (i <= k)
281       {
282         y[i] += Ax[p];            /* scatter A(i,k) into Y (sum duplicates) */
283         int len;
284         for (len = 0; tags[i] != k; i = m_parent[i])
285         {
286           pattern[len++] = i;     /* L(k,i) is nonzero */
287           tags[i] = k;            /* mark i as visited */
288         }
289         while (len > 0)
290           pattern[--top] = pattern[--len];
291       }
292     }
293     /* compute numerical values kth row of L (a sparse triangular solve) */
294     m_diag[k] = y[k];                       /* get D(k,k) and clear Y(k) */
295     y[k] = 0.0;
296     for (; top < size; ++top)
297     {
298       int i = pattern[top];      /* pattern[top:n-1] is pattern of L(:,k) */
299       Scalar yi = y[i];          /* get and clear Y(i) */
300       y[i] = 0.0;
301       int p2 = Lp[i] + m_nonZerosPerCol[i];
302       int p;
303       for (p = Lp[i]; p < p2; ++p)
304         y[Li[p]] -= Lx[p] * yi;
305       Scalar l_ki = yi / m_diag[i];       /* the nonzero entry L(k,i) */
306       m_diag[k] -= l_ki * yi;
307       Li[p] = k;                          /* store L(k,i) in column form of L */
308       Lx[p] = l_ki;
309       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
310     }
311     if (m_diag[k] == 0.0)
312     {
313       ok = false;                         /* failure, D(k,k) is zero */
314       break;
315     }
316   }
317
318   ei_aligned_stack_delete(Scalar, y, size);
319   ei_aligned_stack_delete(int, pattern, size);
320   ei_aligned_stack_delete(int, tags, size);
321
322   return ok;  /* success, diagonal of D is all nonzero */
323 }
324
325 /** Computes b = L^-T L^-1 b */
326 template<typename MatrixType, int Backend>
327 template<typename Derived>
328 bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
329 {
330   const int size = m_matrix.rows();
331   ei_assert(size==b.rows());
332   if (!m_succeeded)
333     return false;
334
335   if (m_matrix.nonZeros()>0) // otherwise L==I
336     m_matrix.solveTriangularInPlace(b);
337   b = b.cwise() / m_diag;
338   // FIXME should be .adjoint() but it fails to compile...
339
340   if (m_matrix.nonZeros()>0) // otherwise L==I
341   m_matrix.transpose().solveTriangularInPlace(b);
342
343   return true;
344 }
345
346 #endif // EIGEN_SPARSELDLT_H