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