From df6561a73fcaa11249d7fe62fa661be99c0afd67 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 15 Jul 2009 17:00:49 +0200 Subject: [PATCH] change the implementation of BandMatrix to follow the BLAS/LAPACK storage scheme --- Eigen/src/Core/BandMatrix.h | 119 +++++++++-------- Eigen/src/Core/util/ForwardDeclarations.h | 2 +- disabled/SkylineMatrix.h | 154 ++++++++++++++++++++++ test/bandmatrix.cpp | 2 +- 4 files changed, 220 insertions(+), 57 deletions(-) create mode 100644 disabled/SkylineMatrix.h diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index 03d17dac2..902b9a826 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -28,28 +28,28 @@ /** \nonstableyet * \class BandMatrix * - * \brief + * \brief * * \param * - * \sa + * \sa */ -template -struct ei_traits > +template +struct ei_traits > { typedef _Scalar Scalar; enum { CoeffReadCost = NumTraits::ReadCost, - RowsAtCompileTime = Size, - ColsAtCompileTime = Size, - MaxRowsAtCompileTime = Size, - MaxColsAtCompileTime = Size, + RowsAtCompileTime = Rows, + ColsAtCompileTime = Cols, + MaxRowsAtCompileTime = Rows, + MaxColsAtCompileTime = Cols, Flags = 0 }; }; -template -class BandMatrix : public MultiplierBase > +template +class BandMatrix : public MultiplierBase > { public: @@ -63,67 +63,79 @@ class BandMatrix : public MultiplierBase }; typedef typename ei_traits::Scalar Scalar; typedef Matrix PlainMatrixType; - + protected: enum { - DataSizeAtCompileTime = ((Size!=Dynamic) && (Supers!=Dynamic) && (Subs!=Dynamic)) - ? Size*(Supers+Subs+1) - (Supers*Supers+Subs*Subs)/2 - : Dynamic + DataRowsAtCompileTime = ((Supers!=Dynamic) && (Subs!=Dynamic)) + ? 1 + Supers + Subs + : Dynamic, + SizeAtCompileTime = EIGEN_ENUM_MIN(Rows,Cols) }; - typedef Matrix DataType; - + typedef Matrix DataType; + public: -// inline BandMatrix() { } - - inline BandMatrix(int size=Size, int supers=Supers, int subs=Subs) - : m_data(size*(supers+subs+1) - (supers*supers+subs*subs)/2), - m_size(size), m_supers(supers), m_subs(subs) + inline BandMatrix(int rows=Rows, int cols=Cols, int supers=Supers, int subs=Subs) + : m_data(1+supers+subs,cols), + m_rows(rows), m_supers(supers), m_subs(subs) { } - inline int rows() const { return m_size.value(); } - inline int cols() const { return m_size.value(); } + /** \returns the number of columns */ + inline int rows() const { return m_rows.value(); } + /** \returns the number of rows */ + inline int cols() const { return m_data.cols(); } + + /** \returns the number of super diagonals */ inline int supers() const { return m_supers.value(); } + + /** \returns the number of sub diagonals */ inline int subs() const { return m_subs.value(); } - inline VectorBlock diagonal() - { return VectorBlock(m_data,0,m_size.value()); } + /** \returns a vector expression of the main diagonal */ + inline Block diagonal() + { return Block(m_data,supers(),0,1,std::min(rows(),cols())); } - inline const VectorBlock diagonal() const - { return VectorBlock(m_data,0,m_size.value()); } + /** \returns a vector expression of the main diagonal (const version) */ + inline const Block diagonal() const + { return Block(m_data,supers(),0,1,std::min(rows(),cols())); } - template - VectorBlock - diagonal() + template struct DiagonalIntReturnType { + enum { + DiagonalSize = RowsAtCompileTime==Dynamic || ColsAtCompileTime==Dynamic + ? Dynamic + : Index<0 + ? EIGEN_ENUM_MIN(ColsAtCompileTime, RowsAtCompileTime + Index) + : EIGEN_ENUM_MIN(RowsAtCompileTime, ColsAtCompileTime - Index) + }; + typedef Block Type; + }; + + /** \returns a vector expression of the \a Index -th sub or super diagonal */ + template inline typename DiagonalIntReturnType::Type diagonal() { - return VectorBlock - (m_data,Index<0 ? subDiagIndex(-Index) : superDiagIndex(Index), m_size.value()-ei_abs(Index)); + return typename DiagonalIntReturnType::Type(m_data, supers()-Index, ei_abs(Index), 1, diagonalLength(Index)); } - template - const VectorBlock - diagonal() const + /** \returns a vector expression of the \a Index -th sub or super diagonal */ + template inline const typename DiagonalIntReturnType::Type diagonal() const { - return VectorBlock - (m_data,Index<0 ? subDiagIndex(-Index) : superDiagIndex(Index), m_size.value()-ei_abs(Index)); + return typename DiagonalIntReturnType::Type(m_data, supers()-Index, ei_abs(Index), 1, diagonalLength(Index)); } - inline VectorBlock diagonal(int index) + /** \returns a vector expression of the \a i -th sub or super diagonal */ + inline Block diagonal(int i) { - ei_assert((index<0 && -index<=subs()) || (index>=0 && index<=supers())); - return VectorBlock(m_data, - index<0 ? subDiagIndex(-index) : superDiagIndex(index), m_size.value()-ei_abs(index)); - } - const VectorBlock diagonal(int index) const - { - ei_assert((index<0 && -index<=subs()) || (index>=0 && index<=supers())); - return VectorBlock(m_data, - index<0 ? subDiagIndex(-index) : superDiagIndex(index), m_size.value()-ei_abs(index)); + ei_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers())); + return Block(m_data, supers()-i, ei_abs(i), 1, diagonalLength(i)); } -// inline VectorBlock subDiagonal() -// { return VectorBlock(m_data,0,m_size.value()); } + /** \returns a vector expression of the \a i -th sub or super diagonal */ + inline const Block diagonal(int i) const + { + ei_assert((i<0 && -i<=subs()) || (i>=0 && i<=supers())); + return Block(m_data, supers()-i, ei_abs(i), 1, diagonalLength(i)); + } PlainMatrixType toDense() const { @@ -139,14 +151,11 @@ class BandMatrix : public MultiplierBase protected: - inline int subDiagIndex(int i) const - { return m_size.value()*(m_supers.value()+i)-(ei_abs2(i-1) + ei_abs2(m_supers.value()))/2; } - - inline int superDiagIndex(int i) const - { return m_size.value()*i-ei_abs2(i-1)/2; } + inline int diagonalLength(int i) const + { return i<0 ? std::min(cols(),rows()+i) : std::min(rows(),cols()-i); } DataType m_data; - ei_int_if_dynamic m_size; + ei_int_if_dynamic m_rows; ei_int_if_dynamic m_supers; ei_int_if_dynamic m_subs; }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 9433a04cb..b4fbae28c 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -66,7 +66,7 @@ template class WithFormat; template struct CommaInitializer; template class ReturnByValue; -template class BandMatrix; +template class BandMatrix; template struct ei_product_mode; diff --git a/disabled/SkylineMatrix.h b/disabled/SkylineMatrix.h new file mode 100644 index 000000000..03d17dac2 --- /dev/null +++ b/disabled/SkylineMatrix.h @@ -0,0 +1,154 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_BANDMATRIX_H +#define EIGEN_BANDMATRIX_H + +/** \nonstableyet + * \class BandMatrix + * + * \brief + * + * \param + * + * \sa + */ +template +struct ei_traits > +{ + typedef _Scalar Scalar; + enum { + CoeffReadCost = NumTraits::ReadCost, + RowsAtCompileTime = Size, + ColsAtCompileTime = Size, + MaxRowsAtCompileTime = Size, + MaxColsAtCompileTime = Size, + Flags = 0 + }; +}; + +template +class BandMatrix : public MultiplierBase > +{ + public: + + enum { + Flags = ei_traits::Flags, + CoeffReadCost = ei_traits::CoeffReadCost, + RowsAtCompileTime = ei_traits::RowsAtCompileTime, + ColsAtCompileTime = ei_traits::ColsAtCompileTime, + MaxRowsAtCompileTime = ei_traits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = ei_traits::MaxColsAtCompileTime + }; + typedef typename ei_traits::Scalar Scalar; + typedef Matrix PlainMatrixType; + + protected: + enum { + DataSizeAtCompileTime = ((Size!=Dynamic) && (Supers!=Dynamic) && (Subs!=Dynamic)) + ? Size*(Supers+Subs+1) - (Supers*Supers+Subs*Subs)/2 + : Dynamic + }; + typedef Matrix DataType; + + public: + +// inline BandMatrix() { } + + inline BandMatrix(int size=Size, int supers=Supers, int subs=Subs) + : m_data(size*(supers+subs+1) - (supers*supers+subs*subs)/2), + m_size(size), m_supers(supers), m_subs(subs) + { } + + inline int rows() const { return m_size.value(); } + inline int cols() const { return m_size.value(); } + + inline int supers() const { return m_supers.value(); } + inline int subs() const { return m_subs.value(); } + + inline VectorBlock diagonal() + { return VectorBlock(m_data,0,m_size.value()); } + + inline const VectorBlock diagonal() const + { return VectorBlock(m_data,0,m_size.value()); } + + template + VectorBlock + diagonal() + { + return VectorBlock + (m_data,Index<0 ? subDiagIndex(-Index) : superDiagIndex(Index), m_size.value()-ei_abs(Index)); + } + + template + const VectorBlock + diagonal() const + { + return VectorBlock + (m_data,Index<0 ? subDiagIndex(-Index) : superDiagIndex(Index), m_size.value()-ei_abs(Index)); + } + + inline VectorBlock diagonal(int index) + { + ei_assert((index<0 && -index<=subs()) || (index>=0 && index<=supers())); + return VectorBlock(m_data, + index<0 ? subDiagIndex(-index) : superDiagIndex(index), m_size.value()-ei_abs(index)); + } + const VectorBlock diagonal(int index) const + { + ei_assert((index<0 && -index<=subs()) || (index>=0 && index<=supers())); + return VectorBlock(m_data, + index<0 ? subDiagIndex(-index) : superDiagIndex(index), m_size.value()-ei_abs(index)); + } + +// inline VectorBlock subDiagonal() +// { return VectorBlock(m_data,0,m_size.value()); } + + PlainMatrixType toDense() const + { + PlainMatrixType res(rows(),cols()); + res.setZero(); + res.diagonal() = diagonal(); + for (int i=1; i<=supers();++i) + res.diagonal(i) = diagonal(i); + for (int i=1; i<=subs();++i) + res.diagonal(-i) = diagonal(-i); + return res; + } + + protected: + + inline int subDiagIndex(int i) const + { return m_size.value()*(m_supers.value()+i)-(ei_abs2(i-1) + ei_abs2(m_supers.value()))/2; } + + inline int superDiagIndex(int i) const + { return m_size.value()*i-ei_abs2(i-1)/2; } + + DataType m_data; + ei_int_if_dynamic m_size; + ei_int_if_dynamic m_supers; + ei_int_if_dynamic m_subs; +}; + +#endif // EIGEN_BANDMATRIX_H diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp index 96ccda2cf..69ab0ed1d 100644 --- a/test/bandmatrix.cpp +++ b/test/bandmatrix.cpp @@ -55,7 +55,7 @@ template void bandmatrix(MatrixType& m) void test_bandmatrix() { for(int i = 0; i < g_repeat ; i++) { - BandMatrix m(6,3,2); + BandMatrix m(6,6,3,2); CALL_SUBTEST( bandmatrix(m) ); } }