From bf09fe55edd7114dfb150a8a9dbb62965e3e3201 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 19 May 2010 16:35:34 +0200 Subject: [PATCH] fix selfadjoint to dense --- Eigen/src/Core/SelfAdjointView.h | 72 ++++++++++++++++++++++++++----- Eigen/src/Core/TriangularMatrix.h | 34 --------------- test/CMakeLists.txt | 1 + test/selfadjoint.cpp | 66 ++++++++++++++++++++++++++++ 4 files changed, 129 insertions(+), 44 deletions(-) create mode 100644 test/selfadjoint.cpp diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 2975bb2bc..f6ae05511 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -39,14 +39,14 @@ * * \sa class TriangularBase, MatrixBase::selfAdjointView() */ -template -struct ei_traits > : ei_traits +template +struct ei_traits > : ei_traits { typedef typename ei_nested::type MatrixTypeNested; typedef typename ei_unref::type _MatrixTypeNested; typedef MatrixType ExpressionType; enum { - Mode = TriangularPart | SelfAdjoint, + Mode = UpLo | SelfAdjoint, Flags = _MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit)), // FIXME these flags should be preserved CoeffReadCost = _MatrixTypeNested::CoeffReadCost @@ -77,7 +77,7 @@ template class SelfAdjointView inline int cols() const { return m_matrix.cols(); } inline int outerStride() const { return m_matrix.outerStride(); } inline int innerStride() const { return m_matrix.innerStride(); } - + /** \sa MatrixBase::coeff() * \warning the coordinates must fit into the referenced triangular part */ @@ -165,8 +165,10 @@ template class SelfAdjointView // return ei_matrix_selfadjoint_product_returntype >(lhs.derived(),rhs); // } +// selfadjoint to dense matrix + template -struct ei_triangular_assignment_selector +struct ei_triangular_assignment_selector { enum { col = (UnrollCount-1) / Derived1::RowsAtCompileTime, @@ -175,7 +177,7 @@ struct ei_triangular_assignment_selector::run(dst, src); + ei_triangular_assignment_selector::run(dst, src); if(row == col) dst.coeffRef(row, col) = ei_real(src.coeff(row, col)); @@ -184,17 +186,67 @@ struct ei_triangular_assignment_selector -struct ei_triangular_assignment_selector +struct ei_triangular_assignment_selector +{ + inline static void run(Derived1 &, const Derived2 &) {} +}; + +template +struct ei_triangular_assignment_selector +{ + enum { + col = (UnrollCount-1) / Derived1::RowsAtCompileTime, + row = (UnrollCount-1) % Derived1::RowsAtCompileTime + }; + + inline static void run(Derived1 &dst, const Derived2 &src) + { + ei_triangular_assignment_selector::run(dst, src); + + if(row == col) + dst.coeffRef(row, col) = ei_real(src.coeff(row, col)); + else if(row > col) + dst.coeffRef(col, row) = ei_conj(dst.coeffRef(row, col) = src.coeff(row, col)); + } +}; + +template +struct ei_triangular_assignment_selector +{ + inline static void run(Derived1 &, const Derived2 &) {} +}; + +template +struct ei_triangular_assignment_selector { inline static void run(Derived1 &dst, const Derived2 &src) { for(int j = 0; j < dst.cols(); ++j) { for(int i = 0; i < j; ++i) - dst.coeffRef(j, i) = ei_conj(dst.coeffRef(i, j) = src.coeff(i, j)); - dst.coeffRef(j, j) = ei_real(src.coeff(j, j)); + { + dst.copyCoeff(i, j, src); + dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j)); + } + dst.copyCoeff(j, j, src); + } + } +}; + +template +struct ei_triangular_assignment_selector +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + for(int i = 0; i < dst.rows(); ++i) + { + for(int j = 0; j < i; ++j) + { + dst.copyCoeff(i, j, src); + dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j)); + } + dst.copyCoeff(i, i, src); } } }; diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 383a3edd3..2230680d1 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -472,40 +472,6 @@ struct ei_triangular_assignment_selector -struct ei_triangular_assignment_selector -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int j = 0; j < dst.cols(); ++j) - { - for(int i = 0; i < j; ++i) - { - dst.copyCoeff(i, j, src); - dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j)); - } - dst.copyCoeff(j, j, src); - } - } -}; - -template -struct ei_triangular_assignment_selector -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - for(int i = 0; i < dst.rows(); ++i) - { - for(int j = 0; j < i; ++j) - { - dst.copyCoeff(i, j, src); - dst.coeffRef(j,i) = ei_conj(dst.coeff(i,j)); - } - dst.copyCoeff(i, i, src); - } - } -}; - // FIXME should we keep that possibility template template diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 82045ff90..4b5ee0d92 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -116,6 +116,7 @@ ei_add_test(array_for_matrix) ei_add_test(array_replicate) ei_add_test(array_reverse) ei_add_test(triangular) +ei_add_test(selfadjoint) ei_add_test(product_selfadjoint) ei_add_test(product_symm) ei_add_test(product_syrk) diff --git a/test/selfadjoint.cpp b/test/selfadjoint.cpp new file mode 100644 index 000000000..2eaf6d459 --- /dev/null +++ b/test/selfadjoint.cpp @@ -0,0 +1,66 @@ +// This file is triangularView of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2010 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" + +// This file tests the basic selfadjointView API, +// the related products and decompositions are tested in specific files. + +template void selfadjoint(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m3(rows, cols); + + m1.diagonal() = m1.diagonal().real().template cast(); + + // check selfadjoint to dense + m3 = m1.template selfadjointView(); + VERIFY_IS_APPROX(MatrixType(m3.template triangularView()), MatrixType(m1.template triangularView())); + VERIFY_IS_APPROX(m3, m3.adjoint()); + + + m3 = m1.template selfadjointView(); + VERIFY_IS_APPROX(MatrixType(m3.template triangularView()), MatrixType(m1.template triangularView())); + VERIFY_IS_APPROX(m3, m3.adjoint()); +} + +void test_selfadjoint() +{ + for(int i = 0; i < g_repeat ; i++) + { + EIGEN_UNUSED int s = ei_random(1,20); + + CALL_SUBTEST_1( selfadjoint(Matrix()) ); + CALL_SUBTEST_2( selfadjoint(Matrix()) ); + CALL_SUBTEST_3( selfadjoint(Matrix3cf()) ); + CALL_SUBTEST_4( selfadjoint(MatrixXcd(s,s)) ); + CALL_SUBTEST_5( selfadjoint(Matrix(s, s)) ); + } +}