BandMatrix.h 14 KB

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