From 289e0bd7a46010b5e8c5be20b3746055029a29f0 Mon Sep 17 00:00:00 2001 From: jenswehner Date: Sat, 20 Mar 2021 21:00:06 +0100 Subject: [PATCH 1/5] working on hermitian matrix --- Eigen/Core | 2 + Eigen/src/Core/HermitianBase.h | 244 +++++++++++++ Eigen/src/Core/HermitianMatrix.h | 590 +++++++++++++++++++++++++++++++ 3 files changed, 836 insertions(+) create mode 100644 Eigen/src/Core/HermitianBase.h create mode 100644 Eigen/src/Core/HermitianMatrix.h diff --git a/Eigen/Core b/Eigen/Core index 5921e15f9..19d7e141c 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -311,6 +311,8 @@ using std::ptrdiff_t; #include "src/Core/IndexedView.h" #include "src/Core/Reshaped.h" #include "src/Core/Transpose.h" +#include "src/Core/HermitianMatrix.h" +#include "src/Core/HermitianBase.h" #include "src/Core/DiagonalMatrix.h" #include "src/Core/Diagonal.h" #include "src/Core/DiagonalProduct.h" diff --git a/Eigen/src/Core/HermitianBase.h b/Eigen/src/Core/HermitianBase.h new file mode 100644 index 000000000..b71b0c6df --- /dev/null +++ b/Eigen/src/Core/HermitianBase.h @@ -0,0 +1,244 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2018 David Tellenbach +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_HERMITIANBASE_H +#define EIGEN_HERMITIANBASE_H +namespace Eigen { + +template +class HermitianBase : public EigenBase +{ + public: +#ifndef EIGEN_PARSED_BY_DOXYGEN + typedef typename internal::traits::NestedExpression NestedExpression; + typedef typename internal::traits::NestedExpression Nested; + typedef typename internal::traits::DenseType DenseType; + typedef typename NestedExpression::Scalar Scalar; + typedef typename internal::conditional::Flags&LvalueBit), const Scalar&, + typename internal::conditional::value, Scalar, const Scalar>::type>::type + CoeffReturnType; +#endif + + enum { + UpLo = internal::traits::UpLo, + RowsAtCompileTime = internal::traits::RowsAtCompileTime, + ColsAtCompileTime = internal::traits::ColsAtCompileTime, + MaxRowsAtCompileTime = internal::traits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime, + SizeAtCompileTime = RowsAtCompileTime * ColsAtCompileTime, + MaxSizeAtCompileTime = MaxRowsAtCompileTime * MaxColsAtCompileTime, + IsRowMajor = NestedExpression::IsRowMajor, + IsVectorAtCompileTime = 0 + }; + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Derived& operator=(const HermitianBase& other) + { + internal::call_assignment(derived(), other.derived()); + return derived(); + } + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Derived& _set(const HermitianBase& other) + { + internal::call_assignment(derived(), other.derived()); + return derived(); + } + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Derived& operator=(const DenseBase& other) + { + internal::call_assignment(derived(), other.derived()); + return derived(); + } + + + /** \returns constant reference to Derived */ + EIGEN_DEVICE_FUNC + inline const Derived& derived() const + { + return *static_cast(this); + } + + /** \returns reference to Derived */ + EIGEN_DEVICE_FUNC + inline Derived& derived() + { + return *static_cast(this); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Scalar coeff(Index row, Index col) const + { + eigen_internal_assert(row >= 0 && row < rows() + && col >= 0 && col < cols()); + return internal::evaluator >(derived()).coeff(row, col); + } + + /** \returns a reference to the coefficient at given the given row and + * column */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + HermitianMatrixCoeffReturnHelper operator()(Index row, Index col) + { + eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols()); + return coeffRef(row, col); + } + + /** Short version: don't use this function, use + * \link operator()(Index,Index) \endlink instead. + * + * Long version: this function is similar to + * \link operator()(Index,Index) \endlink, but without the assertion. + * Use this for limiting the performance cost of debugging code when doing + * repeated coefficient access. Only use this when it is guaranteed that the + * parameters \a row and \a col are in range. + * + * If EIGEN_INTERNAL_DEBUGGING is defined, an assertion will be made, making this + * function equivalent to \link operator()(Index,Index) \endlink. + * + * \sa operator()(Index,Index), coeff(Index, Index) const + */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + HermitianMatrixCoeffReturnHelper coeffRef(Index row, Index col) + { + eigen_internal_assert(row >= 0 && row < rows() && col >= 0 && col < cols()); + return internal::evaluator >(derived()).coeffRef(row, col); + } + + /** \returns a const reference to the nested expression */ + EIGEN_DEVICE_FUNC + const NestedExpression& nestedExpression() const + { + return derived().nestedExpression(); + } + + /** \returns a reference to the nested expression */ + EIGEN_DEVICE_FUNC + NestedExpression& nestedExpression() + { + return derived().nestedExpression(); + } + + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + EIGEN_DEVICE_FUNC + Index cols() + { + return derived().cols(); + } + + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + EIGEN_DEVICE_FUNC + Index cols() const + { + return derived().cols(); + } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + EIGEN_DEVICE_FUNC + Index rows() + { + return derived().rows(); + } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + EIGEN_DEVICE_FUNC + Index rows() const + { + return derived().rows(); + } + + EIGEN_DEVICE_FUNC + DenseType toDenseMatrix() const { return derived(); } + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Derived& operator+=(const HermitianBase& other) + { + call_assignment(nestedExpression(), other.nestedExpression(), + internal::add_assign_op()); + return derived(); + } + + template + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + Derived& operator-=(const HermitianBase& other) + { + call_assignment(nestedExpression(), other.nestedExpression(), + internal::sub_assign_op()); + return derived(); + } + + /** \returns CwiseBinaryOp that represents an abstract sum expression */ + template + EIGEN_STRONG_INLINE const + CwiseBinaryOp::Scalar, + typename internal::traits::Scalar >, const Derived, const OtherDerived> + operator+(const HermitianBase &other) const + { + return CwiseBinaryOp::Scalar, + typename internal::traits::Scalar >, const Derived, const OtherDerived> + (derived(), other.derived()); + } + + /** \returns CwiseBinaryOp that represents an abstract difference expression */ + template + EIGEN_STRONG_INLINE const + CwiseBinaryOp::Scalar, + typename internal::traits::Scalar >, const Derived, const OtherDerived> + operator-(const HermitianBase &other) const + { + return CwiseBinaryOp::Scalar, + typename internal::traits::Scalar >, const Derived, const OtherDerived> + (derived(), other.derived()); + } + + + // hermitian * hermitian + template + EIGEN_STRONG_INLINE const + Product + operator*(const HermitianBase& rhs) const + { + return Product(derived(), rhs.derived()); + } + + // hermitian * dense + template + EIGEN_STRONG_INLINE const + Product + operator*(const MatrixBase &other) const + { + return Product(derived(), other.derived()); + } + + typedef typename internal::add_const_on_value_type::type>::type EvalReturnType; + + /** \returns the matrix or vector obtained by evaluating this expression. + * + * Notice that in the case of a plain matrix or vector (not an expression) this function just returns + * a const reference, in order to avoid a useless copy. + * + * \warning Be careful with eval() and the auto C++ keyword, as detailed in this \link TopicPitfalls_auto_keyword page \endlink. + */ + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE EvalReturnType eval() const + { + // Even though MSVC does not honor strong inlining when the return type + // is a dynamic matrix, we desperately need strong inlining for fixed + // size types on MSVC. + return typename internal::eval::type(derived()); + } +}; + +} // namespace Eigen + +#endif // EIGEN_HERMITIANBASE_H diff --git a/Eigen/src/Core/HermitianMatrix.h b/Eigen/src/Core/HermitianMatrix.h new file mode 100644 index 000000000..475eb8734 --- /dev/null +++ b/Eigen/src/Core/HermitianMatrix.h @@ -0,0 +1,590 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2018 David Tellenbach +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_HERMITIANMATRIX_H +#define EIGEN_HERMITIANMATRIX_H + +namespace Eigen { +namespace internal { + +#ifndef EIGEN_PARSED_BY_DOXYGEN +// Dimension: even +template +struct traits > + : traits > +{ + typedef Matrix<_Scalar, _Dimension + 1, _Dimension / 2, _Storage, _MaxDimension + 1, _MaxDimension / 2> NestedExpression; + typedef _Scalar Scalar; + typedef Matrix<_Scalar, _Dimension, _Dimension, _Storage, _MaxDimension, _MaxDimension> DenseType; + typedef HermitianMatrix<_Scalar, Dynamic, _UpLo, _Storage, _MaxDimension> PlainObject; + typedef HermitianShape StorageKind; + typedef HermitianXpr XprKind; + enum { + Flags = LvalueBit | NoPreferredStorageOrderBit | NestByRefBit, + UpLo = _UpLo, + RowsAtCompileTime = _Dimension, + ColsAtCompileTime = _Dimension, + MaxRowsAtCompileTime = _MaxDimension, + MaxColsAtCompileTime = _MaxDimension + }; +}; + +// Dimension: odd +template +struct traits > + : traits > +{ + typedef Matrix<_Scalar, _Dimension, (_Dimension + 1) / 2, _Storage, _MaxDimension, (_MaxDimension + 1) / 2> NestedExpression; + typedef _Scalar Scalar; + typedef Matrix<_Scalar, _Dimension, _Dimension, _Storage, _MaxDimension, _MaxDimension> DenseType; + typedef HermitianMatrix<_Scalar, Dynamic, _UpLo, _Storage, _MaxDimension> PlainObject; + typedef HermitianShape StorageKind; + typedef HermitianXpr XprKind; + + enum { + Flags = LvalueBit | NoPreferredStorageOrderBit | NestByRefBit, + UpLo = _UpLo, + RowsAtCompileTime = _Dimension, + ColsAtCompileTime = _Dimension, + MaxRowsAtCompileTime = _MaxDimension, + MaxColsAtCompileTime = _MaxDimension + }; +}; + +// Dimension::Eigen::Dynamic +template +struct traits > + : traits > +{ + typedef Matrix<_Scalar, Dynamic, Dynamic, _Storage, _MaxDimension, _MaxDimension> NestedExpression; + typedef Matrix<_Scalar, Dynamic, Dynamic, _Storage, _MaxDimension, _MaxDimension> DenseType; + typedef HermitianMatrix<_Scalar, Dynamic, _UpLo, _Storage, _MaxDimension> PlainObject; + typedef _Scalar Scalar; + typedef HermitianShape StorageKind; + typedef HermitianXpr XprKind; + + enum { + Flags = LvalueBit | NoPreferredStorageOrderBit | NestByRefBit, + UpLo = _UpLo, + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic, + MaxRowsAtCompileTime = _MaxDimension, + MaxColsAtCompileTime = _MaxDimension + }; +}; +#endif + +template +struct hermitian_prod_impl { + + typedef typename internal::traits::Scalar Scalar; + + static void run(DstXprType& dst, const Lhs& lhs, const Rhs& rhs) + { + const Index dimension = dst.cols(); + + const Matrix& tmp = lhs.nestedExpression().block(1, 0, dimension / 2, dimension / 2).template selfadjointView(); + const Matrix& tmp2 = lhs.nestedExpression().block(0, 0, dimension / 2, dimension / 2).template selfadjointView(); + #pragma omp parallel sections + { + #pragma omp section + { + dst.nestedExpression().block(1, 0, dimension / 2, dimension / 2).template triangularView() = tmp * rhs.nestedExpression().block(1, 0, dimension / 2, dimension / 2).template selfadjointView() + + lhs.nestedExpression().block(dimension/2 + 1, 0, dimension / 2, dimension / 2).transpose() + * rhs.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2); + } + #pragma omp section + { + dst.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2) = lhs.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2) + * rhs.nestedExpression().block(1, 0, dimension/2, dimension/2).template selfadjointView() + + lhs.nestedExpression().block(0, 0, dimension/2, dimension/2).template selfadjointView() + * rhs.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2); + } + #pragma omp section + { + dst.nestedExpression().block(0, 0, dimension/2, dimension/2).template triangularView() = lhs.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2) + * rhs.nestedExpression().block(dimension/2 + 1, 0, dimension/2, dimension/2).transpose() + + tmp2 * rhs.nestedExpression().block(0, 0, dimension/2, dimension/2).template selfadjointView(); + } + } + + } +}; +} // namespace internal + +/** \class Hermitian + * \ingroup Core_Module + * + * \brief Represents a hermitian matrix with its storage + * + * \param _Scalar the type of coefficients + * \param _Dimension the dimension of the matrix, or Eigen::Dynamic + * \param _UpLo can take the values Eigen::Upper to store only the upper triagonal part of the matrix or Eigen::Lower + * to store only the lower triangular part respectively. This parameter is optional and defaults to + * Eigen::Upper. + * \param _Storage can take the values Eigen::RowMajor for row-major storage of the underlying Eigen::Matrix or + * Eigen::ColMajor for column-major storage respectively. This parameter is optional and defaults to + * Eigen::ColMajor. + * \param _MaxDimension the dimension of the matrix or Eigen::Dynamic. This parameter is option and defaults to + * _Dimension. + */ +template +class HermitianMatrix + : public HermitianBase > +{ + public: +#ifndef EIGEN_PARSED_BY_DOXYGEN + typedef typename internal::traits::NestedExpression NestedExpression; + typedef typename internal::traits::StorageKind StorageKind; + typedef typename internal::traits::StorageIndex StorageIndex; + typedef HermitianBase > Base; +#endif + + enum { + Dimension = _Dimension, + UpLo = _UpLo + }; + + /** Default constructor without initialization */ + EIGEN_DEVICE_FUNC + inline HermitianMatrix() + : m_dimension(_Dimension == Dynamic ? 0 : _Dimension), + m_odd(_Dimension == Dynamic ? 0 : _Dimension % 2) + {} + + /** Constructs a hermitian matrix with given dimension */ + EIGEN_DEVICE_FUNC + explicit inline HermitianMatrix(Index dim) + : m_odd(dim % 2), m_dimension(dim) + { + if (dim % 2) { + m_nested.resize(dim, (dim + 1) / 2); + } else { + m_nested.resize(dim + 1, dim / 2); + } + } + + /** Constructs a HermitianMatrix by passing a dense matrix */ + EIGEN_DEVICE_FUNC + template + explicit inline HermitianMatrix(const MatrixBase& other) + : m_dimension(other.cols()), m_odd(other.cols() % 2) + { + // Replace with generic template argument check + + // _UpLo must not differ from Eigen::Lower and Eigen::Upper + EIGEN_STATIC_ASSERT(_UpLo == Upper || _UpLo == Lower, + HERMITIANMATRIX_ACCEPTS_UPPER_AND_LOWER_MODE_ONLY); + + // other is quadratic matrix and its size matches the size of this + EIGEN_STATIC_ASSERT(OtherDerived::RowsAtCompileTime == _Dimension + && OtherDerived::RowsAtCompileTime == OtherDerived::ColsAtCompileTime, + YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES); + + if (_Dimension == Dynamic) { + if (m_odd) { + m_nested.resize(m_dimension, (m_dimension + 1) / 2); + } else { + m_nested.resize(m_dimension + 1, m_dimension / 2); + } + } + + // Dimension: Even + // UpLo: Upper + if (_UpLo == Upper && !m_odd) { + for (Index row = 0; row < m_dimension; ++row) { + for (Index col = row; col < m_dimension; ++col) { + if (col < (m_dimension / 2)) { + m_nested(m_dimension / 2 + col + 1, row) = other(row, col); + } else { + m_nested(row, col - m_dimension / 2) = other(row, col); + } + } + } + return; + } + + // Dimension: Even + // UpLo: Lower + if (_UpLo == Lower && !m_odd) { + for (Index row = 0; row < m_dimension; ++row) { + for (Index col = 0; col <= row; ++col) { + if (col < m_dimension / 2) { + m_nested(row +1, col) = other(row, col); + } else { + m_nested(col - m_dimension / 2, row - m_dimension / 2) = other(row, col); + } + } + } + return; + } + + // Dimension: odd + // UpLo: Upper + if (_UpLo == Upper && m_odd) { + for (Index row = 0; row < m_dimension; ++row) { + for (Index col = row; col < m_dimension; ++col) { + if (col < ((m_dimension - 1) / 2)) { + m_nested((m_dimension - 1)/2 + col + 1, row) = other(row, col); + } else { + m_nested(row, col - (m_dimension - 1) / 2) = other(row, col); + } + } + } + return; + } + + // Dimension: Odd + // UpLo: Lower + if (_UpLo == Lower && m_odd) { + for (Index row = 0; row < m_dimension; ++row) { + for (Index col = 0; col <= row; ++col) { + if (col < (other.cols() + 1) / 2) { + m_nested(row, col) = other(row, col); + } else { + m_nested(col - (m_dimension - 1) / 2 - 1, row - (m_dimension + 1) / 2 + 1) = other(row, col); + } + } + } + return; + } + } + + /** Copy constructor */ + template + EIGEN_DEVICE_FUNC + inline HermitianMatrix(const HermitianBase& other) + : m_nested(other.nestedExpression()), m_dimension(other.cols()), m_odd(other.cols() % 2) {} + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */ + explicit inline HermitianMatrix(const HermitianMatrix& other) + : m_nested(other.nestedExpression()), m_dimension(other.m_dimension), + m_odd(other.m_dimension % 2) {} +#endif + + template + explicit inline HermitianMatrix(const CwiseBinaryOp::Scalar >, + const HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, const OtherDerived >& other) + : m_dimension(other.lhs().cols()), m_odd(other.lhs().cols() % 2), + m_nested(other.lhs().nestedExpression() + other.rhs().nestedExpression()){} + + template + explicit inline HermitianMatrix(const CwiseBinaryOp::Scalar >, + const HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, const OtherDerived >& other) + : m_dimension(other.lhs().cols()), m_odd(other.lhs().cols() % 2) + { + m_nested = other.lhs().nestedExpression() - other.rhs().nestedExpression(); + } + + /** Copy assignment constructor */ + template + EIGEN_DEVICE_FUNC + HermitianMatrix& operator=(const HermitianBase& other) + { + return Base::operator=(other); + } + + /** Copy assignment constructor */ + template + EIGEN_DEVICE_FUNC + HermitianMatrix& operator=(const DenseBase& other) + { + return Base::operator=(other); + } + + template + EIGEN_DEVICE_FUNC + HermitianMatrix& + operator=(const CwiseBinaryOp::Scalar>, + const HermitianMatrix, const OtherDerived >& other) + { + m_dimension = other.lhs().cols(); + m_odd = m_dimension % 2; + + // Just add and assign the nested expressions of the left and right hand + // side of the abstract expression + m_nested = other.lhs().nestedExpression() + other.rhs().nestedExpression(); + return *this; + } + + template + EIGEN_DEVICE_FUNC + HermitianMatrix& + operator=(const CwiseBinaryOp::Scalar>, + const HermitianMatrix, const OtherDerived >& other) + { + m_dimension = other.lhs().cols(); + m_odd = m_dimension % 2; + + // Just add and assign the nested expressions of the left and right hand + // side of the abstract expression + m_nested = other.lhs().nestedExpression() - other.rhs().nestedExpression(); + return *this; + } + + template + EIGEN_DEVICE_FUNC + HermitianMatrix& + operator=(const Product, + const OtherDerived, LazyProduct>& product) { + m_dimension = product.lhs().cols(); + m_odd = m_dimension % 2; + + if (_Dimension == Dynamic) { + if (m_odd) { + m_nested.resize(m_dimension, (m_dimension + 1) / 2); + } else { + m_nested.resize(m_dimension + 1, m_dimension / 2); + } + } + // TODO:: Distinguish between dynamic and fixed dimension + // -> Dynamic or fixed dimension for temporaries + // -> Dynamic or fixed form of block expression + + // Even dimension + if (!m_odd) { + if (_UpLo == Lower) { + const typename HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>::NestedExpression& nestedLhs = product.lhs().nestedExpression(); + const typename OtherDerived::NestedExpression& nestedRhs = product.rhs().nestedExpression(); + + const Matrix<_Scalar, Dynamic, Dynamic>& tmp = nestedLhs.block(1, 0, m_dimension/2, m_dimension/2).template selfadjointView(); + const Matrix<_Scalar, Dynamic, Dynamic>& tmp2 = nestedLhs.block(0, 0, m_dimension/2, m_dimension/2).template selfadjointView(); + #pragma omp parallel sections + { + #pragma omp section + { + this->nestedExpression().block(1, 0, m_dimension/2, m_dimension/2).template triangularView() = tmp * nestedRhs.block(1, 0, m_dimension/2, m_dimension/2).template selfadjointView() + + nestedLhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2).transpose() + * nestedRhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2); + } + #pragma omp section + { + this->nestedExpression().block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2) = nestedLhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2) + * nestedRhs.block(1, 0, m_dimension/2, m_dimension/2).template selfadjointView() + + nestedLhs.block(0, 0, m_dimension/2, m_dimension/2).template selfadjointView() + * nestedRhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2); + } + #pragma omp section + { + this->nestedExpression().block(0, 0, m_dimension/2, m_dimension/2).template triangularView() = nestedLhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2) + * nestedRhs.block(m_dimension/2 + 1, 0, m_dimension/2, m_dimension/2).transpose() + + tmp2 + * nestedRhs.block(0, 0, m_dimension/2, m_dimension/2).template selfadjointView(); + } + } + + } else if (_UpLo == Upper) { + + return *this; + } + } + return *this; + + } + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** copy assignment constructor. Prevent a default copy constructor from + * hiding the other templated constructor */ + // EIGEN_DEVICE_FUNC + // HermitianMatrix& operator=(const HermitianMatrix& other) + // { + // m_dimension = other.cols(); + // m_odd = m_dimension % 2; + // m_nested = other.nestedExpression(); + // return *this; + // } +#endif + + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + EIGEN_DEVICE_FUNC + Index cols() + { + return m_dimension; + } + + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + EIGEN_DEVICE_FUNC + Index cols() const + { + return m_dimension; + } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + EIGEN_DEVICE_FUNC + Index rows() + { + return m_dimension; + } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + EIGEN_DEVICE_FUNC + Index rows() const + { + return m_dimension; + } + + /** \returns a const reference to the nested expression */ + EIGEN_DEVICE_FUNC + const NestedExpression& nestedExpression() const + { + return m_nested; + } + + /** \returns a reference to the nested expression */ + EIGEN_DEVICE_FUNC + NestedExpression& nestedExpression() { + return m_nested; + } + + EIGEN_DEVICE_FUNC + inline void resize(Index dim) + { + if (dim % 2) { + m_nested.resize(dim, (dim + 1) / 2); + } else { + m_nested.resize(dim + 1, dim / 2); + } + m_dimension = dim; + } + + /** Sets all coefficients to zero. */ + EIGEN_DEVICE_FUNC + inline void setZero() + { + m_nested.setZero(); + } + + EIGEN_DEVICE_FUNC + inline static HermitianMatrix Random() + { + HermitianMatrix ret; + ret.nestedExpression() = NestedExpression::Random(); + return ret; + } + + EIGEN_DEVICE_FUNC + inline static HermitianMatrix Random(Index dim) + { + HermitianMatrix ret; + ret.nestedExpression() = NestedExpression::Random(dim, dim); + } + + private: + NestedExpression m_nested; + bool m_odd; + Index m_dimension; +}; + +namespace internal { + +template<> struct storage_kind_to_shape { typedef HermitianShape Shape; }; + +struct Hermitian2Dense {}; +struct Hermitian2Hermitian {}; +struct Dense2Hermitian {}; + +template<> struct AssignmentKind { typedef Hermitian2Hermitian Kind; }; +template<> struct AssignmentKind { typedef Hermitian2Dense Kind; }; +template<> struct AssignmentKind { typedef Dense2Hermitian Kind; }; + +// Hermitian to Hermitian assignment +template< typename DstXprType, typename SrcXprType, typename Functor> +struct Assignment +{ + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + Index dstDim = src.rows(); + if(dst.rows() != dstDim) { + dst.resize(dstDim); + } + + eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); + dst.nestedExpression() = src.nestedExpression(); + } + + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) + { dst.nestedExpression() += src.nestedExpression(); } + + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) + { dst.nestedExpression() -= src.nestedExpression(); } +}; + +// Hermitian matrix to Dense assignment +template< typename DstXprType, typename SrcXprType, typename Functor> +struct Assignment +{ + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + const Index dim = src.cols(); + if (!DstXprType::IsRowMajor) { + for (Index col = 0; col < dim; ++col) { + for (Index row = col; row < dim; ++row) { + typename SrcXprType::Scalar val = src.coeff(row, col); + dst.coeffRef(row, col) = val; + dst.coeffRef(col, row) = val; + } + } + } else { + for (Index row = 0; row < dim; ++row) { + for (Index col = 0; col < row; ++col) { + typename SrcXprType::Scalar val = src.coeff(row, col); + dst.coeffRef(row, col) = val; + dst.coeffRef(col, row) = val; + } + } + } + } + + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) + { dst.nestedExpression() += src.nestedExpression(); } + + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) + { dst.nestedExpression() -= src.nestedExpression(); } +}; + +// Hermitian product to Dense assignment +template< typename DstXprType, typename _Scalar, int _Dimension, unsigned int _UpLo, int _Storage, int _MaxDimension, typename Functor> +struct Assignment, HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, LazyProduct> , Functor, Dense2Dense> +{ + typedef Product, HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, LazyProduct > SrcXprType; + + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + dst = src.lhs().toDenseMatrix() * src.rhs().toDenseMatrix(); + } + +}; + +// Hermitian product to hermitian assignment +template< typename DstXprType, typename _Scalar, int _Dimension, unsigned int _UpLo, int _Storage, int _MaxDimension, typename Functor> +struct Assignment, HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, LazyProduct> , Functor, Dense2Hermitian> +{ + typedef Product, HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension>, LazyProduct > SrcXprType; + + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + hermitian_prod_impl, HermitianMatrix<_Scalar, _Dimension, _UpLo, _Storage, _MaxDimension> >::run(dst, src.lhs(), src.rhs()); + } +}; + +} // namespace internal + +// dense * hermitian +template +template +EIGEN_DEVICE_FUNC inline const Product +MatrixBase::operator*(const HermitianBase& hermitian) const +{ + return Product(derived(), hermitian.derived()); +} + +} // namespace Eigen + +#endif // EIGEN_HERMITIANMATRIX_H -- GitLab From c32ad9e6173cde363c32d90dbb1a92fbe136bad5 Mon Sep 17 00:00:00 2001 From: jenswehner Date: Wed, 29 Sep 2021 18:05:30 +0200 Subject: [PATCH 2/5] moved to unsupported --- unsupported/Eigen/HermitianMatrix | 36 +++++++++++++++++++ .../src/HermitianMatrix}/HermitianBase.h | 0 .../src/HermitianMatrix}/HermitianMatrix.h | 2 +- 3 files changed, 37 insertions(+), 1 deletion(-) create mode 100644 unsupported/Eigen/HermitianMatrix rename {Eigen/src/Core => unsupported/Eigen/src/HermitianMatrix}/HermitianBase.h (100%) rename {Eigen/src/Core => unsupported/Eigen/src/HermitianMatrix}/HermitianMatrix.h (99%) diff --git a/unsupported/Eigen/HermitianMatrix b/unsupported/Eigen/HermitianMatrix new file mode 100644 index 000000000..8eb18635a --- /dev/null +++ b/unsupported/Eigen/HermitianMatrix @@ -0,0 +1,36 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2021 Jens Wehner +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_HERMITIANMATRIX_MODULE_H +#define EIGEN_HERMITIANMATRIX_MODULE_H + +#include "../../Eigen/Core" + +#include "../../Eigen/src/Core/util/DisableStupidWarnings.h" + + +/** + * \defgroup HermitianMatrix_Module HermitianMatrix Module + * + * This module contains an implementation of the Packed Format for symmetric and hermitian matrices + * + * \code + * #include + * \endcode + */ + + +#include "src/HermitianMatrix/HermitianBase.h" + +#include "src/HermitianMatrix/HermitianMatrix.h" + + +#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h" + +#endif // EIGEN_HERMITIANMATRIX_MODULE_H diff --git a/Eigen/src/Core/HermitianBase.h b/unsupported/Eigen/src/HermitianMatrix/HermitianBase.h similarity index 100% rename from Eigen/src/Core/HermitianBase.h rename to unsupported/Eigen/src/HermitianMatrix/HermitianBase.h diff --git a/Eigen/src/Core/HermitianMatrix.h b/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h similarity index 99% rename from Eigen/src/Core/HermitianMatrix.h rename to unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h index 475eb8734..c29d30526 100644 --- a/Eigen/src/Core/HermitianMatrix.h +++ b/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h @@ -119,7 +119,7 @@ struct hermitian_prod_impl { } // namespace internal /** \class Hermitian - * \ingroup Core_Module + * \ingroup HermitianMatrix_Module * * \brief Represents a hermitian matrix with its storage * -- GitLab From efc064beed6efc36613804f6179b2b10e96ddaa2 Mon Sep 17 00:00:00 2001 From: jenswehner Date: Wed, 29 Sep 2021 18:06:17 +0200 Subject: [PATCH 3/5] removed include from core --- Eigen/Core | 2 -- 1 file changed, 2 deletions(-) diff --git a/Eigen/Core b/Eigen/Core index 52dc52ecd..264398139 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -313,8 +313,6 @@ using std::ptrdiff_t; #include "src/Core/IndexedView.h" #include "src/Core/Reshaped.h" #include "src/Core/Transpose.h" -#include "src/Core/HermitianMatrix.h" -#include "src/Core/HermitianBase.h" #include "src/Core/DiagonalMatrix.h" #include "src/Core/Diagonal.h" #include "src/Core/DiagonalProduct.h" -- GitLab From 0282a4f62bbbed7ce66c5bc12b15f56d5424cd26 Mon Sep 17 00:00:00 2001 From: jenswehner Date: Wed, 29 Sep 2021 21:17:12 +0200 Subject: [PATCH 4/5] more stuff for hermitian matrix --- .../Eigen/src/HermitianMatrix/Evaluators.h | 408 ++++++++++++++++++ unsupported/test/hermitianmatrix.cpp | 35 ++ 2 files changed, 443 insertions(+) create mode 100644 unsupported/Eigen/src/HermitianMatrix/Evaluators.h create mode 100644 unsupported/test/hermitianmatrix.cpp diff --git a/unsupported/Eigen/src/HermitianMatrix/Evaluators.h b/unsupported/Eigen/src/HermitianMatrix/Evaluators.h new file mode 100644 index 000000000..53d7cfcf5 --- /dev/null +++ b/unsupported/Eigen/src/HermitianMatrix/Evaluators.h @@ -0,0 +1,408 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Benoit Jacob +// Copyright (C) 2011-2014 Gael Guennebaud +// Copyright (C) 2011-2012 Jitse Niesen +// Copyright (C) 2018 David Tellenbach +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + + +#ifndef EIGEN_HERMITIAN_EVALUATORS_H +#define EIGEN_HERMITIAN_EVALUATORS_H + + +namespace Eigen { +// NOTE: struct HermitianMatrixCoeffReturnHelper cannot be member of namespace +// internal since it is the return type of HermitianBase::coeffRef() + +// Scalar is no complex type, i.e., a type for that +// internal::is_complex::value is false +template +struct HermitianMatrixCoeffReturnHelper { + + HermitianMatrixCoeffReturnHelper(int row, int col, const NestedExpression* nested, const bool isStoredValue) + : m_isStoredValue(isStoredValue), m_nested(nested), m_row(row), m_col(col) {} + + operator Scalar() const { + return m_nested->coeff(m_row, m_col); + } + + HermitianMatrixCoeffReturnHelper& operator=(const Scalar& value) { + const_cast(m_nested)->coeffRef(m_row, m_col) = value; + return *this; + } + + template + OtherScalar operator+(const HermitianMatrixCoeffReturnHelper& rhs) const { + return static_cast(*this) + static_cast(rhs); + } + + template + OtherScalar operator-(const HermitianMatrixCoeffReturnHelper& rhs) const { + return static_cast(*this) - static_cast(rhs); + } + + template + OtherScalar operator*(const HermitianMatrixCoeffReturnHelper& rhs) const { + return static_cast(*this) * static_cast(rhs); + } + + template + OtherScalar operator/(const HermitianMatrixCoeffReturnHelper& rhs) const { + return static_cast(*this) / static_cast(rhs); + } + + template + HermitianMatrixCoeffReturnHelper& operator+=(const OtherScalar& other) { + const_cast(m_nested)->coeffRef(m_row, m_col) += other; + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator-=(const OtherScalar& other) { + const_cast(m_nested)->coeffRef(m_row, m_col) -= other; + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator*=(const OtherScalar& other) { + const_cast(m_nested)->coeffRef(m_row, m_col) *= other; + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator/=(const OtherScalar& other) { + const_cast(m_nested)->coeffRef(m_row, m_col) /= other; + return *this; + } + +private: + const bool m_isStoredValue; + const NestedExpression* m_nested; + const int m_row; + const int m_col; +}; + +// Scalar is a complex type, i.e., a type for that +// internal::is_complex::value is true +template +struct HermitianMatrixCoeffReturnHelper : public Scalar { + + HermitianMatrixCoeffReturnHelper(int row, int col, const NestedExpression* nested, const bool isStoredValue) + : Scalar((isStoredValue)?(nested->coeff(row, col)):(std::conj(nested->coeff(row, col)))), + m_isStoredValue(isStoredValue), m_nested(nested), m_row(row), m_col(col) {} + + HermitianMatrixCoeffReturnHelper& operator=(const Scalar& value) { + if (m_isStoredValue) { + const_cast(m_nested)->coeffRef(m_row, m_col) = value; + } else { + Scalar tmp = std::conj(value); + const_cast(m_nested)->coeffRef(m_row, m_col) = tmp; + } + return *this; + } + + operator Scalar() { + if (m_isStoredValue) { + return m_nested->coeffRef(m_row, m_col); + } else { + return std::conj(m_nested->coeffRef(m_row, m_col)); + } + } + + template + Scalar operator+(const HermitianMatrixCoeffReturnHelper& other) const { + return *this + static_cast(other); + } + + template + Scalar operator-(const HermitianMatrixCoeffReturnHelper& other) const { + return *this - static_cast(other); + } + + template + Scalar operator*(const HermitianMatrixCoeffReturnHelper& other) const { + return *this * static_cast(other); + } + + template + Scalar operator/(const HermitianMatrixCoeffReturnHelper& other) const { + return *this / static_cast(other); + } + + template + bool operator==(const HermitianMatrixCoeffReturnHelper& other) const { + return *this == static_cast(other); + } + + template + bool operator!=(const HermitianMatrixCoeffReturnHelper& other) const { + return *this != static_cast(other); + } + + template + HermitianMatrixCoeffReturnHelper& operator+=(const OtherScalar& other) { + if (m_isStoredValue) { + const_cast(m_nested)->coeffRef(m_row, m_col) += other; + } else { + Scalar tmp = std::conj(other); + const_cast(m_nested)->coeffRef(m_row, m_col) += tmp; + } + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator-=(const OtherScalar& other) { + if (m_isStoredValue) { + const_cast(m_nested)->coeffRef(m_row, m_col) -= other; + } else { + Scalar tmp = std::conj(other); + const_cast(m_nested)->coeffRef(m_row, m_col) -= tmp; + } + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator*=(const OtherScalar& other) { + if (m_isStoredValue) { + const_cast(m_nested)->coeffRef(m_row, m_col) *= other; + } else { + Scalar tmp = std::conj(other); + const_cast(m_nested)->coeffRef(m_row, m_col) *= tmp; + } + return *this; + } + + template + HermitianMatrixCoeffReturnHelper& operator/=(const OtherScalar& other) { + if (m_isStoredValue) { + const_cast(m_nested)->coeffRef(m_row, m_col) /= other; + } else { + Scalar tmp = std::conj(other); + const_cast(m_nested)->coeffRef(m_row, m_col) /= tmp; + } + return *this; + } + +private: + const bool m_isStoredValue; + const NestedExpression* m_nested; + const int m_row; + const int m_col; +}; + +namespace internal { + +template +struct evaluator > : evaluator_base > +{ + typedef HermitianBase HermitianType; + typedef typename HermitianType::Scalar Scalar; + typedef typename HermitianType::CoeffReturnType CoeffReturnType; + typedef typename internal::traits::NestedExpression NestedExpression; + + enum { + UpLo = HermitianType::UpLo, + RowsAtCompileTime = HermitianType::RowsAtCompileTime, + ColsAtCompileTime = HermitianType::ColsAtCompileTime, + SizeAtCompileTime = HermitianType::SizeAtCompileTime, + Flags = (unsigned int)(NoPreferredStorageOrderBit & LvalueBit) + }; + +#if EIGEN_HAS_CXX11 + EIGEN_DEVICE_FUNC + evaluator() : m_nested(nullptr) {} +#else + EIGEN_DEVICE_FUNC + evaluator() : m_nested(0) {} +#endif + + EIGEN_DEVICE_FUNC evaluator(const HermitianType& m) + : m_nested(&m.nestedExpression()), m_dimension(m.cols()) {} + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + HermitianMatrixCoeffReturnHelper coeff(Index row, Index col) const + { + // Even dimension and upper triangular storage + if (UpLo == Upper && !(m_dimension % 2)) { + if (row <= col) { + if (col < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(m_dimension / 2 + col + 1, row, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(row, col - m_dimension / 2, m_nested, true); + } + } else { + if (col <= row) { + if (row < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(m_dimension / 2 + row + 1, col, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(col, row - m_dimension / 2, m_nested, false); + } + } + } + } + + // Even dimension and lower triangular storage + if (UpLo == Lower && !(m_dimension % 2)) { + if (row < col) { + if (row < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(col + 1, row, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(row - m_dimension / 2, col - m_dimension / 2, + m_nested, false); + } + } else { + if (col < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(row + 1, col, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(col - m_dimension / 2, row - m_dimension / 2, + m_nested, true); + } + } + } + + // Odd dimension and upper triangular storage + if (UpLo == Upper && (m_dimension % 2)) { + if (row <= col) { + if (col < (m_dimension - 1) / 2) { + return HermitianMatrixCoeffReturnHelper((m_dimension - 1) / 2 + col + 1, row, + m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(row, col - (m_dimension - 1) / 2, + m_nested, true); + } + } else { + if (row < (m_dimension - 1) / 2) { + return HermitianMatrixCoeffReturnHelper((m_dimension - 1) / 2 + row + 1, col, + m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(col, row - (m_dimension - 1) / 2, + m_nested, false); + } + } + } + + // Odd dimension and lower triangular storage + if (UpLo == Lower && (m_dimension % 2)) { + if (row < col) { + if (row < (m_dimension + 1) / 2) { + return HermitianMatrixCoeffReturnHelper(col, row, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(row - (m_dimension - 1) / 2 - 1, + col - (m_dimension + 1) / 2 + 1, + m_nested, false); + } + } else { + if (col < (m_dimension + 1) / 2) { + return HermitianMatrixCoeffReturnHelper(row, col, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(col - (m_dimension - 1) / 2 - 1, + row - (m_dimension + 1) / 2 + 1, + m_nested, true); + } + } + } + // Some compilers might mock without this fallback return value + return HermitianMatrixCoeffReturnHelper(0, 0, m_nested, true); + } + + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE + HermitianMatrixCoeffReturnHelper coeffRef(Index row, Index col) { + // Even dimension and upper triangular storage + if (UpLo == Upper && !(m_dimension % 2)) { + if (row <= col) { + if (col < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(m_dimension / 2 + col + 1, row, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(row, col - m_dimension / 2, m_nested, true); + } + } else { + if (col <= row) { + if (row < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(m_dimension / 2 + row + 1, col, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(col, row - m_dimension / 2, m_nested, false); + } + } + } + } + + // Even dimension and lower triangular storage + if (UpLo == Lower && !(m_dimension % 2)) { + if (row < col) { + if (row < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(col + 1, row, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(row - m_dimension / 2, col - m_dimension / 2, + m_nested, false); + } + } else { + if (col < m_dimension / 2) { + return HermitianMatrixCoeffReturnHelper(row + 1, col, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(col - m_dimension / 2, row - m_dimension / 2, + m_nested, true); + } + } + } + + // Odd dimension and upper triangular storage + if (UpLo == Upper && (m_dimension % 2)) { + if (row <= col) { + if (col < (m_dimension - 1) / 2) { + return HermitianMatrixCoeffReturnHelper((m_dimension - 1) / 2 + col + 1, row, + m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(row, col - (m_dimension - 1) / 2, + m_nested, true); + } + } else { + if (row < (m_dimension - 1) / 2) { + return HermitianMatrixCoeffReturnHelper((m_dimension - 1) / 2 + row + 1, col, + m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(col, row - (m_dimension - 1) / 2, + m_nested, false); + } + } + } + + // Odd dimension and lower triangular storage + if (UpLo == Lower && (m_dimension % 2)) { + if (row < col) { + if (row < (m_dimension + 1) / 2) { + return HermitianMatrixCoeffReturnHelper(col, row, m_nested, false); + } else { + return HermitianMatrixCoeffReturnHelper(row - (m_dimension - 1) / 2 - 1, + col - (m_dimension + 1) / 2 + 1, + m_nested, false); + } + } else { + if (col < (m_dimension + 1) / 2) { + return HermitianMatrixCoeffReturnHelper(row, col, m_nested, true); + } else { + return HermitianMatrixCoeffReturnHelper(col - (m_dimension - 1) / 2 - 1, + row - (m_dimension + 1) / 2 + 1, + m_nested, true); + } + } + } + // Some compilers might mock without this fallback return value + return HermitianMatrixCoeffReturnHelper(0, 0, m_nested, true); + } + +protected: + const NestedExpression* m_nested; + Index m_dimension; +}; + + +} // namespace internal + +} // end namespace Eigen + +#endif // EIGEN_COREEVALUATORS_H \ No newline at end of file diff --git a/unsupported/test/hermitianmatrix.cpp b/unsupported/test/hermitianmatrix.cpp new file mode 100644 index 000000000..4d69cff05 --- /dev/null +++ b/unsupported/test/hermitianmatrix.cpp @@ -0,0 +1,35 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2012 Kolja Brix +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "../../test/sparse_solver.h" +#include + +template +void test_assignment(int size = HermitianMatrix::ColsAtCompileTime){ + using DenseType = typename HermitianMatrix::DenseType; + + DenseType A=DenseType::Random(size,size); + DenseType B=A+A.tranpose(); + + HermitianMatrix H= B; + DenseType C = H; + VERIFY_IS_APPROX(C, B); + +} + + + +EIGEN_DECLARE_TEST(hermitianMatrix) +{ +CALL_SUBTEST_1(test_assignment>(100)); +CALL_SUBTEST_1(test_assignment>(100)); +CALL_SUBTEST_1(test_assignment>>(100)); + +} -- GitLab From b2dcb07dbae3ffc97d3f2edbd806dbfcbbf1a8bc Mon Sep 17 00:00:00 2001 From: jenswehner Date: Mon, 4 Oct 2021 10:08:29 +0200 Subject: [PATCH 5/5] trying to make hermitian matrix work --- Eigen/src/Core/util/Constants.h | 7 +++++ Eigen/src/Core/util/ForwardDeclarations.h | 8 +++++ Eigen/src/Core/util/XprHelper.h | 29 +++++++++++++++++++ unsupported/Eigen/CMakeLists.txt | 1 + .../Eigen/src/HermitianMatrix/HermitianBase.h | 2 ++ .../src/HermitianMatrix/HermitianMatrix.h | 24 +++++++++++++++ unsupported/test/CMakeLists.txt | 2 +- 7 files changed, 72 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 72fe3307c..eeac25de0 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -520,6 +520,12 @@ struct PermutationStorage {}; /** The type used to identify a permutation storage. */ struct TranspositionsStorage {}; +/** The type used to identify a hermitian storage */ +struct HermitianStorage {}; + +/** The type used to identify a hermitian expression */ +struct HermitianXpr {}; + /** The type used to identify a matrix expression */ struct MatrixXpr {}; @@ -532,6 +538,7 @@ struct SolverShape { static std::string debugName() { return "SolverS struct HomogeneousShape { static std::string debugName() { return "HomogeneousShape"; } }; struct DiagonalShape { static std::string debugName() { return "DiagonalShape"; } }; struct BandShape { static std::string debugName() { return "BandShape"; } }; +struct HermitianShape { static std::string debugName() { return "HermitianShape"; } }; struct TriangularShape { static std::string debugName() { return "TriangularShape"; } }; struct SelfAdjointShape { static std::string debugName() { return "SelfAdjointShape"; } }; struct PermutationShape { static std::string debugName() { return "PermutationShape"; } }; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 6b0ac502b..9e36c3756 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -310,6 +310,14 @@ template class MatrixLogarithmReturnValue; template class MatrixPowerReturnValue; template class MatrixComplexPowerReturnValue; +// Hermitian Matrix module +template class HermitianBase; +template class HermitianMatrix; +template::IsComplex> struct HermitianMatrixCoeffReturnHelper; + namespace internal { template struct stem_function diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 918d48881..49c3ee6bd 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -305,6 +305,11 @@ template struct plain_matrix_type typedef typename T::PlainObject type; }; +template struct plain_matrix_type +{ + typedef typename internal::traits::PlainObject type; +}; + template struct plain_matrix_type_dense { typedef Matrix::Scalar, @@ -351,6 +356,11 @@ template struct eval typedef typename plain_matrix_type::type type; }; +template struct eval +{ + typedef typename plain_matrix_type::type type; +}; + // for matrices, no need to evaluate, just use a const reference to avoid a useless copy template struct eval, Dense> @@ -499,9 +509,26 @@ struct dense_xpr_base typedef ArrayBase type; }; +// Hermitian +template::XprKind> +struct hermitian_xpr_base {}; + +template +struct hermitian_xpr_base +{ + typedef HermitianBase type; +}; + template::XprKind, typename StorageKind = typename traits::StorageKind> struct generic_xpr_base; +template +struct generic_xpr_base +{ + typedef typename hermitian_xpr_base::type type; +}; + + template struct generic_xpr_base { @@ -594,6 +621,8 @@ template struct product_promote_storage_type struct product_promote_storage_type { typedef Dense ret; }; template struct product_promote_storage_type { typedef Dense ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; + /** \internal gives the plain matrix or array type to store a row/column/diagonal of a matrix type. * \tparam Scalar optional parameter allowing to pass a different scalar type than the one of the MatrixType. */ diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt index 631a06014..7882fb64c 100644 --- a/unsupported/Eigen/CMakeLists.txt +++ b/unsupported/Eigen/CMakeLists.txt @@ -6,6 +6,7 @@ set(Eigen_HEADERS BVH EulerAngles FFT + HermitianMatrix IterativeSolvers KroneckerProduct LevenbergMarquardt diff --git a/unsupported/Eigen/src/HermitianMatrix/HermitianBase.h b/unsupported/Eigen/src/HermitianMatrix/HermitianBase.h index b71b0c6df..1448d3c4e 100644 --- a/unsupported/Eigen/src/HermitianMatrix/HermitianBase.h +++ b/unsupported/Eigen/src/HermitianMatrix/HermitianBase.h @@ -11,6 +11,8 @@ #define EIGEN_HERMITIANBASE_H namespace Eigen { + + template class HermitianBase : public EigenBase { diff --git a/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h b/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h index c29d30526..54b8886e9 100644 --- a/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h +++ b/unsupported/Eigen/src/HermitianMatrix/HermitianMatrix.h @@ -11,6 +11,30 @@ #define EIGEN_HERMITIANMATRIX_H namespace Eigen { + + +template +class ProductImpl + : public internal::hermitian_product_base +{ + private: + typedef typename Lhs::NestedExpression NestedExpression; + + public: + typedef typename internal::hermitian_product_base Base; + + EIGEN_DEVICE_FUNC + inline const NestedExpression nestedExpression() const + { + // In the case of a product we just want a brand new nested expression + NestedExpression nested; + return nested; + } + + +}; + + namespace internal { #ifndef EIGEN_PARSED_BY_DOXYGEN diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index fffd0019e..8158c9aa3 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -91,7 +91,7 @@ if(EIGEN_TEST_OPENGL) else() ei_add_property(EIGEN_MISSING_BACKENDS "OpenGL, ") endif() - +ei_add_test(hermitianmatrix) ei_add_test(polynomialsolver) ei_add_test(polynomialutils) ei_add_test(splines) -- GitLab