SparseCore/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 // 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>
25 
26 // forward substitution, row-major
27 template<typename Lhs, typename Rhs, int Mode>
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>
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>
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 
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 
178 
179  if (copy)
180  other = otherCopy;
181 }
182 
183 template<typename ExpressionType,int Mode>
184 template<typename OtherDerived>
187 {
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)>
205 
206 // forward substitution, col-major
207 template<typename Lhs, typename Rhs, int Mode, int UpLo>
209 {
210  typedef typename Rhs::Scalar Scalar;
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>
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 
303 
304 // if (copy)
305 // other = otherCopy;
306 }
307 
308 #ifdef EIGEN2_SUPPORT
309 
310 // deprecated stuff:
311 
313 template<typename Derived>
314 template<typename OtherDerived>
316 {
317  this->template triangular<Flags&(Upper|Lower)>().solveInPlace(other);
318 }
319 
321 template<typename Derived>
322 template<typename OtherDerived>
325 {
327  derived().solveTriangularInPlace(res);
328  return res;
329 }
330 #endif // EIGEN2_SUPPORT
331 
332 } // end namespace Eigen
333 
334 #endif // EIGEN_SPARSETRIANGULARSOLVER_H
void setBounds(Index start, Index end)
Definition: AmbiVector.h:42
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
const unsigned int RowMajorBit
internal::plain_matrix_type_column_major< OtherDerived >::type solve(const MatrixBase< OtherDerived > &other) const
Base class of any sparse matrices or sparse expressions.
const Derived & derived() const
void solveInPlace(MatrixBase< OtherDerived > &other) const
promote_index_type< typename traits< Lhs >::Index, typename traits< Rhs >::Index >::type Index
ColXpr col(Index i)
Definition: BlockMethods.h:708
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
#define eigen_assert(x)
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:14