BandMatrix.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_BANDMATRIX_H
00011 #define EIGEN_BANDMATRIX_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Derived>
00018 class BandMatrixBase : public EigenBase<Derived>
00019 {
00020   public:
00021 
00022     enum {
00023       Flags = internal::traits<Derived>::Flags,
00024       CoeffReadCost = internal::traits<Derived>::CoeffReadCost,
00025       RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime,
00026       ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime,
00027       MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime,
00028       MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime,
00029       Supers = internal::traits<Derived>::Supers,
00030       Subs   = internal::traits<Derived>::Subs,
00031       Options = internal::traits<Derived>::Options
00032     };
00033     typedef typename internal::traits<Derived>::Scalar Scalar;
00034     typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> DenseMatrixType;
00035     typedef typename DenseMatrixType::Index Index;
00036     typedef typename internal::traits<Derived>::CoefficientsType CoefficientsType;
00037     typedef EigenBase<Derived> Base;
00038 
00039   protected:
00040     enum {
00041       DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic))
00042                             ? 1 + Supers + Subs
00043                             : Dynamic,
00044       SizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime)
00045     };
00046 
00047   public:
00048     
00049     using Base::derived ;
00050     using Base::rows ;
00051     using Base::cols ;
00052 
00053     /** \returns the number of super diagonals */
00054     inline Index supers() const { return derived ().supers(); }
00055 
00056     /** \returns the number of sub diagonals */
00057     inline Index subs() const { return derived ().subs(); }
00058     
00059     /** \returns an expression of the underlying coefficient matrix */
00060     inline const CoefficientsType& coeffs() const { return derived ().coeffs(); }
00061     
00062     /** \returns an expression of the underlying coefficient matrix */
00063     inline CoefficientsType& coeffs() { return derived ().coeffs(); }
00064 
00065     /** \returns a vector expression of the \a i -th column,
00066       * only the meaningful part is returned.
00067       * \warning the internal storage must be column major. */
00068     inline Block<CoefficientsType,Dynamic,1> col(Index i)
00069     {
00070       EIGEN_STATIC_ASSERT((Options&RowMajor)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00071       Index start = 0;
00072       Index len = coeffs().rows();
00073       if (i<=supers())
00074       {
00075         start = supers()-i;
00076         len = (std::min)(rows (),std::max<Index>(0,coeffs().rows() - (supers()-i)));
00077       }
00078       else if (i>=rows ()-subs())
00079         len = std::max<Index>(0,coeffs().rows() - (i + 1 - rows () + subs()));
00080       return Block<CoefficientsType,Dynamic,1>(coeffs(), start, i, len, 1);
00081     }
00082 
00083     /** \returns a vector expression of the main diagonal */
00084     inline Block<CoefficientsType,1,SizeAtCompileTime> diagonal()
00085     { return Block<CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows (),cols ())); }
00086 
00087     /** \returns a vector expression of the main diagonal (const version) */
00088     inline const Block<const CoefficientsType,1,SizeAtCompileTime> diagonal() const
00089     { return Block<const CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows (),cols ())); }
00090 
00091     template<int Index> struct DiagonalIntReturnType {
00092       enum {
00093         ReturnOpposite = (Options&SelfAdjoint) && (((Index)>0 && Supers==0) || ((Index)<0 && Subs==0)),
00094         Conjugate = ReturnOpposite && NumTraits<Scalar>::IsComplex,
00095         ActualIndex = ReturnOpposite ? -Index : Index,
00096         DiagonalSize = (RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic)
00097                      ? Dynamic
00098                      : (ActualIndex<0
00099                      ? EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime, RowsAtCompileTime + ActualIndex)
00100                      : EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime, ColsAtCompileTime - ActualIndex))
00101       };
00102       typedef Block<CoefficientsType,1, DiagonalSize> BuildType;
00103       typedef typename internal::conditional<Conjugate,
00104                  CwiseUnaryOp<internal::scalar_conjugate_op<Scalar>,BuildType >,
00105                  BuildType>::type Type;
00106     };
00107 
00108     /** \returns a vector expression of the \a N -th sub or super diagonal */
00109     template<int N> inline typename DiagonalIntReturnType<N>::Type diagonal()
00110     {
00111       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
00112     }
00113 
00114     /** \returns a vector expression of the \a N -th sub or super diagonal */
00115     template<int N> inline const typename DiagonalIntReturnType<N>::Type diagonal() const
00116     {
00117       return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
00118     }
00119 
00120     /** \returns a vector expression of the \a i -th sub or super diagonal */
00121     inline Block<CoefficientsType,1,Dynamic> diagonal(Index i)
00122     {
00123       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
00124       return Block<CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
00125     }
00126 
00127     /** \returns a vector expression of the \a i -th sub or super diagonal */
00128     inline const Block<const CoefficientsType,1,Dynamic> diagonal(Index i) const
00129     {
00130       eigen_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers()));
00131       return Block<const CoefficientsType,1,Dynamic>(coeffs(), supers()-i, std::max<Index>(0,i), 1, diagonalLength(i));
00132     }
00133     
00134     template<typename Dest> inline void evalTo(Dest& dst) const
00135     {
00136       dst.resize(rows(),cols());
00137       dst.setZero();
00138       dst.diagonal() = diagonal();
00139       for (Index i=1; i<=supers();++i)
00140         dst.diagonal(i) = diagonal(i);
00141       for (Index i=1; i<=subs();++i)
00142         dst.diagonal(-i) = diagonal(-i);
00143     }
00144 
00145     DenseMatrixType toDenseMatrix() const
00146     {
00147       DenseMatrixType res(rows(),cols());
00148       evalTo(res);
00149       return res;
00150     }
00151 
00152   protected:
00153 
00154     inline Index diagonalLength(Index i) const
00155     { return i<0 ? (std::min)(cols(),rows()+i) : (std::min)(rows(),cols()-i); }
00156 };
00157 
00158 /**
00159   * \class BandMatrix
00160   * \ingroup Core_Module
00161   *
00162   * \brief Represents a rectangular matrix with a banded storage
00163   *
00164   * \param _Scalar Numeric type, i.e. float, double, int
00165   * \param Rows Number of rows, or \b Dynamic
00166   * \param Cols Number of columns, or \b Dynamic
00167   * \param Supers Number of super diagonal
00168   * \param Subs Number of sub diagonal
00169   * \param _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
00170   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to
00171   *                 column-major. The latter controls whether the matrix represents a selfadjoint 
00172   *                 matrix in which case either Supers of Subs have to be null.
00173   *
00174   * \sa class TridiagonalMatrix
00175   */
00176 
00177 template<typename _Scalar, int _Rows, int _Cols, int _Supers, int _Subs, int _Options>
00178 struct traits<BandMatrix<_Scalar,_Rows,_Cols,_Supers,_Subs,_Options> >
00179 {
00180   typedef _Scalar Scalar;
00181   typedef Dense StorageKind;
00182   typedef DenseIndex Index;
00183   enum {
00184     CoeffReadCost = NumTraits<Scalar>::ReadCost,
00185     RowsAtCompileTime = _Rows,
00186     ColsAtCompileTime = _Cols,
00187     MaxRowsAtCompileTime = _Rows,
00188     MaxColsAtCompileTime = _Cols,
00189     Flags = LvalueBit,
00190     Supers = _Supers,
00191     Subs = _Subs,
00192     Options = _Options,
00193     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
00194   };
00195   typedef Matrix<Scalar,DataRowsAtCompileTime,ColsAtCompileTime,Options&RowMajor?RowMajor:ColMajor> CoefficientsType;
00196 };
00197 
00198 template<typename _Scalar, int Rows, int Cols, int Supers, int Subs, int Options>
00199 class BandMatrix : public BandMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs,Options> >
00200 {
00201   public:
00202 
00203     typedef typename internal::traits<BandMatrix>::Scalar Scalar;
00204     typedef typename internal::traits<BandMatrix>::Index Index;
00205     typedef typename internal::traits<BandMatrix>::CoefficientsType CoefficientsType;
00206 
00207     inline BandMatrix(Index rows=Rows, Index cols=Cols, Index supers=Supers, Index subs=Subs)
00208       : m_coeffs(1+supers+subs,cols),
00209         m_rows(rows), m_supers(supers), m_subs(subs)
00210     {
00211     }
00212 
00213     /** \returns the number of columns */
00214     inline Index rows () const { return m_rows.value(); }
00215 
00216     /** \returns the number of rows */
00217     inline Index cols () const { return m_coeffs.cols(); }
00218 
00219     /** \returns the number of super diagonals */
00220     inline Index supers () const { return m_supers.value(); }
00221 
00222     /** \returns the number of sub diagonals */
00223     inline Index subs () const { return m_subs.value(); }
00224 
00225     inline const CoefficientsType& coeffs() const { return m_coeffs; }
00226     inline CoefficientsType& coeffs() { return m_coeffs; }
00227 
00228   protected:
00229 
00230     CoefficientsType m_coeffs;
00231     internal::variable_if_dynamic<Index, Rows>   m_rows;
00232     internal::variable_if_dynamic<Index, Supers> m_supers;
00233     internal::variable_if_dynamic<Index, Subs>   m_subs;
00234 };
00235 
00236 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
00237 class BandMatrixWrapper;
00238 
00239 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
00240 struct traits<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
00241 {
00242   typedef typename _CoefficientsType::Scalar Scalar;
00243   typedef typename _CoefficientsType::StorageKind StorageKind;
00244   typedef typename _CoefficientsType::Index Index;
00245   enum {
00246     CoeffReadCost = internal::traits<_CoefficientsType>::CoeffReadCost,
00247     RowsAtCompileTime = _Rows,
00248     ColsAtCompileTime = _Cols,
00249     MaxRowsAtCompileTime = _Rows,
00250     MaxColsAtCompileTime = _Cols,
00251     Flags = LvalueBit,
00252     Supers = _Supers,
00253     Subs = _Subs,
00254     Options = _Options,
00255     DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) ? 1 + Supers + Subs : Dynamic
00256   };
00257   typedef _CoefficientsType CoefficientsType;
00258 };
00259 
00260 template<typename _CoefficientsType,int _Rows, int _Cols, int _Supers, int _Subs,int _Options>
00261 class BandMatrixWrapper : public BandMatrixBase<BandMatrixWrapper<_CoefficientsType,_Rows,_Cols,_Supers,_Subs,_Options> >
00262 {
00263   public:
00264 
00265     typedef typename internal::traits<BandMatrixWrapper>::Scalar Scalar;
00266     typedef typename internal::traits<BandMatrixWrapper>::CoefficientsType CoefficientsType;
00267     typedef typename internal::traits<BandMatrixWrapper>::Index Index;
00268 
00269     inline BandMatrixWrapper(const CoefficientsType& coeffs, Index rows=_Rows, Index cols=_Cols, Index supers=_Supers, Index subs=_Subs)
00270       : m_coeffs(coeffs),
00271         m_rows(rows), m_supers(supers), m_subs(subs)
00272     {
00273       EIGEN_UNUSED_VARIABLE(cols);
00274       //internal::assert(coeffs.cols()==cols() && (supers()+subs()+1)==coeffs.rows());
00275     }
00276 
00277     /** \returns the number of columns */
00278     inline Index rows() const { return m_rows.value(); }
00279 
00280     /** \returns the number of rows */
00281     inline Index cols() const { return m_coeffs.cols(); }
00282 
00283     /** \returns the number of super diagonals */
00284     inline Index supers() const { return m_supers.value(); }
00285 
00286     /** \returns the number of sub diagonals */
00287     inline Index subs() const { return m_subs.value(); }
00288 
00289     inline const CoefficientsType& coeffs() const { return m_coeffs; }
00290 
00291   protected:
00292 
00293     const CoefficientsType& m_coeffs;
00294     internal::variable_if_dynamic<Index, _Rows>   m_rows;
00295     internal::variable_if_dynamic<Index, _Supers> m_supers;
00296     internal::variable_if_dynamic<Index, _Subs>   m_subs;
00297 };
00298 
00299 /**
00300   * \class TridiagonalMatrix
00301   * \ingroup Core_Module
00302   *
00303   * \brief Represents a tridiagonal matrix with a compact banded storage
00304   *
00305   * \param _Scalar Numeric type, i.e. float, double, int
00306   * \param Size Number of rows and cols, or \b Dynamic
00307   * \param _Options Can be 0 or \b SelfAdjoint
00308   *
00309   * \sa class BandMatrix
00310   */
00311 template<typename Scalar, int Size, int Options>
00312 class TridiagonalMatrix : public BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>
00313 {
00314     typedef BandMatrix<Scalar,Size,Size,Options&SelfAdjoint?0:1,1,Options|RowMajor>  Base;
00315     typedef typename Base::Index Index;
00316   public:
00317     TridiagonalMatrix(Index size = Size) : Base(size,size,Options&SelfAdjoint?0:1,1) {}
00318 
00319     inline typename Base::template DiagonalIntReturnType<1>::Type super()
00320     { return Base::template diagonal<1>(); }
00321     inline const typename Base::template DiagonalIntReturnType<1>::Type super() const
00322     { return Base::template diagonal<1>(); }
00323     inline typename Base::template DiagonalIntReturnType<-1>::Type sub()
00324     { return Base::template diagonal<-1>(); }
00325     inline const typename Base::template DiagonalIntReturnType<-1>::Type sub() const
00326     { return Base::template diagonal<-1>(); }
00327   protected:
00328 };
00329 
00330 } // end namespace internal
00331 
00332 } // end namespace Eigen
00333 
00334 #endif // EIGEN_BANDMATRIX_H