Eigen  3.2.91
SparseTriangularView.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11 #ifndef EIGEN_SPARSE_TRIANGULARVIEW_H
12 #define EIGEN_SPARSE_TRIANGULARVIEW_H
13 
14 namespace Eigen {
15 
25 template<typename MatrixType, unsigned int Mode> class TriangularViewImpl<MatrixType,Mode,Sparse>
26  : public SparseMatrixBase<TriangularView<MatrixType,Mode> >
27 {
28  enum { SkipFirst = ((Mode&Lower) && !(MatrixType::Flags&RowMajorBit))
29  || ((Mode&Upper) && (MatrixType::Flags&RowMajorBit)),
30  SkipLast = !SkipFirst,
31  SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
32  HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
33  };
34 
36 
37 protected:
38  // dummy solve function to make TriangularView happy.
39  void solve() const;
40 
41  public:
42 
43  EIGEN_SPARSE_PUBLIC_INTERFACE(TriangularViewType)
44 
45  class InnerIterator;
46  class ReverseInnerIterator;
47 
48  typedef typename MatrixType::Nested MatrixTypeNested;
49  typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
50  typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
51 
52  template<typename RhsType, typename DstType>
53  EIGEN_DEVICE_FUNC
54  EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const {
55  if(!(internal::is_same<RhsType,DstType>::value && internal::extract_data(dst) == internal::extract_data(rhs)))
56  dst = rhs;
57  this->solveInPlace(dst);
58  }
59 
60  template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
61  template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
62 
63 };
64 
65 template<typename MatrixType, unsigned int Mode>
66 class TriangularViewImpl<MatrixType,Mode,Sparse>::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator
67 {
68  typedef typename MatrixTypeNestedCleaned::InnerIterator Base;
69  public:
70 
71  EIGEN_STRONG_INLINE InnerIterator(const TriangularViewImpl& view, Index outer)
72  : Base(view.derived().nestedExpression(), outer), m_returnOne(false)
73  {
74  if(SkipFirst)
75  {
76  while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()<outer))
77  Base::operator++();
78  if(HasUnitDiag)
79  m_returnOne = true;
80  }
81  else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
82  {
83  if((!SkipFirst) && Base::operator bool())
84  Base::operator++();
85  m_returnOne = true;
86  }
87  }
88 
89  EIGEN_STRONG_INLINE InnerIterator& operator++()
90  {
91  if(HasUnitDiag && m_returnOne)
92  m_returnOne = false;
93  else
94  {
95  Base::operator++();
96  if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer()))
97  {
98  if((!SkipFirst) && Base::operator bool())
99  Base::operator++();
100  m_returnOne = true;
101  }
102  }
103  return *this;
104  }
105 
106  inline Index row() const { return (MatrixType::Flags&RowMajorBit ? Base::outer() : this->index()); }
107  inline Index col() const { return (MatrixType::Flags&RowMajorBit ? this->index() : Base::outer()); }
108  inline StorageIndex index() const
109  {
110  if(HasUnitDiag && m_returnOne) return Base::outer();
111  else return Base::index();
112  }
113  inline Scalar value() const
114  {
115  if(HasUnitDiag && m_returnOne) return Scalar(1);
116  else return Base::value();
117  }
118 
119  EIGEN_STRONG_INLINE operator bool() const
120  {
121  if(HasUnitDiag && m_returnOne)
122  return true;
123  if(SkipFirst) return Base::operator bool();
124  else
125  {
126  if (SkipDiag) return (Base::operator bool() && this->index() < this->outer());
127  else return (Base::operator bool() && this->index() <= this->outer());
128  }
129  }
130  protected:
131  bool m_returnOne;
132 };
133 
134 template<typename MatrixType, unsigned int Mode>
135 class TriangularViewImpl<MatrixType,Mode,Sparse>::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator
136 {
137  typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base;
138  public:
139 
140  EIGEN_STRONG_INLINE ReverseInnerIterator(const TriangularViewType& view, Index outer)
141  : Base(view.derived().nestedExpression(), outer)
142  {
143  eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal");
144  if(SkipLast) {
145  while((*this) && (SkipDiag ? this->index()>=outer : this->index()>outer))
146  --(*this);
147  }
148  }
149 
150  EIGEN_STRONG_INLINE ReverseInnerIterator& operator--()
151  { Base::operator--(); return *this; }
152 
153  inline Index row() const { return Base::row(); }
154  inline Index col() const { return Base::col(); }
155 
156  EIGEN_STRONG_INLINE operator bool() const
157  {
158  if (SkipLast) return Base::operator bool() ;
159  else
160  {
161  if(SkipDiag) return (Base::operator bool() && this->index() > this->outer());
162  else return (Base::operator bool() && this->index() >= this->outer());
163  }
164  }
165 };
166 
167 namespace internal {
168 
169 template<typename ArgType, unsigned int Mode>
170 struct unary_evaluator<TriangularView<ArgType,Mode>, IteratorBased>
171  : evaluator_base<TriangularView<ArgType,Mode> >
172 {
173  typedef TriangularView<ArgType,Mode> XprType;
174 
175 protected:
176 
177  typedef typename XprType::Scalar Scalar;
178  typedef typename XprType::StorageIndex StorageIndex;
179  typedef typename evaluator<ArgType>::InnerIterator EvalIterator;
180 
181  enum { SkipFirst = ((Mode&Lower) && !(ArgType::Flags&RowMajorBit))
182  || ((Mode&Upper) && (ArgType::Flags&RowMajorBit)),
183  SkipLast = !SkipFirst,
184  SkipDiag = (Mode&ZeroDiag) ? 1 : 0,
185  HasUnitDiag = (Mode&UnitDiag) ? 1 : 0
186  };
187 
188 public:
189 
190  enum {
191  CoeffReadCost = evaluator<ArgType>::CoeffReadCost,
192  Flags = XprType::Flags
193  };
194 
195  explicit unary_evaluator(const XprType &xpr) : m_argImpl(xpr.nestedExpression()) {}
196 
197  inline Index nonZerosEstimate() const {
198  return m_argImpl.nonZerosEstimate();
199  }
200 
201  class InnerIterator : public EvalIterator
202  {
203  typedef EvalIterator Base;
204  public:
205 
206  EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& xprEval, Index outer)
207  : Base(xprEval.m_argImpl,outer), m_returnOne(false)
208  {
209  if(SkipFirst)
210  {
211  while((*this) && ((HasUnitDiag||SkipDiag) ? this->index()<=outer : this->index()<outer))
212  Base::operator++();
213  if(HasUnitDiag)
214  m_returnOne = true;
215  }
216  else if(HasUnitDiag && ((!Base::operator bool()) || Base::index()>=Base::outer()))
217  {
218  if((!SkipFirst) && Base::operator bool())
219  Base::operator++();
220  m_returnOne = true; // FIXME check innerSize()>outer();
221  }
222  }
223 
224  EIGEN_STRONG_INLINE InnerIterator& operator++()
225  {
226  if(HasUnitDiag && m_returnOne)
227  m_returnOne = false;
228  else
229  {
230  Base::operator++();
231  if(HasUnitDiag && (!SkipFirst) && ((!Base::operator bool()) || Base::index()>=Base::outer()))
232  {
233  if((!SkipFirst) && Base::operator bool())
234  Base::operator++();
235  m_returnOne = true; // FIXME check innerSize()>outer();
236  }
237  }
238  return *this;
239  }
240 
241  EIGEN_STRONG_INLINE operator bool() const
242  {
243  if(HasUnitDiag && m_returnOne)
244  return true;
245  if(SkipFirst) return Base::operator bool();
246  else
247  {
248  if (SkipDiag) return (Base::operator bool() && this->index() < this->outer());
249  else return (Base::operator bool() && this->index() <= this->outer());
250  }
251  }
252 
253 // inline Index row() const { return (ArgType::Flags&RowMajorBit ? Base::outer() : this->index()); }
254 // inline Index col() const { return (ArgType::Flags&RowMajorBit ? this->index() : Base::outer()); }
255  inline StorageIndex index() const
256  {
257  if(HasUnitDiag && m_returnOne) return internal::convert_index<StorageIndex>(Base::outer());
258  else return Base::index();
259  }
260  inline Scalar value() const
261  {
262  if(HasUnitDiag && m_returnOne) return Scalar(1);
263  else return Base::value();
264  }
265 
266  protected:
267  bool m_returnOne;
268  private:
269  Scalar& valueRef();
270  };
271 
272 protected:
273  evaluator<ArgType> m_argImpl;
274 };
275 
276 } // end namespace internal
277 
278 template<typename Derived>
279 template<int Mode>
280 inline const TriangularView<const Derived, Mode>
281 SparseMatrixBase<Derived>::triangularView() const
282 {
283  return TriangularView<const Derived, Mode>(derived());
284 }
285 
286 } // end namespace Eigen
287 
288 #endif // EIGEN_SPARSE_TRIANGULARVIEW_H
Definition: Constants.h:196
Definition: Constants.h:202
Definition: LDLT.h:16
const unsigned int RowMajorBit
Definition: Constants.h:53
Definition: Constants.h:198
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
Definition: Constants.h:485
Definition: Constants.h:200
Definition: Eigen_Colamd.h:54
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48