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 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_SPARSETRIANGULARSOLVER_H
11 #define EIGEN_SPARSETRIANGULARSOLVER_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template<typename Lhs, typename Rhs, int Mode,
18   int UpLo = (Mode & Lower)
19            ? Lower
20            : (Mode & Upper)
21            ? Upper
22            : -1,
23   int StorageOrder = int(traits<Lhs>::Flags) & RowMajorBit>
24 struct sparse_solve_triangular_selector;
25 
26 // forward substitution, row-major
27 template<typename Lhs, typename Rhs, int Mode>
28 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,RowMajor>
29 {
30   typedef typename Rhs::Scalar Scalar;
31   static void run(const Lhs& lhs, Rhs& other)
32   {
33     for(int col=0 ; col<other.cols() ; ++col)
34     {
35       for(int i=0; i<lhs.rows(); ++i)
36       {
37         Scalar tmp = other.coeff(i,col);
38         Scalar lastVal(0);
39         int lastIndex = 0;
40         for(typename Lhs::InnerIterator it(lhs, i); it; ++it)
41         {
42           lastVal = it.value();
43           lastIndex = it.index();
44           if(lastIndex==i)
45             break;
46           tmp -= lastVal * other.coeff(lastIndex,col);
47         }
48         if (Mode & UnitDiag)
49           other.coeffRef(i,col) = tmp;
50         else
51         {
52           eigen_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, int Mode>
62 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,RowMajor>
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         Scalar l_ii(0);
73         typename Lhs::InnerIterator it(lhs, i);
74         while(it && it.index()<i)
75           ++it;
76         if(!(Mode & UnitDiag))
77         {
78           eigen_assert(it && it.index()==i);
79           l_ii = it.value();
80           ++it;
81         }
82         else if (it && it.index() == i)
83           ++it;
84         for(; it; ++it)
85         {
86           tmp -= it.value() * other.coeff(it.index(),col);
87         }
88 
89         if (Mode & UnitDiag)
90           other.coeffRef(i,col) = tmp;
91         else
92           other.coeffRef(i,col) = tmp/l_ii;
93       }
94     }
95   }
96 };
97 
98 // forward substitution, col-major
99 template<typename Lhs, typename Rhs, int Mode>
100 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
101 {
102   typedef typename Rhs::Scalar Scalar;
103   static void run(const Lhs& lhs, Rhs& other)
104   {
105     for(int col=0 ; col<other.cols() ; ++col)
106     {
107       for(int i=0; i<lhs.cols(); ++i)
108       {
109         Scalar& tmp = other.coeffRef(i,col);
110         if (tmp!=Scalar(0)) // optimization when other is actually sparse
111         {
112           typename Lhs::InnerIterator it(lhs, i);
113           while(it && it.index()<i)
114             ++it;
115           if(!(Mode & UnitDiag))
116           {
117             eigen_assert(it && it.index()==i);
118             tmp /= it.value();
119           }
120           if (it && it.index()==i)
121             ++it;
122           for(; it; ++it)
123             other.coeffRef(it.index(), col) -= tmp * it.value();
124         }
125       }
126     }
127   }
128 };
129 
130 // backward substitution, col-major
131 template<typename Lhs, typename Rhs, int Mode>
132 struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
133 {
134   typedef typename Rhs::Scalar Scalar;
135   static void run(const Lhs& lhs, Rhs& other)
136   {
137     for(int col=0 ; col<other.cols() ; ++col)
138     {
139       for(int i=lhs.cols()-1; i>=0; --i)
140       {
141         Scalar& tmp = other.coeffRef(i,col);
142         if (tmp!=Scalar(0)) // optimization when other is actually sparse
143         {
144           if(!(Mode & UnitDiag))
145           {
146             // TODO replace this by a binary search. make sure the binary search is safe for partially sorted elements
147             typename Lhs::ReverseInnerIterator it(lhs, i);
148             while(it && it.index()!=i)
149               --it;
150             eigen_assert(it && it.index()==i);
151             other.coeffRef(i,col) /= it.value();
152           }
153           typename Lhs::InnerIterator it(lhs, i);
154           for(; it && it.index()<i; ++it)
155             other.coeffRef(it.index(), col) -= tmp * it.value();
156         }
157       }
158     }
159   }
160 };
161 
162 } // end namespace internal
163 
164 template<typename ExpressionType,int Mode>
165 template<typename OtherDerived>
166 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDerived>& other) const
167 {
168   eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
169   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
170 
171   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
172 
173   typedef typename internal::conditional<copy,
174     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
175   OtherCopy otherCopy(other.derived());
176 
177   internal::sparse_solve_triangular_selector<ExpressionType, typename internal::remove_reference<OtherCopy>::type, Mode>::run(m_matrix, otherCopy);
178 
179   if (copy)
180     other = otherCopy;
181 }
182 
183 template<typename ExpressionType,int Mode>
184 template<typename OtherDerived>
185 typename internal::plain_matrix_type_column_major<OtherDerived>::type
186 SparseTriangularView<ExpressionType,Mode>::solve(const MatrixBase<OtherDerived>& other) const
187 {
188   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
189   solveInPlace(res);
190   return res;
191 }
192 
193 // pure sparse path
194 
195 namespace internal {
196 
197 template<typename Lhs, typename Rhs, int Mode,
198   int UpLo = (Mode & Lower)
199            ? Lower
200            : (Mode & Upper)
201            ? Upper
202            : -1,
203   int StorageOrder = int(Lhs::Flags) & (RowMajorBit)>
204 struct sparse_solve_triangular_sparse_selector;
205 
206 // forward substitution, col-major
207 template<typename Lhs, typename Rhs, int Mode, int UpLo>
208 struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
209 {
210   typedef typename Rhs::Scalar Scalar;
211   typedef typename promote_index_type<typename traits<Lhs>::Index,
212                                          typename traits<Rhs>::Index>::type Index;
213   static void run(const Lhs& lhs, Rhs& other)
214   {
215     const bool IsLower = (UpLo==Lower);
216     AmbiVector<Scalar,Index> tempVector(other.rows()*2);
217     tempVector.setBounds(0,other.rows());
218 
219     Rhs res(other.rows(), other.cols());
220     res.reserve(other.nonZeros());
221 
222     for(int col=0 ; col<other.cols() ; ++col)
223     {
224       // FIXME estimate number of non zeros
225       tempVector.init(.99/*float(other.col(col).nonZeros())/float(other.rows())*/);
226       tempVector.setZero();
227       tempVector.restart();
228       for (typename Rhs::InnerIterator rhsIt(other, col); rhsIt; ++rhsIt)
229       {
230         tempVector.coeffRef(rhsIt.index()) = rhsIt.value();
231       }
232 
233       for(int i=IsLower?0:lhs.cols()-1;
234           IsLower?i<lhs.cols():i>=0;
235           i+=IsLower?1:-1)
236       {
237         tempVector.restart();
238         Scalar& ci = tempVector.coeffRef(i);
239         if (ci!=Scalar(0))
240         {
241           // find
242           typename Lhs::InnerIterator it(lhs, i);
243           if(!(Mode & UnitDiag))
244           {
245             if (IsLower)
246             {
247               eigen_assert(it.index()==i);
248               ci /= it.value();
249             }
250             else
251               ci /= lhs.coeff(i,i);
252           }
253           tempVector.restart();
254           if (IsLower)
255           {
256             if (it.index()==i)
257               ++it;
258             for(; it; ++it)
259               tempVector.coeffRef(it.index()) -= ci * it.value();
260           }
261           else
262           {
263             for(; it && it.index()<i; ++it)
264               tempVector.coeffRef(it.index()) -= ci * it.value();
265           }
266         }
267       }
268 
269 
270       int count = 0;
271       // FIXME compute a reference value to filter zeros
272       for (typename AmbiVector<Scalar,Index>::Iterator it(tempVector/*,1e-12*/); it; ++it)
273       {
274         ++ count;
275 //         std::cerr << "fill " << it.index() << ", " << col << "\n";
276 //         std::cout << it.value() << "  ";
277         // FIXME use insertBack
278         res.insert(it.index(), col) = it.value();
279       }
280 //       std::cout << "tempVector.nonZeros() == " << int(count) << " / " << (other.rows()) << "\n";
281     }
282     res.finalize();
283     other = res.markAsRValue();
284   }
285 };
286 
287 } // end namespace internal
288 
289 template<typename ExpressionType,int Mode>
290 template<typename OtherDerived>
291 void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<OtherDerived>& other) const
292 {
293   eigen_assert(m_matrix.cols() == m_matrix.rows() && m_matrix.cols() == other.rows());
294   eigen_assert( (!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
295 
296 //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
297 
298 //   typedef typename internal::conditional<copy,
299 //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
300 //   OtherCopy otherCopy(other.derived());
301 
302   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(m_matrix, other.derived());
303 
304 //   if (copy)
305 //     other = otherCopy;
306 }
307 
308 #ifdef EIGEN2_SUPPORT
309 
310 // deprecated stuff:
311 
312 /** \deprecated */
313 template<typename Derived>
314 template<typename OtherDerived>
315 void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const
316 {
317   this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other);
318 }
319 
320 /** \deprecated */
321 template<typename Derived>
322 template<typename OtherDerived>
323 typename internal::plain_matrix_type_column_major<OtherDerived>::type
324 SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
325 {
326   typename internal::plain_matrix_type_column_major<OtherDerived>::type res(other);
327   derived().solveTriangularInPlace(res);
328   return res;
329 }
330 #endif // EIGEN2_SUPPORT
331 
332 } // end namespace Eigen
333 
334 #endif // EIGEN_SPARSETRIANGULARSOLVER_H
335