Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / CholmodSupport.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-2009 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_CHOLMODSUPPORT_H
26 #define EIGEN_CHOLMODSUPPORT_H
27
28 template<typename Scalar, typename CholmodType>
29 void ei_cholmod_configure_matrix(CholmodType& mat)
30 {
31   if (ei_is_same_type<Scalar,float>::ret)
32   {
33     mat.xtype = CHOLMOD_REAL;
34     mat.dtype = 1;
35   }
36   else if (ei_is_same_type<Scalar,double>::ret)
37   {
38     mat.xtype = CHOLMOD_REAL;
39     mat.dtype = 0;
40   }
41   else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
42   {
43     mat.xtype = CHOLMOD_COMPLEX;
44     mat.dtype = 1;
45   }
46   else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
47   {
48     mat.xtype = CHOLMOD_COMPLEX;
49     mat.dtype = 0;
50   }
51   else
52   {
53     ei_assert(false && "Scalar type not supported by CHOLMOD");
54   }
55 }
56
57 template<typename Derived>
58 cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
59 {
60   typedef typename Derived::Scalar Scalar;
61   cholmod_sparse res;
62   res.nzmax   = nonZeros();
63   res.nrow    = rows();;
64   res.ncol    = cols();
65   res.p       = derived()._outerIndexPtr();
66   res.i       = derived()._innerIndexPtr();
67   res.x       = derived()._valuePtr();
68   res.xtype = CHOLMOD_REAL;
69   res.itype = CHOLMOD_INT;
70   res.sorted = 1;
71   res.packed = 1;
72   res.dtype = 0;
73   res.stype = -1;
74
75   ei_cholmod_configure_matrix<Scalar>(res);
76
77   if (Derived::Flags & SelfAdjoint)
78   {
79     if (Derived::Flags & UpperTriangular)
80       res.stype = 1;
81     else if (Derived::Flags & LowerTriangular)
82       res.stype = -1;
83     else
84       res.stype = 0;
85   }
86   else
87     res.stype = 0;
88
89   return res;
90 }
91
92 template<typename Derived>
93 cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
94 {
95   EIGEN_STATIC_ASSERT((ei_traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
96   typedef typename Derived::Scalar Scalar;
97
98   cholmod_dense res;
99   res.nrow   = mat.rows();
100   res.ncol   = mat.cols();
101   res.nzmax  = res.nrow * res.ncol;
102   res.d      = mat.derived().stride();
103   res.x      = mat.derived().data();
104   res.z      = 0;
105
106   ei_cholmod_configure_matrix<Scalar>(res);
107
108   return res;
109 }
110
111 template<typename Scalar, int Flags>
112 MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(cholmod_sparse& cm)
113 {
114   m_innerSize = cm.nrow;
115   m_outerSize = cm.ncol;
116   m_outerIndex = reinterpret_cast<int*>(cm.p);
117   m_innerIndices = reinterpret_cast<int*>(cm.i);
118   m_values = reinterpret_cast<Scalar*>(cm.x);
119   m_nnz = m_outerIndex[cm.ncol];
120 }
121
122 template<typename MatrixType>
123 class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
124 {
125   protected:
126     typedef SparseLLT<MatrixType> Base;
127     typedef typename Base::Scalar Scalar;
128     typedef typename Base::RealScalar RealScalar;
129     using Base::MatrixLIsDirty;
130     using Base::SupernodalFactorIsDirty;
131     using Base::m_flags;
132     using Base::m_matrix;
133     using Base::m_status;
134
135   public:
136
137     SparseLLT(int flags = 0)
138       : Base(flags), m_cholmodFactor(0)
139     {
140       cholmod_start(&m_cholmod);
141     }
142
143     SparseLLT(const MatrixType& matrix, int flags = 0)
144       : Base(flags), m_cholmodFactor(0)
145     {
146       cholmod_start(&m_cholmod);
147       compute(matrix);
148     }
149
150     ~SparseLLT()
151     {
152       if (m_cholmodFactor)
153         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
154       cholmod_finish(&m_cholmod);
155     }
156
157     inline const typename Base::CholMatrixType& matrixL(void) const;
158
159     template<typename Derived>
160     void solveInPlace(MatrixBase<Derived> &b) const;
161
162     void compute(const MatrixType& matrix);
163
164   protected:
165     mutable cholmod_common m_cholmod;
166     cholmod_factor* m_cholmodFactor;
167 };
168
169 template<typename MatrixType>
170 void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
171 {
172   if (m_cholmodFactor)
173   {
174     cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
175     m_cholmodFactor = 0;
176   }
177
178   cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
179   m_cholmod.supernodal = CHOLMOD_AUTO;
180   // TODO
181   if (m_flags&IncompleteFactorization)
182   {
183     m_cholmod.nmethods = 1;
184     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
185     m_cholmod.postorder = 0;
186   }
187   else
188   {
189     m_cholmod.nmethods = 1;
190     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
191     m_cholmod.postorder = 0;
192   }
193   m_cholmod.final_ll = 1;
194   m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
195   cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
196
197   m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
198 }
199
200 template<typename MatrixType>
201 inline const typename SparseLLT<MatrixType>::CholMatrixType&
202 SparseLLT<MatrixType,Cholmod>::matrixL() const
203 {
204   if (m_status & MatrixLIsDirty)
205   {
206     ei_assert(!(m_status & SupernodalFactorIsDirty));
207
208     cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
209     const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
210     free(cmRes);
211
212     m_status = (m_status & ~MatrixLIsDirty);
213   }
214   return m_matrix;
215 }
216
217 template<typename MatrixType>
218 template<typename Derived>
219 void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
220 {
221   const int size = m_cholmodFactor->n;
222   ei_assert(size==b.rows());
223
224   // this uses Eigen's triangular sparse solver
225 //   if (m_status & MatrixLIsDirty)
226 //     matrixL();
227 //   Base::solveInPlace(b);
228   // as long as our own triangular sparse solver is not fully optimal,
229   // let's use CHOLMOD's one:
230   cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b);
231   cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
232   b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
233   cholmod_free_dense(&x, &m_cholmod);
234 }
235
236 #endif // EIGEN_CHOLMODSUPPORT_H