1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2007-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@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_DIAGONAL_H
12 #define EIGEN_DIAGONAL_H
13 
14 namespace Eigen {
15 
16 /** \class Diagonal
17   * \ingroup Core_Module
18   *
19   * \brief Expression of a diagonal/subdiagonal/superdiagonal in a matrix
20   *
21   * \param MatrixType the type of the object in which we are taking a sub/main/super diagonal
22   * \param DiagIndex the index of the sub/super diagonal. The default is 0 and it means the main diagonal.
23   *              A positive value means a superdiagonal, a negative value means a subdiagonal.
24   *              You can also use DynamicIndex so the index can be set at runtime.
25   *
26   * The matrix is not required to be square.
27   *
28   * This class represents an expression of the main diagonal, or any sub/super diagonal
29   * of a square matrix. It is the return type of MatrixBase::diagonal() and MatrixBase::diagonal(Index) and most of the
30   * time this is the only way it is used.
31   *
32   * \sa MatrixBase::diagonal(), MatrixBase::diagonal(Index)
33   */
34 
35 namespace internal {
36 template<typename MatrixType, int DiagIndex>
37 struct traits<Diagonal<MatrixType,DiagIndex> >
38  : traits<MatrixType>
39 {
40   typedef typename ref_selector<MatrixType>::type MatrixTypeNested;
41   typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested;
42   typedef typename MatrixType::StorageKind StorageKind;
43   enum {
44     RowsAtCompileTime = (int(DiagIndex) == DynamicIndex || int(MatrixType::SizeAtCompileTime) == Dynamic) ? Dynamic
45                       : (EIGEN_PLAIN_ENUM_MIN(MatrixType::RowsAtCompileTime - EIGEN_PLAIN_ENUM_MAX(-DiagIndex, 0),
46                                               MatrixType::ColsAtCompileTime - EIGEN_PLAIN_ENUM_MAX( DiagIndex, 0))),
47     ColsAtCompileTime = 1,
48     MaxRowsAtCompileTime = int(MatrixType::MaxSizeAtCompileTime) == Dynamic ? Dynamic
49                          : DiagIndex == DynamicIndex ? EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::MaxRowsAtCompileTime,
50                                                                               MatrixType::MaxColsAtCompileTime)
51                          : (EIGEN_PLAIN_ENUM_MIN(MatrixType::MaxRowsAtCompileTime - EIGEN_PLAIN_ENUM_MAX(-DiagIndex, 0),
52                                                  MatrixType::MaxColsAtCompileTime - EIGEN_PLAIN_ENUM_MAX( DiagIndex, 0))),
53     MaxColsAtCompileTime = 1,
54     MaskLvalueBit = is_lvalue<MatrixType>::value ? LvalueBit : 0,
55     Flags = (unsigned int)_MatrixTypeNested::Flags & (RowMajorBit | MaskLvalueBit | DirectAccessBit) & ~RowMajorBit, // FIXME DirectAccessBit should not be handled by expressions
56     MatrixTypeOuterStride = outer_stride_at_compile_time<MatrixType>::ret,
57     InnerStrideAtCompileTime = MatrixTypeOuterStride == Dynamic ? Dynamic : MatrixTypeOuterStride+1,
58     OuterStrideAtCompileTime = 0
59   };
60 };
61 }
62 
63 template<typename MatrixType, int _DiagIndex> class Diagonal
64    : public internal::dense_xpr_base< Diagonal<MatrixType,_DiagIndex> >::type
65 {
66   public:
67 
68     enum { DiagIndex = _DiagIndex };
69     typedef typename internal::dense_xpr_base<Diagonal>::type Base;
70     EIGEN_DENSE_PUBLIC_INTERFACE(Diagonal)
71 
72     EIGEN_DEVICE_FUNC
73     explicit inline Diagonal(MatrixType& matrix, Index a_index = DiagIndex) : m_matrix(matrix), m_index(a_index) {}
74 
75     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Diagonal)
76 
77     EIGEN_DEVICE_FUNC
78     inline Index rows() const
79     {
80       return m_index.value()<0 ? numext::mini<Index>(m_matrix.cols(),m_matrix.rows()+m_index.value())
81                                : numext::mini<Index>(m_matrix.rows(),m_matrix.cols()-m_index.value());
82     }
83 
84     EIGEN_DEVICE_FUNC
85     inline Index cols() const { return 1; }
86 
87     EIGEN_DEVICE_FUNC
88     inline Index innerStride() const
89     {
90       return m_matrix.outerStride() + 1;
91     }
92 
93     EIGEN_DEVICE_FUNC
94     inline Index outerStride() const
95     {
96       return 0;
97     }
98 
99     typedef typename internal::conditional<
100                        internal::is_lvalue<MatrixType>::value,
101                        Scalar,
102                        const Scalar
103                      >::type ScalarWithConstIfNotLvalue;
104 
105     EIGEN_DEVICE_FUNC
106     inline ScalarWithConstIfNotLvalue* data() { return &(m_matrix.coeffRef(rowOffset(), colOffset())); }
107     EIGEN_DEVICE_FUNC
108     inline const Scalar* data() const { return &(m_matrix.coeffRef(rowOffset(), colOffset())); }
109 
110     EIGEN_DEVICE_FUNC
111     inline Scalar& coeffRef(Index row, Index)
112     {
113       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
114       return m_matrix.coeffRef(row+rowOffset(), row+colOffset());
115     }
116 
117     EIGEN_DEVICE_FUNC
118     inline const Scalar& coeffRef(Index row, Index) const
119     {
120       return m_matrix.coeffRef(row+rowOffset(), row+colOffset());
121     }
122 
123     EIGEN_DEVICE_FUNC
124     inline CoeffReturnType coeff(Index row, Index) const
125     {
126       return m_matrix.coeff(row+rowOffset(), row+colOffset());
127     }
128 
129     EIGEN_DEVICE_FUNC
130     inline Scalar& coeffRef(Index idx)
131     {
132       EIGEN_STATIC_ASSERT_LVALUE(MatrixType)
133       return m_matrix.coeffRef(idx+rowOffset(), idx+colOffset());
134     }
135 
136     EIGEN_DEVICE_FUNC
137     inline const Scalar& coeffRef(Index idx) const
138     {
139       return m_matrix.coeffRef(idx+rowOffset(), idx+colOffset());
140     }
141 
142     EIGEN_DEVICE_FUNC
143     inline CoeffReturnType coeff(Index idx) const
144     {
145       return m_matrix.coeff(idx+rowOffset(), idx+colOffset());
146     }
147 
148     EIGEN_DEVICE_FUNC
149     inline const typename internal::remove_all<typename MatrixType::Nested>::type&
150     nestedExpression() const
151     {
152       return m_matrix;
153     }
154 
155     EIGEN_DEVICE_FUNC
156     inline Index index() const
157     {
158       return m_index.value();
159     }
160 
161   protected:
162     typename internal::ref_selector<MatrixType>::non_const_type m_matrix;
163     const internal::variable_if_dynamicindex<Index, DiagIndex> m_index;
164 
165   private:
166     // some compilers may fail to optimize std::max etc in case of compile-time constants...
167     EIGEN_DEVICE_FUNC
168     EIGEN_STRONG_INLINE Index absDiagIndex() const { return m_index.value()>0 ? m_index.value() : -m_index.value(); }
169     EIGEN_DEVICE_FUNC
170     EIGEN_STRONG_INLINE Index rowOffset() const { return m_index.value()>0 ? 0 : -m_index.value(); }
171     EIGEN_DEVICE_FUNC
172     EIGEN_STRONG_INLINE Index colOffset() const { return m_index.value()>0 ? m_index.value() : 0; }
173     // trigger a compile-time error if someone try to call packet
174     template<int LoadMode> typename MatrixType::PacketReturnType packet(Index) const;
175     template<int LoadMode> typename MatrixType::PacketReturnType packet(Index,Index) const;
176 };
177 
178 /** \returns an expression of the main diagonal of the matrix \c *this
179   *
180   * \c *this is not required to be square.
181   *
182   * Example: \include MatrixBase_diagonal.cpp
183   * Output: \verbinclude MatrixBase_diagonal.out
184   *
185   * \sa class Diagonal */
186 template<typename Derived>
187 inline typename MatrixBase<Derived>::DiagonalReturnType
188 MatrixBase<Derived>::diagonal()
189 {
190   return DiagonalReturnType(derived());
191 }
192 
193 /** This is the const version of diagonal(). */
194 template<typename Derived>
195 inline typename MatrixBase<Derived>::ConstDiagonalReturnType
196 MatrixBase<Derived>::diagonal() const
197 {
198   return ConstDiagonalReturnType(derived());
199 }
200 
201 /** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this
202   *
203   * \c *this is not required to be square.
204   *
205   * The template parameter \a DiagIndex represent a super diagonal if \a DiagIndex > 0
206   * and a sub diagonal otherwise. \a DiagIndex == 0 is equivalent to the main diagonal.
207   *
208   * Example: \include MatrixBase_diagonal_int.cpp
209   * Output: \verbinclude MatrixBase_diagonal_int.out
210   *
211   * \sa MatrixBase::diagonal(), class Diagonal */
212 template<typename Derived>
213 inline typename MatrixBase<Derived>::DiagonalDynamicIndexReturnType
214 MatrixBase<Derived>::diagonal(Index index)
215 {
216   return DiagonalDynamicIndexReturnType(derived(), index);
217 }
218 
219 /** This is the const version of diagonal(Index). */
220 template<typename Derived>
221 inline typename MatrixBase<Derived>::ConstDiagonalDynamicIndexReturnType
222 MatrixBase<Derived>::diagonal(Index index) const
223 {
224   return ConstDiagonalDynamicIndexReturnType(derived(), index);
225 }
226 
227 /** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this
228   *
229   * \c *this is not required to be square.
230   *
231   * The template parameter \a DiagIndex represent a super diagonal if \a DiagIndex > 0
232   * and a sub diagonal otherwise. \a DiagIndex == 0 is equivalent to the main diagonal.
233   *
234   * Example: \include MatrixBase_diagonal_template_int.cpp
235   * Output: \verbinclude MatrixBase_diagonal_template_int.out
236   *
237   * \sa MatrixBase::diagonal(), class Diagonal */
238 template<typename Derived>
239 template<int Index_>
240 inline typename MatrixBase<Derived>::template DiagonalIndexReturnType<Index_>::Type
241 MatrixBase<Derived>::diagonal()
242 {
243   return typename DiagonalIndexReturnType<Index_>::Type(derived());
244 }
245 
246 /** This is the const version of diagonal<int>(). */
247 template<typename Derived>
248 template<int Index_>
249 inline typename MatrixBase<Derived>::template ConstDiagonalIndexReturnType<Index_>::Type
250 MatrixBase<Derived>::diagonal() const
251 {
252   return typename ConstDiagonalIndexReturnType<Index_>::Type(derived());
253 }
254 
255 } // end namespace Eigen
256 
257 #endif // EIGEN_DIAGONAL_H
258