Merge of itasc branch. Project files, scons and cmake should be working. Makefile...
[blender.git] / extern / Eigen2 / Eigen / src / Sparse / TriangularSolver.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_SPARSETRIANGULARSOLVER_H
26 #define EIGEN_SPARSETRIANGULARSOLVER_H
27
28 // forward substitution, row-major
29 template<typename Lhs, typename Rhs>
30 struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse>
31 {
32   typedef typename Rhs::Scalar Scalar;
33   static void run(const Lhs& lhs, Rhs& other)
34   {
35     for(int col=0 ; col<other.cols() ; ++col)
36     {
37       for(int i=0; i<lhs.rows(); ++i)
38       {
39         Scalar tmp = other.coeff(i,col);
40         Scalar lastVal = 0;
41         int lastIndex = 0;
42         for(typename Lhs::InnerIterator it(lhs, i); it; ++it)
43         {
44           lastVal = it.value();
45           lastIndex = it.index();
46           tmp -= lastVal * other.coeff(lastIndex,col);
47         }
48         if (Lhs::Flags & UnitDiagBit)
49           other.coeffRef(i,col) = tmp;
50         else
51         {
52           ei_assert(lastIndex==i);
53           other.coeffRef(i,col) = tmp/lastVal;
54         }
55       }
56     }
57   }
58 };
59
60 // backward substitution, row-major
61 template<typename Lhs, typename Rhs>
62 struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,RowMajor|IsSparse>
63 {
64   typedef typename Rhs::Scalar Scalar;
65   static void run(const Lhs& lhs, Rhs& other)
66   {
67     for(int col=0 ; col<other.cols() ; ++col)
68     {
69       for(int i=lhs.rows()-1 ; i>=0 ; --i)
70       {
71         Scalar tmp = other.coeff(i,col);
72         typename Lhs::InnerIterator it(lhs, i);
73         if (it.index() == i)
74           ++it;
75         for(; it; ++it)
76         {
77           tmp -= it.value() * other.coeff(it.index(),col);
78         }
79
80         if (Lhs::Flags & UnitDiagBit)
81           other.coeffRef(i,col) = tmp;
82         else
83         {
84           typename Lhs::InnerIterator it(lhs, i);
85           ei_assert(it.index() == i);
86           other.coeffRef(i,col) = tmp/it.value();
87         }
88       }
89     }
90   }
91 };
92
93 // forward substitution, col-major
94 template<typename Lhs, typename Rhs>
95 struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,ColMajor|IsSparse>
96 {
97   typedef typename Rhs::Scalar Scalar;
98   static void run(const Lhs& lhs, Rhs& other)
99   {
100     for(int col=0 ; col<other.cols() ; ++col)
101     {
102       for(int i=0; i<lhs.cols(); ++i)
103       {
104         typename Lhs::InnerIterator it(lhs, i);
105         if(!(Lhs::Flags & UnitDiagBit))
106         {
107           // std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
108           ei_assert(it.index()==i);
109           other.coeffRef(i,col) /= it.value();
110         }
111         Scalar tmp = other.coeffRef(i,col);
112         if (it.index()==i)
113           ++it;
114         for(; it; ++it)
115           other.coeffRef(it.index(), col) -= tmp * it.value();
116       }
117     }
118   }
119 };
120
121 // backward substitution, col-major
122 template<typename Lhs, typename Rhs>
123 struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,ColMajor|IsSparse>
124 {
125   typedef typename Rhs::Scalar Scalar;
126   static void run(const Lhs& lhs, Rhs& other)
127   {
128     for(int col=0 ; col<other.cols() ; ++col)
129     {
130       for(int i=lhs.cols()-1; i>=0; --i)
131       {
132         if(!(Lhs::Flags & UnitDiagBit))
133         {
134           // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the
135           // last element of the column !
136           other.coeffRef(i,col) /= lhs.coeff(i,i);
137         }
138         Scalar tmp = other.coeffRef(i,col);
139         typename Lhs::InnerIterator it(lhs, i);
140         for(; it && it.index()<i; ++it)
141           other.coeffRef(it.index(), col) -= tmp * it.value();
142       }
143     }
144   }
145 };
146
147 template<typename Derived>
148 template<typename OtherDerived>
149 void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const
150 {
151   ei_assert(derived().cols() == derived().rows());
152   ei_assert(derived().cols() == other.rows());
153   ei_assert(!(Flags & ZeroDiagBit));
154   ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit));
155
156   enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit };
157
158   typedef typename ei_meta_if<copy,
159     typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy;
160   OtherCopy otherCopy(other.derived());
161
162   ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy);
163
164   if (copy)
165     other = otherCopy;
166 }
167
168 template<typename Derived>
169 template<typename OtherDerived>
170 typename ei_plain_matrix_type_column_major<OtherDerived>::type
171 SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
172 {
173   typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other);
174   solveTriangularInPlace(res);
175   return res;
176 }
177
178 #endif // EIGEN_SPARSETRIANGULARSOLVER_H