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 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_BANDMATRIX_H
11 #define EIGEN_BANDMATRIX_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template<typename Derived>
18 class BandMatrixBase : public EigenBase<Derived>
19 {
20   public:
21 
22     enum {
23       Flags = internal::traits<Derived>::Flags,
24       CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
25       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
26       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
27       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
28       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
29       Supers = internal::traits<Derived>::Supers,
30       Subs   = internal::traits<Derived>::Subs,
31       Options = internal::traits<Derived>::Options
32     };
33     typedef typename internal::traits<Derived>::Scalar Scalar;
34     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> DenseMatrixType;
35     typedef typename DenseMatrixType::Index Index;
36     typedef typename internal::traits<Derived>::CoefficientsType CoefficientsType;
37     typedef EigenBase<Derived> Base;
38 
39   protected:
40     enum {
41       DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic))
42                             ? 1 + Supers + Subs
43                             : Dynamic,
44       SizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime)
45     };
46 
47   public:
48 
49     using Base::derived;
50     using Base::rows;
51     using Base::cols;
52 
53     /** \returns the number of super diagonals */
supers()54     inline Index supers() const { return derived().supers(); }
55 
56     /** \returns the number of sub diagonals */
subs()57     inline Index subs() const { return derived().subs(); }
58 
59     /** \returns an expression of the underlying coefficient matrix */
coeffs()60     inline const CoefficientsType& coeffs() const { return derived().coeffs(); }
61 
62     /** \returns an expression of the underlying coefficient matrix */
coeffs()63     inline CoefficientsType& coeffs() { return derived().coeffs(); }
64 
65     /** \returns a vector expression of the \a i -th column,
66       * only the meaningful part is returned.
67       * \warning the internal storage must be column major. */
col(Index i)68     inline Block<CoefficientsType,Dynamic,1> col(Index i)
69     {
70       EIGEN_STATIC_ASSERT((Options&RowMajor)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
71       Index start = 0;
72       Index len = coeffs().rows();
73       if (i<=supers())
74       {
75         start = supers()-i;
76         len = (std::min)(rows(),std::max<Index>(0,coeffs().rows() - (supers()-i)));
77       }
78       else if (i>=rows()-subs())
79         len = std::max<Index>(0,coeffs().rows() - (i + 1 - rows() + subs()));
80       return Block<CoefficientsType,Dynamic,1>(coeffs(), start, i, len, 1);
81     }
82 
83     /** \returns a vector expression of the main diagonal */
diagonal()84     inline Block<CoefficientsType,1,SizeAtCompileTime> diagonal()
85     { return Block<CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
86 
87     /** \returns a vector expression of the main diagonal (const version) */
diagonal()88     inline const Block<const CoefficientsType,1,SizeAtCompileTime> diagonal() const
89     { return Block<const CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
90 
91     template<int Index> struct DiagonalIntReturnType {
92       enum {
93         ReturnOpposite = (Options&SelfAdjoint) && (((Index)>0 && Supers==0) || ((Index)<0 && Subs==0)),
94         Conjugate = ReturnOpposite && NumTraits<Scalar>::IsComplex,
95         ActualIndex = ReturnOpposite ? -Index : Index,
96         DiagonalSize = (RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic)
97                      ? Dynamic
98                      : (ActualIndex<0
99                      ? EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime, RowsAtCompileTime + ActualIndex)
100                      : EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime - ActualIndex))
101       };
102       typedef Block<CoefficientsType,1, DiagonalSize> BuildType;
103       typedef typename internal::conditional<Conjugate,
104                  CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>,BuildType >,
105                  BuildType>::type Type;
106     };
107 
108     /** \returns a vector expression of the \a N -th sub or super diagonal */
diagonal()109     template<int N> inline typename DiagonalIntReturnType<N>::Type diagonal()
110     {
111       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
112     }
113 
114     /** \returns a vector expression of the \a N -th sub or super diagonal */
diagonal()115     template<int N> inline const typename DiagonalIntReturnType<N>::Type diagonal() const
116     {
117       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
118     }
119 
120     /** \returns a vector expression of the \a i -th sub or super diagonal */
diagonal(Index i)121     inline Block<CoefficientsType,1,Dynamic> diagonal(Index i)
122     {
123       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
124       return Block<CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
125     }
126 
127     /** \returns a vector expression of the \a i -th sub or super diagonal */
diagonal(Index i)128     inline const Block<const CoefficientsType,1,Dynamic> diagonal(Index i) const
129     {
130       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
131       return Block<const CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
132     }
133 
evalTo(Dest & dst)134     template<typename Dest> inline void evalTo(Dest& dst) const
135     {
136       dst.resize(rows(),cols());
137       dst.setZero();
138       dst.diagonal() = diagonal();
139       for (Index i=1; i<=supers();++i)
140         dst.diagonal(i) = diagonal(i);
141       for (Index i=1; i<=subs();++i)
142         dst.diagonal(-i) = diagonal(-i);
143     }
144 
toDenseMatrix()145     DenseMatrixType toDenseMatrix() const
146     {
147       DenseMatrixType res(rows(),cols());
148       evalTo(res);
149       return res;
150     }
151 
152   protected:
153 
diagonalLength(Index i)154     inline Index diagonalLength(Index i) const
155     { return i<0 ? (std::min)(cols(),rows()+i) : (std::min)(rows(),cols()-i); }
156 };
157 
158 /**
159   * \class BandMatrix
160   * \ingroup Core_Module
161   *
162   * \brief Represents a rectangular matrix with a banded storage
163   *
164   * \param _Scalar Numeric type, i.e. float, double, int
165   * \param Rows Number of rows, or \b Dynamic
166   * \param Cols Number of columns, or \b Dynamic
167   * \param Supers Number of super diagonal
168   * \param Subs Number of sub diagonal
169   * \param _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
170   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to
171   *                 column-major. The latter controls whether the matrix represents a selfadjoint
172   *                 matrix in which case either Supers of Subs have to be null.
173   *
174   * \sa class TridiagonalMatrix
175   */
176 
177 template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
178 struct traits<BandMatrix<_Scalar,_Rows,_Cols,_Supers,_Subs,_Options> >
179 {
180   typedef _Scalar Scalar;
181   typedef Dense StorageKind;
182   typedef DenseIndex Index;
183   enum {
184     CoeffReadCost = NumTraits<Scalar>::ReadCost,
185     RowsAtCompileTime = _Rows,
186     ColsAtCompileTime = _Cols,
187     MaxRowsAtCompileTime = _Rows,
188     MaxColsAtCompileTime = _Cols,
189     Flags = LvalueBit,
190     Supers = _Supers,
191     Subs = _Subs,
192     Options = _Options,
193     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
194   };
195   typedef Matrix<Scalar,DataRowsAtCompileTime,ColsAtCompileTime,Options&RowMajor?RowMajor:ColMajor> CoefficientsType;
196 };
197 
198 template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
199 class BandMatrix : public BandMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
200 {
201   public:
202 
203     typedef typename internal::traits<BandMatrix>::Scalar Scalar;
204     typedef typename internal::traits<BandMatrix>::Index Index;
205     typedef typename internal::traits<BandMatrix>::CoefficientsType CoefficientsType;
206 
207     inline BandMatrix(Index rows=Rows, Index cols=Cols, Index supers=Supers, Index subs=Subs)
208       : m_coeffs(1+supers+subs,cols),
209         m_rows(rows), m_supers(supers), m_subs(subs)
210     {
211     }
212 
213     /** \returns the number of columns */
214     inline Index rows() const { return m_rows.value(); }
215 
216     /** \returns the number of rows */
217     inline Index cols() const { return m_coeffs.cols(); }
218 
219     /** \returns the number of super diagonals */
220     inline Index supers() const { return m_supers.value(); }
221 
222     /** \returns the number of sub diagonals */
223     inline Index subs() const { return m_subs.value(); }
224 
225     inline const CoefficientsType& coeffs() const { return m_coeffs; }
226     inline CoefficientsType& coeffs() { return m_coeffs; }
227 
228   protected:
229 
230     CoefficientsType m_coeffs;
231     internal::variable_if_dynamic<Index, Rows>   m_rows;
232     internal::variable_if_dynamic<Index, Supers> m_supers;
233     internal::variable_if_dynamic<Index, Subs>   m_subs;
234 };
235 
236 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
237 class BandMatrixWrapper;
238 
239 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
240 struct traits<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
241 {
242   typedef typename _CoefficientsType::Scalar Scalar;
243   typedef typename _CoefficientsType::StorageKind StorageKind;
244   typedef typename _CoefficientsType::Index Index;
245   enum {
246     CoeffReadCost = internal::traits<_CoefficientsType>::CoeffReadCost,
247     RowsAtCompileTime = _Rows,
248     ColsAtCompileTime = _Cols,
249     MaxRowsAtCompileTime = _Rows,
250     MaxColsAtCompileTime = _Cols,
251     Flags = LvalueBit,
252     Supers = _Supers,
253     Subs = _Subs,
254     Options = _Options,
255     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
256   };
257   typedef _CoefficientsType CoefficientsType;
258 };
259 
260 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
261 class BandMatrixWrapper : public BandMatrixBase<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
262 {
263   public:
264 
265     typedef typename internal::traits<BandMatrixWrapper>::Scalar Scalar;
266     typedef typename internal::traits<BandMatrixWrapper>::CoefficientsType CoefficientsType;
267     typedef typename internal::traits<BandMatrixWrapper>::Index Index;
268 
269     inline BandMatrixWrapper(const CoefficientsType& coeffs, Index rows=_Rows, Index cols=_Cols, Index supers=_Supers, Index subs=_Subs)
270       : m_coeffs(coeffs),
271         m_rows(rows), m_supers(supers), m_subs(subs)
272     {
273       EIGEN_UNUSED_VARIABLE(cols);
274       //internal::assert(coeffs.cols()==cols() && (supers()+subs()+1)==coeffs.rows());
275     }
276 
277     /** \returns the number of columns */
278     inline Index rows() const { return m_rows.value(); }
279 
280     /** \returns the number of rows */
281     inline Index cols() const { return m_coeffs.cols(); }
282 
283     /** \returns the number of super diagonals */
284     inline Index supers() const { return m_supers.value(); }
285 
286     /** \returns the number of sub diagonals */
287     inline Index subs() const { return m_subs.value(); }
288 
289     inline const CoefficientsType& coeffs() const { return m_coeffs; }
290 
291   protected:
292 
293     const CoefficientsType& m_coeffs;
294     internal::variable_if_dynamic<Index, _Rows>   m_rows;
295     internal::variable_if_dynamic<Index, _Supers> m_supers;
296     internal::variable_if_dynamic<Index, _Subs>   m_subs;
297 };
298 
299 /**
300   * \class TridiagonalMatrix
301   * \ingroup Core_Module
302   *
303   * \brief Represents a tridiagonal matrix with a compact banded storage
304   *
305   * \param _Scalar Numeric type, i.e. float, double, int
306   * \param Size Number of rows and cols, or \b Dynamic
307   * \param _Options Can be 0 or \b SelfAdjoint
308   *
309   * \sa class BandMatrix
310   */
311 template<typename Scalar, int Size, int Options>
312 class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>
313 {
314     typedef BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor> Base;
315     typedef typename Base::Index Index;
316   public:
317     TridiagonalMatrix(Index size = Size) : Base(size,size,Options&SelfAdjoint?0:1,1) {}
318 
319     inline typename Base::template DiagonalIntReturnType<1>::Type super()
320     { return Base::template diagonal<1>(); }
321     inline const typename Base::template DiagonalIntReturnType<1>::Type super() const
322     { return Base::template diagonal<1>(); }
323     inline typename Base::template DiagonalIntReturnType<-1>::Type sub()
324     { return Base::template diagonal<-1>(); }
325     inline const typename Base::template DiagonalIntReturnType<-1>::Type sub() const
326     { return Base::template diagonal<-1>(); }
327   protected:
328 };
329 
330 } // end namespace internal
331 
332 } // end namespace Eigen
333 
334 #endif // EIGEN_BANDMATRIX_H
335