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