Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseLLT.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 #ifndef EIGEN_SPARSELLT_H
26 #define EIGEN_SPARSELLT_H
27
28 /** \ingroup Sparse_Module
29   *
30   * \class SparseLLT
31   *
32   * \brief LLT Cholesky decomposition of a sparse matrix and associated features
33   *
34   * \param MatrixType the type of the matrix of which we are computing the LLT Cholesky decomposition
35   *
36   * \sa class LLT, class LDLT
37   */
38 template<typename MatrixType, int Backend = DefaultBackend>
39 class SparseLLT
40 {
41   protected:
42     typedef typename MatrixType::Scalar Scalar;
43     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
44     typedef SparseMatrix<Scalar,LowerTriangular> CholMatrixType;
45
46     enum {
47       SupernodalFactorIsDirty      = 0x10000,
48       MatrixLIsDirty               = 0x20000
49     };
50
51   public:
52
53     /** Creates a dummy LLT factorization object with flags \a flags. */
54     SparseLLT(int flags = 0)
55       : m_flags(flags), m_status(0)
56     {
57       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
58     }
59
60     /** Creates a LLT object and compute the respective factorization of \a matrix using
61       * flags \a flags. */
62     SparseLLT(const MatrixType& matrix, int flags = 0)
63       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
64     {
65       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
66       compute(matrix);
67     }
68
69     /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
70       *
71       * Setting a value greater than zero speeds up computation, and yields to an imcomplete
72       * factorization with fewer non zero coefficients. Such approximate factors are especially
73       * useful to initialize an iterative solver.
74       *
75       * \warning if precision is greater that zero, the LLT factorization is not guaranteed to succeed
76       * even if the matrix is positive definite.
77       *
78       * Note that the exact meaning of this parameter might depends on the actual
79       * backend. Moreover, not all backends support this feature.
80       *
81       * \sa precision() */
82     void setPrecision(RealScalar v) { m_precision = v; }
83
84     /** \returns the current precision.
85       *
86       * \sa setPrecision() */
87     RealScalar precision() const { return m_precision; }
88
89     /** Sets the flags. Possible values are:
90       *  - CompleteFactorization
91       *  - IncompleteFactorization
92       *  - MemoryEfficient          (hint to use the memory most efficient method offered by the backend)
93       *  - SupernodalMultifrontal   (implies a complete factorization if supported by the backend,
94       *                              overloads the MemoryEfficient flags)
95       *  - SupernodalLeftLooking    (implies a complete factorization  if supported by the backend,
96       *                              overloads the MemoryEfficient flags)
97       *
98       * \sa flags() */
99     void setFlags(int f) { m_flags = f; }
100     /** \returns the current flags */
101     int flags() const { return m_flags; }
102
103     /** Computes/re-computes the LLT factorization */
104     void compute(const MatrixType& matrix);
105
106     /** \returns the lower triangular matrix L */
107     inline const CholMatrixType& matrixL(void) const { return m_matrix; }
108
109     template<typename Derived>
110     bool solveInPlace(MatrixBase<Derived> &b) const;
111
112     /** \returns true if the factorization succeeded */
113     inline bool succeeded(void) const { return m_succeeded; }
114
115   protected:
116     CholMatrixType m_matrix;
117     RealScalar m_precision;
118     int m_flags;
119     mutable int m_status;
120     bool m_succeeded;
121 };
122
123 /** Computes / recomputes the LLT decomposition of matrix \a a
124   * using the default algorithm.
125   */
126 template<typename MatrixType, int Backend>
127 void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
128 {
129   assert(a.rows()==a.cols());
130   const int size = a.rows();
131   m_matrix.resize(size, size);
132
133   // allocate a temporary vector for accumulations
134   AmbiVector<Scalar> tempVector(size);
135   RealScalar density = a.nonZeros()/RealScalar(size*size);
136
137   // TODO estimate the number of non zeros
138   m_matrix.startFill(a.nonZeros()*2);
139   for (int j = 0; j < size; ++j)
140   {
141     Scalar x = ei_real(a.coeff(j,j));
142
143     // TODO better estimate of the density !
144     tempVector.init(density>0.001? IsDense : IsSparse);
145     tempVector.setBounds(j+1,size);
146     tempVector.setZero();
147     // init with current matrix a
148     {
149       typename MatrixType::InnerIterator it(a,j);
150       ++it; // skip diagonal element
151       for (; it; ++it)
152         tempVector.coeffRef(it.index()) = it.value();
153     }
154     for (int k=0; k<j+1; ++k)
155     {
156       typename CholMatrixType::InnerIterator it(m_matrix, k);
157       while (it && it.index()<j)
158         ++it;
159       if (it && it.index()==j)
160       {
161         Scalar y = it.value();
162         x -= ei_abs2(y);
163         ++it; // skip j-th element, and process remaining column coefficients
164         tempVector.restart();
165         for (; it; ++it)
166         {
167           tempVector.coeffRef(it.index()) -= it.value() * y;
168         }
169       }
170     }
171     // copy the temporary vector to the respective m_matrix.col()
172     // while scaling the result by 1/real(x)
173     RealScalar rx = ei_sqrt(ei_real(x));
174     m_matrix.fill(j,j) = rx;
175     Scalar y = Scalar(1)/rx;
176     for (typename AmbiVector<Scalar>::Iterator it(tempVector, m_precision*rx); it; ++it)
177     {
178       m_matrix.fill(it.index(), j) = it.value() * y;
179     }
180   }
181   m_matrix.endFill();
182 }
183
184 /** Computes b = L^-T L^-1 b */
185 template<typename MatrixType, int Backend>
186 template<typename Derived>
187 bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
188 {
189   const int size = m_matrix.rows();
190   ei_assert(size==b.rows());
191
192   m_matrix.solveTriangularInPlace(b);
193   // FIXME should be simply .adjoint() but it fails to compile...
194   if (NumTraits<Scalar>::IsComplex)
195   {
196     CholMatrixType aux = m_matrix.conjugate();
197     aux.transpose().solveTriangularInPlace(b);
198   }
199   else
200     m_matrix.transpose().solveTriangularInPlace(b);
201
202   return true;
203 }
204
205 #endif // EIGEN_SPARSELLT_H