src/Core/BandMatrix.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_BANDMATRIX_H
0011 #define EIGEN_BANDMATRIX_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 template<typename Derived>
0018 class BandMatrixBase : public EigenBase<Derived>
0019 {
0020   public:
0021
0022     enum {
0023       Flags = internal::traits<Derived>::Flags,
0024       CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
0025       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
0026       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
0027       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
0028       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
0029       Supers = internal::traits<Derived>::Supers,
0030       Subs   = internal::traits<Derived>::Subs,
0031       Options = internal::traits<Derived>::Options
0032     };
0033     typedef typename internal::traits<Derived>::Scalar Scalar;
0034     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> DenseMatrixType;
0035     typedef typename DenseMatrixType::StorageIndex StorageIndex;
0036     typedef typename internal::traits<Derived>::CoefficientsType CoefficientsType;
0037     typedef EigenBase<Derived> Base;
0038
0039   protected:
0040     enum {
0041       DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic))
0042                             ? 1 + Supers + Subs
0043                             : Dynamic,
0044       SizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime)
0045     };
0046
0047   public:
0048
0049     using Base::derived;
0050     using Base::rows;
0051     using Base::cols;
0052
0053     /** \returns the number of super diagonals */
0054     inline Index supers() const { return derived().supers(); }
0055
0056     /** \returns the number of sub diagonals */
0057     inline Index subs() const { return derived().subs(); }
0058
0059     /** \returns an expression of the underlying coefficient matrix */
0060     inline const CoefficientsType& coeffs() const { return derived().coeffs(); }
0061
0062     /** \returns an expression of the underlying coefficient matrix */
0063     inline CoefficientsType& coeffs() { return derived().coeffs(); }
0064
0065     /** \returns a vector expression of the \a i -th column,
0066       * only the meaningful part is returned.
0067       * \warning the internal storage must be column major. */
0068     inline Block<CoefficientsType,Dynamic,1> col(Index i)
0069     {
0070       EIGEN_STATIC_ASSERT((int(Options) & int(RowMajor)) == 0, THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
0071       Index start = 0;
0072       Index len = coeffs().rows();
0073       if (i<=supers())
0074       {
0075         start = supers()-i;
0076         len = (std::min)(rows(),std::max<Index>(0,coeffs().rows() - (supers()-i)));
0077       }
0078       else if (i>=rows()-subs())
0079         len = std::max<Index>(0,coeffs().rows() - (i + 1 - rows() + subs()));
0080       return Block<CoefficientsType,Dynamic,1>(coeffs(), start, i, len, 1);
0081     }
0082
0083     /** \returns a vector expression of the main diagonal */
0084     inline Block<CoefficientsType,1,SizeAtCompileTime> diagonal()
0085     { return Block<CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
0086
0087     /** \returns a vector expression of the main diagonal (const version) */
0088     inline const Block<const CoefficientsType,1,SizeAtCompileTime> diagonal() const
0089     { return Block<const CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
0090
0091     template<int Index> struct DiagonalIntReturnType {
0092       enum {
0093         ReturnOpposite = (int(Options) & int(SelfAdjoint)) && (((Index) > 0 && Supers == 0) || ((Index) < 0 && Subs == 0)),
0094         Conjugate = ReturnOpposite && NumTraits<Scalar>::IsComplex,
0095         ActualIndex = ReturnOpposite ? -Index : Index,
0096         DiagonalSize = (RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic)
0097                      ? Dynamic
0098                      : (ActualIndex<0
0099                      ? EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime, RowsAtCompileTime + ActualIndex)
0100                      : EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime - ActualIndex))
0101       };
0102       typedef Block<CoefficientsType,1, DiagonalSize> BuildType;
0103       typedef typename internal::conditional<Conjugate,
0104                  CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>,BuildType >,
0105                  BuildType>::type Type;
0106     };
0107
0108     /** \returns a vector expression of the \a N -th sub or super diagonal */
0109     template<int N> inline typename DiagonalIntReturnType<N>::Type diagonal()
0110     {
0111       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
0112     }
0113
0114     /** \returns a vector expression of the \a N -th sub or super diagonal */
0115     template<int N> inline const typename DiagonalIntReturnType<N>::Type diagonal() const
0116     {
0117       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
0118     }
0119
0120     /** \returns a vector expression of the \a i -th sub or super diagonal */
0121     inline Block<CoefficientsType,1,Dynamic> diagonal(Index i)
0122     {
0123       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
0124       return Block<CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
0125     }
0126
0127     /** \returns a vector expression of the \a i -th sub or super diagonal */
0128     inline const Block<const CoefficientsType,1,Dynamic> diagonal(Index i) const
0129     {
0130       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
0131       return Block<const CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
0132     }
0133
0134     template<typename Dest> inline void evalTo(Dest& dst) const
0135     {
0136       dst.resize(rows(),cols());
0137       dst.setZero();
0138       dst.diagonal() = diagonal();
0139       for (Index i=1; i<=supers();++i)
0140         dst.diagonal(i) = diagonal(i);
0141       for (Index i=1; i<=subs();++i)
0142         dst.diagonal(-i) = diagonal(-i);
0143     }
0144
0145     DenseMatrixType toDenseMatrix() const
0146     {
0147       DenseMatrixType res(rows(),cols());
0148       evalTo(res);
0149       return res;
0150     }
0151
0152   protected:
0153
0154     inline Index diagonalLength(Index i) const
0155     { return i<0 ? (std::min)(cols(),rows()+i) : (std::min)(rows(),cols()-i); }
0156 };
0157
0158 /**
0159   * \class BandMatrix
0160   * \ingroup Core_Module
0161   *
0162   * \brief Represents a rectangular matrix with a banded storage
0163   *
0164   * \tparam _Scalar Numeric type, i.e. float, double, int
0165   * \tparam _Rows Number of rows, or \b Dynamic
0166   * \tparam _Cols Number of columns, or \b Dynamic
0167   * \tparam _Supers Number of super diagonal
0168   * \tparam _Subs Number of sub diagonal
0169   * \tparam _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
0170   *                  The former controls \ref TopicStorageOrders "storage order", and defaults to
0171   *                  column-major. The latter controls whether the matrix represents a selfadjoint
0172   *                  matrix in which case either Supers of Subs have to be null.
0173   *
0174   * \sa class TridiagonalMatrix
0175   */
0176
0177 template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
0178 struct traits<BandMatrix<_Scalar,_Rows,_Cols,_Supers,_Subs,_Options> >
0179 {
0180   typedef _Scalar Scalar;
0181   typedef Dense StorageKind;
0182   typedef Eigen::Index StorageIndex;
0183   enum {
0184     CoeffReadCost = NumTraits<Scalar>::ReadCost,
0185     RowsAtCompileTime = _Rows,
0186     ColsAtCompileTime = _Cols,
0187     MaxRowsAtCompileTime = _Rows,
0188     MaxColsAtCompileTime = _Cols,
0189     Flags = LvalueBit,
0190     Supers = _Supers,
0191     Subs = _Subs,
0192     Options = _Options,
0193     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
0194   };
0195   typedef Matrix<Scalar, DataRowsAtCompileTime, ColsAtCompileTime, int(Options) & int(RowMajor) ? RowMajor : ColMajor> CoefficientsType;
0196 };
0197
0198 template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
0199 class BandMatrix : public BandMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
0200 {
0201   public:
0202
0203     typedef typename internal::traits<BandMatrix>::Scalar Scalar;
0204     typedef typename internal::traits<BandMatrix>::StorageIndex StorageIndex;
0205     typedef typename internal::traits<BandMatrix>::CoefficientsType CoefficientsType;
0206
0207     explicit inline BandMatrix(Index rows=Rows, Index cols=Cols, Index supers=Supers, Index subs=Subs)
0208       : m_coeffs(1+supers+subs,cols),
0209         m_rows(rows), m_supers(supers), m_subs(subs)
0210     {
0211     }
0212
0213     /** \returns the number of columns */
0214     inline EIGEN_CONSTEXPR Index rows() const { return m_rows.value(); }
0215
0216     /** \returns the number of rows */
0217     inline EIGEN_CONSTEXPR Index cols() const { return m_coeffs.cols(); }
0218
0219     /** \returns the number of super diagonals */
0220     inline EIGEN_CONSTEXPR Index supers() const { return m_supers.value(); }
0221
0222     /** \returns the number of sub diagonals */
0223     inline EIGEN_CONSTEXPR Index subs() const { return m_subs.value(); }
0224
0225     inline const CoefficientsType& coeffs() const { return m_coeffs; }
0226     inline CoefficientsType& coeffs() { return m_coeffs; }
0227
0228   protected:
0229
0230     CoefficientsType m_coeffs;
0231     internal::variable_if_dynamic<Index, Rows>   m_rows;
0232     internal::variable_if_dynamic<Index, Supers> m_supers;
0233     internal::variable_if_dynamic<Index, Subs>   m_subs;
0234 };
0235
0236 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
0237 class BandMatrixWrapper;
0238
0239 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
0240 struct traits<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
0241 {
0242   typedef typename _CoefficientsType::Scalar Scalar;
0243   typedef typename _CoefficientsType::StorageKind StorageKind;
0244   typedef typename _CoefficientsType::StorageIndex StorageIndex;
0245   enum {
0246     CoeffReadCost = internal::traits<_CoefficientsType>::CoeffReadCost,
0247     RowsAtCompileTime = _Rows,
0248     ColsAtCompileTime = _Cols,
0249     MaxRowsAtCompileTime = _Rows,
0250     MaxColsAtCompileTime = _Cols,
0251     Flags = LvalueBit,
0252     Supers = _Supers,
0253     Subs = _Subs,
0254     Options = _Options,
0255     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
0256   };
0257   typedef _CoefficientsType CoefficientsType;
0258 };
0259
0260 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
0261 class BandMatrixWrapper : public BandMatrixBase<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
0262 {
0263   public:
0264
0265     typedef typename internal::traits<BandMatrixWrapper>::Scalar Scalar;
0266     typedef typename internal::traits<BandMatrixWrapper>::CoefficientsType CoefficientsType;
0267     typedef typename internal::traits<BandMatrixWrapper>::StorageIndex StorageIndex;
0268
0269     explicit inline BandMatrixWrapper(const CoefficientsType& coeffs, Index rows=_Rows, Index cols=_Cols, Index supers=_Supers, Index subs=_Subs)
0270       : m_coeffs(coeffs),
0271         m_rows(rows), m_supers(supers), m_subs(subs)
0272     {
0273       EIGEN_UNUSED_VARIABLE(cols);
0274       //internal::assert(coeffs.cols()==cols() && (supers()+subs()+1)==coeffs.rows());
0275     }
0276
0277     /** \returns the number of columns */
0278     inline EIGEN_CONSTEXPR Index rows() const { return m_rows.value(); }
0279
0280     /** \returns the number of rows */
0281     inline EIGEN_CONSTEXPR Index cols() const { return m_coeffs.cols(); }
0282
0283     /** \returns the number of super diagonals */
0284     inline EIGEN_CONSTEXPR Index supers() const { return m_supers.value(); }
0285
0286     /** \returns the number of sub diagonals */
0287     inline EIGEN_CONSTEXPR Index subs() const { return m_subs.value(); }
0288
0289     inline const CoefficientsType& coeffs() const { return m_coeffs; }
0290
0291   protected:
0292
0293     const CoefficientsType& m_coeffs;
0294     internal::variable_if_dynamic<Index, _Rows>   m_rows;
0295     internal::variable_if_dynamic<Index, _Supers> m_supers;
0296     internal::variable_if_dynamic<Index, _Subs>   m_subs;
0297 };
0298
0299 /**
0300   * \class TridiagonalMatrix
0301   * \ingroup Core_Module
0302   *
0303   * \brief Represents a tridiagonal matrix with a compact banded storage
0304   *
0305   * \tparam Scalar Numeric type, i.e. float, double, int
0306   * \tparam Size Number of rows and cols, or \b Dynamic
0307   * \tparam Options Can be 0 or \b SelfAdjoint
0308   *
0309   * \sa class BandMatrix
0310   */
0311 template<typename Scalar, int Size, int Options>
0312 class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>
0313 {
0314     typedef BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor> Base;
0315     typedef typename Base::StorageIndex StorageIndex;
0316   public:
0317     explicit TridiagonalMatrix(Index size = Size) : Base(size,size,Options&SelfAdjoint?0:1,1) {}
0318
0319     inline typename Base::template DiagonalIntReturnType<1>::Type super()
0320     { return Base::template diagonal<1>(); }
0321     inline const typename Base::template DiagonalIntReturnType<1>::Type super() const
0322     { return Base::template diagonal<1>(); }
0323     inline typename Base::template DiagonalIntReturnType<-1>::Type sub()
0324     { return Base::template diagonal<-1>(); }
0325     inline const typename Base::template DiagonalIntReturnType<-1>::Type sub() const
0326     { return Base::template diagonal<-1>(); }
0327   protected:
0328 };
0329
0330
0331 struct BandShape {};
0332
0333 template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
0334 struct evaluator_traits<BandMatrix<_Scalar,_Rows,_Cols,_Supers,_Subs,_Options> >
0335   : public evaluator_traits_base<BandMatrix<_Scalar,_Rows,_Cols,_Supers,_Subs,_Options> >
0336 {
0337   typedef BandShape Shape;
0338 };
0339
0340 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
0341 struct evaluator_traits<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
0342   : public evaluator_traits_base<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
0343 {
0344   typedef BandShape Shape;
0345 };
0346
0347 template<> struct AssignmentKind<DenseShape,BandShape> { typedef EigenBase2EigenBase Kind; };
0348
0349 } // end namespace internal
0350
0351 } // end namespace Eigen
0352
0353 #endif // EIGEN_BANDMATRIX_H