diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index b1b4f9e32..4787bcbd8 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -120,23 +120,25 @@ template class SelfAdjointView /** Perform a symmetric rank 2 update of the selfadjoint matrix \c *this: * \f$ this = this + \alpha ( u v^* + v u^*) \f$ + * \returns a reference to \c *this * * The vectors \a u and \c v \b must be column vectors, however they can be * a adjoint expression without any overhead. Only the meaningful triangular * part of the matrix is updated, the rest is left unchanged. */ template - void rank2update(const MatrixBase& u, const MatrixBase& v, Scalar alpha = Scalar(1)); + SelfAdjointView& rank2update(const MatrixBase& u, const MatrixBase& v, Scalar alpha = Scalar(1)); /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: - * \f$ this = this + \alpha ( u u^* ) \f$ - * where \a u is a vector or matrix. + * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. + * + * \returns a reference to \c *this * * Note that to perform \f$ this = this + \alpha ( u^* u ) \f$ you can simply * call this function with u.adjoint(). */ template - void rankKupdate(const MatrixBase& u, Scalar alpha = Scalar(1)); + SelfAdjointView& rankKupdate(const MatrixBase& u, Scalar alpha = Scalar(1)); /////////// Cholesky module /////////// diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h index b4ca4e786..08df9e15c 100644 --- a/Eigen/src/Core/products/SelfadjointProduct.h +++ b/Eigen/src/Core/products/SelfadjointProduct.h @@ -126,7 +126,7 @@ struct ei_selfadjoint_product template template -void SelfAdjointView +SelfAdjointView& SelfAdjointView ::rankKupdate(const MatrixBase& u, Scalar alpha) { typedef ei_blas_traits UBlasTraits; @@ -144,6 +144,8 @@ void SelfAdjointView !UBlasTraits::NeedToConjugate, UpLo> ::run(_expression().cols(), actualU.cols(), &actualU.coeff(0,0), actualU.stride(), const_cast(_expression().data()), _expression().stride(), actualAlpha); + + return *this; } diff --git a/Eigen/src/Core/products/SelfadjointRank2Update.h b/Eigen/src/Core/products/SelfadjointRank2Update.h index edb57ecd5..6c8a28f65 100644 --- a/Eigen/src/Core/products/SelfadjointRank2Update.h +++ b/Eigen/src/Core/products/SelfadjointRank2Update.h @@ -69,7 +69,7 @@ template struct ei_conj_expr_if template template -void SelfAdjointView +SelfAdjointView& SelfAdjointView ::rank2update(const MatrixBase& u, const MatrixBase& v, Scalar alpha) { typedef ei_blas_traits UBlasTraits; @@ -91,6 +91,8 @@ void SelfAdjointView typename ei_conj_expr_if::ret, (IsRowMajor ? (UpLo==UpperTriangular ? LowerTriangular : UpperTriangular) : UpLo)> ::run(const_cast(_expression().data()),_expression().stride(),actualU,actualV,actualAlpha); + + return *this; } #endif // EIGEN_SELFADJOINTRANK2UPTADE_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5e1adcbf4..d1c4d49e2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -97,8 +97,10 @@ ei_add_test(cwiseop) ei_add_test(redux) ei_add_test(product_small) ei_add_test(product_large ${EI_OFLAG}) -ei_add_test(product_selfadjoint) -ei_add_test(product_extra) +ei_add_test(product_extra ${EI_OFLAG}) +ei_add_test(product_selfadjoint ${EI_OFLAG}) +ei_add_test(product_symm ${EI_OFLAG}) +ei_add_test(product_syrk ${EI_OFLAG}) ei_add_test(diagonalmatrices) ei_add_test(adjoint) ei_add_test(submatrices) diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp index 2cacc8e5e..efa487ab1 100644 --- a/test/product_selfadjoint.cpp +++ b/test/product_selfadjoint.cpp @@ -94,65 +94,6 @@ template void product_selfadjoint(const MatrixType& m) } } -template void symm(const MatrixType& m) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - typedef Matrix Rhs1; - typedef Matrix Rhs2; - typedef Matrix Rhs3; - - int rows = m.rows(); - int cols = m.cols(); - - MatrixType m1 = MatrixType::Random(rows, cols), - m2 = MatrixType::Random(rows, cols); - - m1 = (m1+m1.adjoint()).eval(); - - Rhs1 rhs1 = Rhs1::Random(cols, ei_random(1,320)), rhs12, rhs13; - Rhs2 rhs2 = Rhs2::Random(ei_random(1,320), rows), rhs22, rhs23; - Rhs3 rhs3 = Rhs3::Random(cols, ei_random(1,320)), rhs32, rhs33; - - Scalar s1 = ei_random(), - s2 = ei_random(); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView() * (s2*rhs1), - rhs13 = (s1*m1) * (s2*rhs1)); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView() * (s2*rhs1), - rhs13 = (s1*m1) * (s2*rhs1)); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView() * (s2*rhs2.adjoint()), - rhs23 = (s1*m1) * (s2*rhs2.adjoint())); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView() * (s2*rhs2.adjoint()), - rhs23 = (s1*m1) * (s2*rhs2.adjoint())); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs22 = (s1*m2.adjoint()).template selfadjointView() * (s2*rhs2.adjoint()), - rhs23 = (s1*m1.adjoint()) * (s2*rhs2.adjoint())); - - // test row major = <...> - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs32 = (s1*m2).template selfadjointView() * (s2*rhs3), - rhs33 = (s1*m1) * (s2 * rhs3)); - - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs32 = (s1*m2.adjoint()).template selfadjointView() * (s2*rhs3).conjugate(), - rhs33 = (s1*m1.adjoint()) * (s2*rhs3).conjugate()); - - // test matrix * selfadjoint - m2 = m1.template triangularView(); - VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView(), - rhs23 = (rhs2) * (m1)); - VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView(), - rhs23 = (s2*rhs2) * (s1*m1)); -} void test_product_selfadjoint() { for(int i = 0; i < g_repeat ; i++) { @@ -165,13 +106,4 @@ void test_product_selfadjoint() CALL_SUBTEST( product_selfadjoint(Matrix(17,17)) ); CALL_SUBTEST( product_selfadjoint(Matrix,Dynamic,Dynamic,RowMajor>(19, 19)) ); } - - for(int i = 0; i < g_repeat ; i++) - { - int s; - s = ei_random(10,320); - CALL_SUBTEST( symm(MatrixXf(s, s)) ); - s = ei_random(10,320); - CALL_SUBTEST( symm(MatrixXcd(s, s)) ); - } } diff --git a/test/product_symm.cpp b/test/product_symm.cpp new file mode 100644 index 000000000..3a0cd94d0 --- /dev/null +++ b/test/product_symm.cpp @@ -0,0 +1,96 @@ +// This file is part 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 symm(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix Rhs1; + typedef Matrix Rhs2; + typedef Matrix Rhs3; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols); + + m1 = (m1+m1.adjoint()).eval(); + + Rhs1 rhs1 = Rhs1::Random(cols, ei_random(1,320)), rhs12, rhs13; + Rhs2 rhs2 = Rhs2::Random(ei_random(1,320), rows), rhs22, rhs23; + Rhs3 rhs3 = Rhs3::Random(cols, ei_random(1,320)), rhs32, rhs33; + + Scalar s1 = ei_random(), + s2 = ei_random(); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView() * (s2*rhs1), + rhs13 = (s1*m1) * (s2*rhs1)); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView() * (s2*rhs1), + rhs13 = (s1*m1) * (s2*rhs1)); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView() * (s2*rhs2.adjoint()), + rhs23 = (s1*m1) * (s2*rhs2.adjoint())); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView() * (s2*rhs2.adjoint()), + rhs23 = (s1*m1) * (s2*rhs2.adjoint())); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = (s1*m2.adjoint()).template selfadjointView() * (s2*rhs2.adjoint()), + rhs23 = (s1*m1.adjoint()) * (s2*rhs2.adjoint())); + + // test row major = <...> + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs32 = (s1*m2).template selfadjointView() * (s2*rhs3), + rhs33 = (s1*m1) * (s2 * rhs3)); + + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs32 = (s1*m2.adjoint()).template selfadjointView() * (s2*rhs3).conjugate(), + rhs33 = (s1*m1.adjoint()) * (s2*rhs3).conjugate()); + + // test matrix * selfadjoint + m2 = m1.template triangularView(); + VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView(), + rhs23 = (rhs2) * (m1)); + VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView(), + rhs23 = (s2*rhs2) * (s1*m1)); +} +void test_product_symm() +{ + for(int i = 0; i < g_repeat ; i++) + { + int s; + s = ei_random(10,320); + CALL_SUBTEST( symm(MatrixXf(s, s)) ); + s = ei_random(10,320); + CALL_SUBTEST( symm(MatrixXcd(s, s)) ); + } +} diff --git a/test/product_syrk.cpp b/test/product_syrk.cpp new file mode 100644 index 000000000..ff8bf933d --- /dev/null +++ b/test/product_syrk.cpp @@ -0,0 +1,83 @@ +// This file is part 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 syrk(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix Rhs1; + typedef Matrix Rhs2; + typedef Matrix Rhs3; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols); + + Rhs1 rhs1 = Rhs1::Random(ei_random(1,320), cols); + Rhs2 rhs2 = Rhs2::Random(rows, ei_random(1,320)); + Rhs3 rhs3 = Rhs3::Random(ei_random(1,320), rows); + + Scalar s1 = ei_random(), + s2 = ei_random(); + + m2.setZero(); + VERIFY_IS_APPROX((m2.template selfadjointView().rankKupdate(rhs2,s1)._expression()), + ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView().toDense())); + + m2.setZero(); + VERIFY_IS_APPROX(m2.template selfadjointView().rankKupdate(rhs2,s1)._expression(), + (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView().toDense()); + + m2.setZero(); + VERIFY_IS_APPROX(m2.template selfadjointView().rankKupdate(rhs1.adjoint(),s1)._expression(), + (s1 * rhs1.adjoint() * rhs1).eval().template triangularView().toDense()); + + m2.setZero(); + VERIFY_IS_APPROX(m2.template selfadjointView().rankKupdate(rhs1.adjoint(),s1)._expression(), + (s1 * rhs1.adjoint() * rhs1).eval().template triangularView().toDense()); + + m2.setZero(); + VERIFY_IS_APPROX(m2.template selfadjointView().rankKupdate(rhs3.adjoint(),s1)._expression(), + (s1 * rhs3.adjoint() * rhs3).eval().template triangularView().toDense()); + + m2.setZero(); + VERIFY_IS_APPROX(m2.template selfadjointView().rankKupdate(rhs3.adjoint(),s1)._expression(), + (s1 * rhs3.adjoint() * rhs3).eval().template triangularView().toDense()); +} + +void test_product_syrk() +{ + for(int i = 0; i < g_repeat ; i++) + { + int s; + s = ei_random(10,320); + CALL_SUBTEST( syrk(MatrixXf(s, s)) ); + s = ei_random(10,320); + CALL_SUBTEST( syrk(MatrixXcd(s, s)) ); + } +}