1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 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 
16 namespace internal {
17 
18 template<typename MatrixType, int Mode>
19 struct traits<SparseTriangularView<MatrixType,Mode> >
20 : public traits<MatrixType>
21 {};
22 
23 } // namespace internal
24 
25 template<typename MatrixType, int Mode> class SparseTriangularView
26   : public SparseMatrixBase<SparseTriangularView<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 
35   public:
36 
37     EIGEN_SPARSE_PUBLIC_INTERFACE(SparseTriangularView)
38 
39     class InnerIterator;
40     class ReverseInnerIterator;
41 
42     inline Index rows() const { return m_matrix.rows(); }
43     inline Index cols() const { return m_matrix.cols(); }
44 
45     typedef typename MatrixType::Nested MatrixTypeNested;
46     typedef typename internal::remove_reference<MatrixTypeNested>::type MatrixTypeNestedNonRef;
47     typedef typename internal::remove_all<MatrixTypeNested>::type MatrixTypeNestedCleaned;
48 
49     inline SparseTriangularView(const MatrixType& matrix) : m_matrix(matrix) {}
50 
51     /** \internal */
52     inline const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; }
53 
54     template<typename OtherDerived>
55     typename internal::plain_matrix_type_column_major<OtherDerived>::type
56     solve(const MatrixBase<OtherDerived>& other) const;
57 
58     template<typename OtherDerived> void solveInPlace(MatrixBase<OtherDerived>& other) const;
59     template<typename OtherDerived> void solveInPlace(SparseMatrixBase<OtherDerived>& other) const;
60 
61   protected:
62     MatrixTypeNested m_matrix;
63 };
64 
65 template<typename MatrixType, int Mode>
66 class SparseTriangularView<MatrixType,Mode>::InnerIterator : public MatrixTypeNestedCleaned::InnerIterator
67 {
68     typedef typename MatrixTypeNestedCleaned::InnerIterator Base;
69     typedef typename SparseTriangularView::Index Index;
70   public:
71 
72     EIGEN_STRONG_INLINE InnerIterator(const SparseTriangularView& view, Index outer)
73       : Base(view.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 Index 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, int Mode>
136 class SparseTriangularView<MatrixType,Mode>::ReverseInnerIterator : public MatrixTypeNestedCleaned::ReverseInnerIterator
137 {
138     typedef typename MatrixTypeNestedCleaned::ReverseInnerIterator Base;
139     typedef typename SparseTriangularView::Index Index;
140   public:
141 
142     EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTriangularView& view, Index outer)
143       : Base(view.nestedExpression(), outer)
144     {
145       eigen_assert((!HasUnitDiag) && "ReverseInnerIterator does not support yet triangular views with a unit diagonal");
146       if(SkipLast) {
147         while((*this) && (SkipDiag ? this->index()>=outer : this->index()>outer))
148           --(*this);
149       }
150     }
151 
152     EIGEN_STRONG_INLINE ReverseInnerIterator& operator--()
153     { Base::operator--(); return *this; }
154 
155     inline Index row() const { return Base::row(); }
156     inline Index col() const { return Base::col(); }
157 
158     EIGEN_STRONG_INLINE operator bool() const
159     {
160       if (SkipLast) return Base::operator bool() ;
161       else
162       {
163         if(SkipDiag) return (Base::operator bool() && this->index() > this->outer());
164         else return (Base::operator bool() && this->index() >= this->outer());
165       }
166     }
167 };
168 
169 template<typename Derived>
170 template<int Mode>
171 inline const SparseTriangularView<Derived, Mode>
172 SparseMatrixBase<Derived>::triangularView() const
173 {
174   return derived();
175 }
176 
177 } // end namespace Eigen
178 
179 #endif // EIGEN_SPARSE_TRIANGULARVIEW_H
180