From 587029a612df2a98a545702e90dc40feb7a62d3e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 14 Jul 2009 23:27:37 +0200 Subject: [PATCH] started an implementation of BandMatrix: at least the read/write access to the main/sub/super diagonals seems to work well. --- Eigen/Core | 1 + Eigen/src/Core/BandMatrix.h | 154 ++++++++++++++++++++++ Eigen/src/Core/util/ForwardDeclarations.h | 2 + test/CMakeLists.txt | 1 + test/bandmatrix.cpp | 61 +++++++++ 5 files changed, 219 insertions(+) create mode 100644 Eigen/src/Core/BandMatrix.h create mode 100644 test/bandmatrix.cpp diff --git a/Eigen/Core b/Eigen/Core index 2bec68f9c..ee9bfe325 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -185,6 +185,7 @@ namespace Eigen { #include "src/Core/SolveTriangular.h" #include "src/Core/products/SelfadjointRank2Update.h" #include "src/Core/products/TriangularMatrixVector.h" +#include "src/Core/BandMatrix.h" } // namespace Eigen diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h new file mode 100644 index 000000000..03d17dac2 --- /dev/null +++ b/Eigen/src/Core/BandMatrix.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/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 363972b60..9433a04cb 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -66,6 +66,8 @@ template class WithFormat; template struct CommaInitializer; template class ReturnByValue; +template class BandMatrix; + template struct ei_product_mode; template::value> struct ProductReturnType; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 886731d45..5e1adcbf4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -111,6 +111,7 @@ ei_add_test(array_replicate) ei_add_test(array_reverse) ei_add_test(triangular) ei_add_test(product_triangular) +ei_add_test(bandmatrix) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) ei_add_test(determinant) diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp new file mode 100644 index 000000000..96ccda2cf --- /dev/null +++ b/test/bandmatrix.cpp @@ -0,0 +1,61 @@ +// This file is triangularView of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-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 . + +#include "main.h" + +template void bandmatrix(MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix DenseMatrixType; + + int size = m.rows(); + + DenseMatrixType dm1(size,size); + dm1.setZero(); + + m.diagonal().setConstant(123); + dm1.diagonal().setConstant(123); + for (int i=1; i<=m.supers();++i) + { + m.diagonal(i).setConstant(i); + dm1.diagonal(i).setConstant(i); + } + for (int i=1; i<=m.subs();++i) + { + m.diagonal(-i).setConstant(-i); + dm1.diagonal(-i).setConstant(-i); + } + std::cerr << m.toDense() << "\n\n" << dm1 << "\n\n"; + VERIFY_IS_APPROX(dm1,m.toDense()); + +} + +void test_bandmatrix() +{ + for(int i = 0; i < g_repeat ; i++) { + BandMatrix m(6,3,2); + CALL_SUBTEST( bandmatrix(m) ); + } +}