aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/Core/SelfAdjointView.h10
-rw-r--r--Eigen/src/Core/products/SelfadjointProduct.h4
-rw-r--r--Eigen/src/Core/products/SelfadjointRank2Update.h4
-rw-r--r--test/CMakeLists.txt6
-rw-r--r--test/product_selfadjoint.cpp68
-rw-r--r--test/product_symm.cpp96
-rw-r--r--test/product_syrk.cpp83
7 files changed, 195 insertions, 76 deletions
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<typename MatrixType, unsigned int UpLo> 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<typename DerivedU, typename DerivedV>
- void rank2update(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha = Scalar(1));
+ SelfAdjointView& rank2update(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& 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<typename DerivedU>
- void rankKupdate(const MatrixBase<DerivedU>& u, Scalar alpha = Scalar(1));
+ SelfAdjointView& rankKupdate(const MatrixBase<DerivedU>& 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<Scalar,MatStorageOrder, ColMajor, AAT, UpLo>
template<typename MatrixType, unsigned int UpLo>
template<typename DerivedU>
-void SelfAdjointView<MatrixType,UpLo>
+SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
::rankKupdate(const MatrixBase<DerivedU>& u, Scalar alpha)
{
typedef ei_blas_traits<DerivedU> UBlasTraits;
@@ -144,6 +144,8 @@ void SelfAdjointView<MatrixType,UpLo>
!UBlasTraits::NeedToConjugate, UpLo>
::run(_expression().cols(), actualU.cols(), &actualU.coeff(0,0), actualU.stride(),
const_cast<Scalar*>(_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<bool Cond, typename T> struct ei_conj_expr_if
template<typename MatrixType, unsigned int UpLo>
template<typename DerivedU, typename DerivedV>
-void SelfAdjointView<MatrixType,UpLo>
+SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo>
::rank2update(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, Scalar alpha)
{
typedef ei_blas_traits<DerivedU> UBlasTraits;
@@ -91,6 +91,8 @@ void SelfAdjointView<MatrixType,UpLo>
typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::ret,
(IsRowMajor ? (UpLo==UpperTriangular ? LowerTriangular : UpperTriangular) : UpLo)>
::run(const_cast<Scalar*>(_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<typename MatrixType> void product_selfadjoint(const MatrixType& m)
}
}
-template<typename MatrixType> void symm(const MatrixType& m)
-{
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1;
- typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2;
- typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> 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<int>(1,320)), rhs12, rhs13;
- Rhs2 rhs2 = Rhs2::Random(ei_random<int>(1,320), rows), rhs22, rhs23;
- Rhs3 rhs3 = Rhs3::Random(cols, ei_random<int>(1,320)), rhs32, rhs33;
-
- Scalar s1 = ei_random<Scalar>(),
- s2 = ei_random<Scalar>();
-
- m2 = m1.template triangularView<LowerTriangular>();
- VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs1),
- rhs13 = (s1*m1) * (s2*rhs1));
-
- m2 = m1.template triangularView<UpperTriangular>();
- VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs1),
- rhs13 = (s1*m1) * (s2*rhs1));
-
- m2 = m1.template triangularView<LowerTriangular>();
- VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()),
- rhs23 = (s1*m1) * (s2*rhs2.adjoint()));
-
- m2 = m1.template triangularView<UpperTriangular>();
- VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs2.adjoint()),
- rhs23 = (s1*m1) * (s2*rhs2.adjoint()));
-
- m2 = m1.template triangularView<UpperTriangular>();
- VERIFY_IS_APPROX(rhs22 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()),
- rhs23 = (s1*m1.adjoint()) * (s2*rhs2.adjoint()));
-
- // test row major = <...>
- m2 = m1.template triangularView<LowerTriangular>();
- VERIFY_IS_APPROX(rhs32 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs3),
- rhs33 = (s1*m1) * (s2 * rhs3));
-
- m2 = m1.template triangularView<UpperTriangular>();
- VERIFY_IS_APPROX(rhs32 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate(),
- rhs33 = (s1*m1.adjoint()) * (s2*rhs3).conjugate());
-
- // test matrix * selfadjoint
- m2 = m1.template triangularView<LowerTriangular>();
- VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<LowerTriangular>(),
- rhs23 = (rhs2) * (m1));
- VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<LowerTriangular>(),
- 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<float,Dynamic,Dynamic,RowMajor>(17,17)) );
CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
}
-
- for(int i = 0; i < g_repeat ; i++)
- {
- int s;
- s = ei_random<int>(10,320);
- CALL_SUBTEST( symm(MatrixXf(s, s)) );
- s = ei_random<int>(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 <gael.guennebaud@gmail.com>
+//
+// 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 <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename MatrixType> void symm(const MatrixType& m)
+{
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1;
+ typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> 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<int>(1,320)), rhs12, rhs13;
+ Rhs2 rhs2 = Rhs2::Random(ei_random<int>(1,320), rows), rhs22, rhs23;
+ Rhs3 rhs3 = Rhs3::Random(cols, ei_random<int>(1,320)), rhs32, rhs33;
+
+ Scalar s1 = ei_random<Scalar>(),
+ s2 = ei_random<Scalar>();
+
+ m2 = m1.template triangularView<LowerTriangular>();
+ VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs1),
+ rhs13 = (s1*m1) * (s2*rhs1));
+
+ m2 = m1.template triangularView<UpperTriangular>();
+ VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs1),
+ rhs13 = (s1*m1) * (s2*rhs1));
+
+ m2 = m1.template triangularView<LowerTriangular>();
+ VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()),
+ rhs23 = (s1*m1) * (s2*rhs2.adjoint()));
+
+ m2 = m1.template triangularView<UpperTriangular>();
+ VERIFY_IS_APPROX(rhs22 = (s1*m2).template selfadjointView<UpperTriangular>() * (s2*rhs2.adjoint()),
+ rhs23 = (s1*m1) * (s2*rhs2.adjoint()));
+
+ m2 = m1.template triangularView<UpperTriangular>();
+ VERIFY_IS_APPROX(rhs22 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs2.adjoint()),
+ rhs23 = (s1*m1.adjoint()) * (s2*rhs2.adjoint()));
+
+ // test row major = <...>
+ m2 = m1.template triangularView<LowerTriangular>();
+ VERIFY_IS_APPROX(rhs32 = (s1*m2).template selfadjointView<LowerTriangular>() * (s2*rhs3),
+ rhs33 = (s1*m1) * (s2 * rhs3));
+
+ m2 = m1.template triangularView<UpperTriangular>();
+ VERIFY_IS_APPROX(rhs32 = (s1*m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate(),
+ rhs33 = (s1*m1.adjoint()) * (s2*rhs3).conjugate());
+
+ // test matrix * selfadjoint
+ m2 = m1.template triangularView<LowerTriangular>();
+ VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<LowerTriangular>(),
+ rhs23 = (rhs2) * (m1));
+ VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<LowerTriangular>(),
+ rhs23 = (s2*rhs2) * (s1*m1));
+}
+void test_product_symm()
+{
+ for(int i = 0; i < g_repeat ; i++)
+ {
+ int s;
+ s = ei_random<int>(10,320);
+ CALL_SUBTEST( symm(MatrixXf(s, s)) );
+ s = ei_random<int>(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 <gael.guennebaud@gmail.com>
+//
+// 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 <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename MatrixType> void syrk(const MatrixType& m)
+{
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1;
+ typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2;
+ typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> 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<int>(1,320), cols);
+ Rhs2 rhs2 = Rhs2::Random(rows, ei_random<int>(1,320));
+ Rhs3 rhs3 = Rhs3::Random(ei_random<int>(1,320), rows);
+
+ Scalar s1 = ei_random<Scalar>(),
+ s2 = ei_random<Scalar>();
+
+ m2.setZero();
+ VERIFY_IS_APPROX((m2.template selfadjointView<LowerTriangular>().rankKupdate(rhs2,s1)._expression()),
+ ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<LowerTriangular>().toDense()));
+
+ m2.setZero();
+ VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankKupdate(rhs2,s1)._expression(),
+ (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<UpperTriangular>().toDense());
+
+ m2.setZero();
+ VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankKupdate(rhs1.adjoint(),s1)._expression(),
+ (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<LowerTriangular>().toDense());
+
+ m2.setZero();
+ VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankKupdate(rhs1.adjoint(),s1)._expression(),
+ (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<UpperTriangular>().toDense());
+
+ m2.setZero();
+ VERIFY_IS_APPROX(m2.template selfadjointView<LowerTriangular>().rankKupdate(rhs3.adjoint(),s1)._expression(),
+ (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<LowerTriangular>().toDense());
+
+ m2.setZero();
+ VERIFY_IS_APPROX(m2.template selfadjointView<UpperTriangular>().rankKupdate(rhs3.adjoint(),s1)._expression(),
+ (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<UpperTriangular>().toDense());
+}
+
+void test_product_syrk()
+{
+ for(int i = 0; i < g_repeat ; i++)
+ {
+ int s;
+ s = ei_random<int>(10,320);
+ CALL_SUBTEST( syrk(MatrixXf(s, s)) );
+ s = ei_random<int>(10,320);
+ CALL_SUBTEST( syrk(MatrixXcd(s, s)) );
+ }
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 10:44:31 +0000