Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / TaucsSupport.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_TAUCSSUPPORT_H
26 #define EIGEN_TAUCSSUPPORT_H
27
28 template<typename Derived>
29 taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix()
30 {
31   taucs_ccs_matrix res;
32   res.n         = cols();
33   res.m         = rows();
34   res.flags     = 0;
35   res.colptr    = derived()._outerIndexPtr();
36   res.rowind    = derived()._innerIndexPtr();
37   res.values.v  = derived()._valuePtr();
38   if (ei_is_same_type<Scalar,int>::ret)
39     res.flags |= TAUCS_INT;
40   else if (ei_is_same_type<Scalar,float>::ret)
41     res.flags |= TAUCS_SINGLE;
42   else if (ei_is_same_type<Scalar,double>::ret)
43     res.flags |= TAUCS_DOUBLE;
44   else if (ei_is_same_type<Scalar,std::complex<float> >::ret)
45     res.flags |= TAUCS_SCOMPLEX;
46   else if (ei_is_same_type<Scalar,std::complex<double> >::ret)
47     res.flags |= TAUCS_DCOMPLEX;
48   else
49   {
50     ei_assert(false && "Scalar type not supported by TAUCS");
51   }
52
53   if (Flags & UpperTriangular)
54     res.flags |= TAUCS_UPPER;
55   if (Flags & LowerTriangular)
56     res.flags |= TAUCS_LOWER;
57   if (Flags & SelfAdjoint)
58     res.flags |= (NumTraits<Scalar>::IsComplex ? TAUCS_HERMITIAN : TAUCS_SYMMETRIC);
59   else if ((Flags & UpperTriangular) || (Flags & LowerTriangular))
60     res.flags |= TAUCS_TRIANGULAR;
61
62   return res;
63 }
64
65 template<typename Scalar, int Flags>
66 MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(taucs_ccs_matrix& taucsMat)
67 {
68   m_innerSize = taucsMat.m;
69   m_outerSize = taucsMat.n;
70   m_outerIndex = taucsMat.colptr;
71   m_innerIndices = taucsMat.rowind;
72   m_values = reinterpret_cast<Scalar*>(taucsMat.values.v);
73   m_nnz = taucsMat.colptr[taucsMat.n];
74 }
75
76 template<typename MatrixType>
77 class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
78 {
79   protected:
80     typedef SparseLLT<MatrixType> Base;
81     typedef typename Base::Scalar Scalar;
82     typedef typename Base::RealScalar RealScalar;
83     using Base::MatrixLIsDirty;
84     using Base::SupernodalFactorIsDirty;
85     using Base::m_flags;
86     using Base::m_matrix;
87     using Base::m_status;
88
89   public:
90
91     SparseLLT(int flags = 0)
92       : Base(flags), m_taucsSupernodalFactor(0)
93     {
94     }
95
96     SparseLLT(const MatrixType& matrix, int flags = 0)
97       : Base(flags), m_taucsSupernodalFactor(0)
98     {
99       compute(matrix);
100     }
101
102     ~SparseLLT()
103     {
104       if (m_taucsSupernodalFactor)
105         taucs_supernodal_factor_free(m_taucsSupernodalFactor);
106     }
107
108     inline const typename Base::CholMatrixType& matrixL(void) const;
109
110     template<typename Derived>
111     void solveInPlace(MatrixBase<Derived> &b) const;
112
113     void compute(const MatrixType& matrix);
114
115   protected:
116     void* m_taucsSupernodalFactor;
117 };
118
119 template<typename MatrixType>
120 void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
121 {
122   if (m_taucsSupernodalFactor)
123   {
124     taucs_supernodal_factor_free(m_taucsSupernodalFactor);
125     m_taucsSupernodalFactor = 0;
126   }
127
128   if (m_flags & IncompleteFactorization)
129   {
130     taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
131     taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
132     // the matrix returned by Taucs is not necessarily sorted,
133     // so let's copy it in two steps
134     DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsRes);
135     m_matrix = tmp;
136     free(taucsRes);
137     m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
138              | IncompleteFactorization
139              | SupernodalFactorIsDirty;
140   }
141   else
142   {
143     taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
144     if ( (m_flags & SupernodalLeftLooking)
145       || ((!(m_flags & SupernodalMultifrontal)) && (m_flags & MemoryEfficient)) )
146     {
147       m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA);
148     }
149     else
150     {
151       // use the faster Multifrontal routine
152       m_taucsSupernodalFactor = taucs_ccs_factor_llt_mf(&taucsMatA);
153     }
154     m_status = (m_status & ~IncompleteFactorization) | CompleteFactorization | MatrixLIsDirty;
155   }
156 }
157
158 template<typename MatrixType>
159 inline const typename SparseLLT<MatrixType>::CholMatrixType&
160 SparseLLT<MatrixType,Taucs>::matrixL() const
161 {
162   if (m_status & MatrixLIsDirty)
163   {
164     ei_assert(!(m_status & SupernodalFactorIsDirty));
165
166     taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor);
167
168     // the matrix returned by Taucs is not necessarily sorted,
169     // so let's copy it in two steps
170     DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsL);
171     const_cast<typename Base::CholMatrixType&>(m_matrix) = tmp;
172     free(taucsL);
173     m_status = (m_status & ~MatrixLIsDirty);
174   }
175   return m_matrix;
176 }
177
178 template<typename MatrixType>
179 template<typename Derived>
180 void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
181 {
182   bool inputIsCompatibleWithTaucs = (Derived::Flags&RowMajorBit)==0;
183
184   if (!inputIsCompatibleWithTaucs)
185   {
186     matrixL();
187     Base::solveInPlace(b);
188   }
189   else if (m_flags & IncompleteFactorization)
190   {
191     taucs_ccs_matrix taucsLLT = const_cast<typename Base::CholMatrixType&>(m_matrix).asTaucsMatrix();
192     typename ei_plain_matrix_type<Derived>::type x(b.rows());
193     for (int j=0; j<b.cols(); ++j)
194     {
195       taucs_ccs_solve_llt(&taucsLLT,x.data(),&b.col(j).coeffRef(0));
196       b.col(j) = x;
197     }
198   }
199   else
200   {
201     typename ei_plain_matrix_type<Derived>::type x(b.rows());
202     for (int j=0; j<b.cols(); ++j)
203     {
204       taucs_supernodal_solve_llt(m_taucsSupernodalFactor,x.data(),&b.col(j).coeffRef(0));
205       b.col(j) = x;
206     }
207   }
208 }
209
210 #endif // EIGEN_TAUCSSUPPORT_H