Another set of UI messages fixes and tweaks! No functional changes.
[blender.git] / extern / Eigen3 / Eigen / src / Sparse / TriangularSolver.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.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 namespace internal {
29
30 template<typename Lhs, typename Rhs, int Mode,
31   int UpLo = (Mode & Lower)
32            ? Lower
33            : (Mode & Upper)
34            ? Upper
35            : -1,
36   int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit>
37 struct sparse_solve_triangular_selector;
38
39 // forward substitution, row-major
40 template<typename Lhs, typename Rhs, int Mode>
41 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
42 {
43   typedef typename Rhs::Scalar Scalar;
44   static void run(const Lhs& lhs, Rhs& other)
45   {
46     for(int col=0 ; col<other.cols() ; ++col)
47     {
48       for(int i=0; i<lhs.rows(); ++i)
49       {
50         Scalar tmp = other.coeff(i,col);
51         Scalar lastVal = 0;
52         int lastIndex = 0;
53         for(typename Lhs::InnerIterator it(lhs, i); it; ++it)
54         {
55           lastVal = it.value();
56           lastIndex = it.index();
57           if(lastIndex==i)
58             break;
59           tmp -= lastVal * other.coeff(lastIndex,col);
60         }
61         if (Mode & UnitDiag)
62           other.coeffRef(i,col) = tmp;
63         else
64         {
65           eigen_assert(lastIndex==i);
66           other.coeffRef(i,col) = tmp/lastVal;
67         }
68       }
69     }
70   }
71 };
72
73 // backward substitution, row-major
74 template<typename Lhs, typename Rhs, int Mode>
75 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
76 {
77   typedef typename Rhs::Scalar Scalar;
78   static void run(const Lhs& lhs, Rhs& other)
79   {
80     for(int col=0 ; col<other.cols() ; ++col)
81     {
82       for(int i=lhs.rows()-1 ; i>=0 ; --i)
83       {
84         Scalar tmp = other.coeff(i,col);
85         typename Lhs::InnerIterator it(lhs, i);
86         if (it && it.index() == i)
87           ++it;
88         for(; it; ++it)
89         {
90           tmp -= it.value() * other.coeff(it.index(),col);
91         }
92
93         if (Mode & UnitDiag)
94           other.coeffRef(i,col) = tmp;
95         else
96         {
97           typename Lhs::InnerIterator it(lhs, i);
98           eigen_assert(it && it.index() == i);
99           other.coeffRef(i,col) = tmp/it.value();
100         }
101       }
102     }
103   }
104 };
105
106 // forward substitution, col-major
107 template<typename Lhs, typename Rhs, int Mode>
108 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
109 {
110   typedef typename Rhs::Scalar Scalar;
111   static void run(const Lhs& lhs, Rhs& other)
112   {
113     for(int col=0 ; col<other.cols() ; ++col)
114     {
115       for(int i=0; i<lhs.cols(); ++i)
116       {
117         Scalar& tmp = other.coeffRef(i,col);
118         if (tmp!=Scalar(0)) // optimization when other is actually sparse
119         {
120           typename Lhs::InnerIterator it(lhs, i);
121           if(!(Mode & UnitDiag))
122           {
123             eigen_assert(it.index()==i);
124             tmp /= it.value();
125           }
126           if (it && it.index()==i)
127             ++it;
128           for(; it; ++it)
129             other.coeffRef(it.index(), col) -= tmp * it.value();
130         }
131       }
132     }
133   }
134 };
135
136 // backward substitution, col-major
137 template<typename Lhs, typename Rhs, int Mode>
138 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
139 {
140   typedef typename Rhs::Scalar Scalar;
141   static void run(const Lhs& lhs, Rhs& other)
142   {
143     for(int col=0 ; col<other.cols() ; ++col)
144     {
145       for(int i=lhs.cols()-1; i>=0; --i)
146       {
147         Scalar& tmp = other.coeffRef(i,col);
148         if (tmp!=Scalar(0)) // optimization when other is actually sparse
149         {
150           if(!(Mode & UnitDiag))
151           {
152             // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the
153             // last element of the column !
154             other.coeffRef(i,col) /= lhs.innerVector(i).lastCoeff();
155           }
156           typename Lhs::InnerIterator it(lhs, i);
157           for(; it && it.index()<i; ++it)
158             other.coeffRef(it.index(), col) -= tmp * it.value();
159         }
160       }
161     }
162   }
163 };
164
165 } // end namespace internal
166
167 template<typename ExpressionType,int Mode>
168 template<typename OtherDerived>
169 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
170 {
171   eigen_assert(m_matrix.cols() == m_matrix.rows());
172   eigen_assert(m_matrix.cols() == other.rows());
173   eigen_assert(!(Mode & ZeroDiag));
174   eigen_assert(Mode & (Upper|Lower));
175
176   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
177
178   typedef typename internal::conditional<copy,
179     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
180   OtherCopy otherCopy(other.derived());
181
182   internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy);
183
184   if (copy)
185     other = otherCopy;
186 }
187
188 template<typename ExpressionType,int Mode>
189 template<typename OtherDerived>
190 typename internal::plain_matrix_type_column_major<OtherDerived>::type
191 SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const
192 {
193   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
194   solveInPlace(res);
195   return res;
196 }
197
198 // pure sparse path
199
200 namespace internal {
201
202 template<typename Lhs, typename Rhs, int Mode,
203   int UpLo = (Mode & Lower)
204            ? Lower
205            : (Mode & Upper)
206            ? Upper
207            : -1,
208   int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
209 struct sparse_solve_triangular_sparse_selector;
210
211 // forward substitution, col-major
212 template<typename Lhs, typename Rhs, int Mode, int UpLo>
213 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
214 {
215   typedef typename Rhs::Scalar Scalar;
216   typedef typename promote_index_type<typename traits<Lhs>::Index,
217                                          typename traits<Rhs>::Index>::type Index;
218   static void run(const Lhs& lhs, Rhs& other)
219   {
220     const bool IsLower = (UpLo==Lower);
221     AmbiVector<Scalar,Index> tempVector(other.rows()*2);
222     tempVector.setBounds(0,other.rows());
223
224     Rhs res(other.rows(), other.cols());
225     res.reserve(other.nonZeros());
226
227     for(int col=0 ; col<other.cols() ; ++col)
228     {
229       // FIXME estimate number of non zeros
230       tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
231       tempVector.setZero();
232       tempVector.restart();
233       for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
234       {
235         tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
236       }
237
238       for(int i=IsLower?0:lhs.cols()-1;
239           IsLower?i<lhs.cols():i>=0;
240           i+=IsLower?1:-1)
241       {
242         tempVector.restart();
243         Scalar& ci = tempVector.coeffRef(i);
244         if (ci!=Scalar(0))
245         {
246           // find
247           typename Lhs::InnerIterator it(lhs, i);
248           if(!(Mode & UnitDiag))
249           {
250             if (IsLower)
251             {
252               eigen_assert(it.index()==i);
253               ci /= it.value();
254             }
255             else
256               ci /= lhs.coeff(i,i);
257           }
258           tempVector.restart();
259           if (IsLower)
260           {
261             if (it.index()==i)
262               ++it;
263             for(; it; ++it)
264               tempVector.coeffRef(it.index()) -= ci * it.value();
265           }
266           else
267           {
268             for(; it && it.index()<i; ++it)
269               tempVector.coeffRef(it.index()) -= ci * it.value();
270           }
271         }
272       }
273
274
275       int count = 0;
276       // FIXME compute a reference value to filter zeros
277       for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it)
278       {
279         ++ count;
280 //         std::cerr << "fill " << it.index() << ", " << col << "\n";
281 //         std::cout << it.value() << "  ";
282         // FIXME use insertBack
283         res.insert(it.index(), col) = it.value();
284       }
285 //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
286     }
287     res.finalize();
288     other = res.markAsRValue();
289   }
290 };
291
292 } // end namespace internal
293
294 template<typename ExpressionType,int Mode>
295 template<typename OtherDerived>
296 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
297 {
298   eigen_assert(m_matrix.cols() == m_matrix.rows());
299   eigen_assert(m_matrix.cols() == other.rows());
300   eigen_assert(!(Mode & ZeroDiag));
301   eigen_assert(Mode & (Upper|Lower));
302
303 //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
304
305 //   typedef typename internal::conditional<copy,
306 //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
307 //   OtherCopy otherCopy(other.derived());
308
309   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived());
310
311 //   if (copy)
312 //     other = otherCopy;
313 }
314
315 #ifdef EIGEN2_SUPPORT
316
317 // deprecated stuff:
318
319 /** \deprecated */
320 template<typename Derived>
321 template<typename OtherDerived>
322 void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const
323 {
324   this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other);
325 }
326
327 /** \deprecated */
328 template<typename Derived>
329 template<typename OtherDerived>
330 typename internal::plain_matrix_type_column_major<OtherDerived>::type
331 SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
332 {
333   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
334   derived().solveTriangularInPlace(res);
335   return res;
336 }
337 #endif // EIGEN2_SUPPORT
338
339 #endif // EIGEN_SPARSETRIANGULARSOLVER_H