Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / SparseLU.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_SPARSELU_H
26 #define EIGEN_SPARSELU_H
27
28 /** \ingroup Sparse_Module
29   *
30   * \class SparseLU
31   *
32   * \brief LU decomposition of a sparse matrix and associated features
33   *
34   * \param MatrixType the type of the matrix of which we are computing the LU factorization
35   *
36   * \sa class LU, class SparseLLT
37   */
38 template<typename MatrixType, int Backend = DefaultBackend>
39 class SparseLU
40 {
41   protected:
42     typedef typename MatrixType::Scalar Scalar;
43     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
44     typedef SparseMatrix<Scalar,LowerTriangular> LUMatrixType;
45
46     enum {
47       MatrixLUIsDirty             = 0x10000
48     };
49
50   public:
51
52     /** Creates a dummy LU factorization object with flags \a flags. */
53     SparseLU(int flags = 0)
54       : m_flags(flags), m_status(0)
55     {
56       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
57     }
58
59     /** Creates a LU object and compute the respective factorization of \a matrix using
60       * flags \a flags. */
61     SparseLU(const MatrixType& matrix, int flags = 0)
62       : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0)
63     {
64       m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
65       compute(matrix);
66     }
67
68     /** Sets the relative threshold value used to prune zero coefficients during the decomposition.
69       *
70       * Setting a value greater than zero speeds up computation, and yields to an imcomplete
71       * factorization with fewer non zero coefficients. Such approximate factors are especially
72       * useful to initialize an iterative solver.
73       *
74       * Note that the exact meaning of this parameter might depends on the actual
75       * backend. Moreover, not all backends support this feature.
76       *
77       * \sa precision() */
78     void setPrecision(RealScalar v) { m_precision = v; }
79
80     /** \returns the current precision.
81       *
82       * \sa setPrecision() */
83     RealScalar precision() const { return m_precision; }
84
85     /** Sets the flags. Possible values are:
86       *  - CompleteFactorization
87       *  - IncompleteFactorization
88       *  - MemoryEfficient
89       *  - one of the ordering methods
90       *  - etc...
91       *
92       * \sa flags() */
93     void setFlags(int f) { m_flags = f; }
94     /** \returns the current flags */
95     int flags() const { return m_flags; }
96
97     void setOrderingMethod(int m)
98     {
99       ei_assert(m&~OrderingMask == 0 && m!=0 && "invalid ordering method");
100       m_flags = m_flags&~OrderingMask | m&OrderingMask;
101     }
102
103     int orderingMethod() const
104     {
105       return m_flags&OrderingMask;
106     }
107
108     /** Computes/re-computes the LU factorization */
109     void compute(const MatrixType& matrix);
110
111     /** \returns the lower triangular matrix L */
112     //inline const MatrixType& matrixL() const { return m_matrixL; }
113
114     /** \returns the upper triangular matrix U */
115     //inline const MatrixType& matrixU() const { return m_matrixU; }
116
117     template<typename BDerived, typename XDerived>
118     bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
119
120     /** \returns true if the factorization succeeded */
121     inline bool succeeded(void) const { return m_succeeded; }
122
123   protected:
124     RealScalar m_precision;
125     int m_flags;
126     mutable int m_status;
127     bool m_succeeded;
128 };
129
130 /** Computes / recomputes the LU decomposition of matrix \a a
131   * using the default algorithm.
132   */
133 template<typename MatrixType, int Backend>
134 void SparseLU<MatrixType,Backend>::compute(const MatrixType& a)
135 {
136   ei_assert(false && "not implemented yet");
137 }
138
139 /** Computes *x = U^-1 L^-1 b */
140 template<typename MatrixType, int Backend>
141 template<typename BDerived, typename XDerived>
142 bool SparseLU<MatrixType,Backend>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const
143 {
144   ei_assert(false && "not implemented yet");
145   return false;
146 }
147
148 #endif // EIGEN_SPARSELU_H