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
14namespace Eigen {
15
16namespace internal {
17
18template<typename MatrixType, int Mode>
19struct traits<SparseTriangularView<MatrixType,Mode> >
20: public traits<MatrixType>
21{};
22
23} // namespace internal
24
25template<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
65template<typename MatrixType, int Mode>
66class 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
135template<typename MatrixType, int Mode>
136class 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
169template<typename Derived>
170template<int Mode>
171inline const SparseTriangularView<Derived, Mode>
172SparseMatrixBase<Derived>::triangularView() const
173{
174  return derived();
175}
176
177} // end namespace Eigen
178
179#endif // EIGEN_SPARSE_TRIANGULARVIEW_H
180