Add an efficient rank2 update function (like the level2 blas xSYR2 routine).

Note that it is already used in Tridiagonalization.

2 files changed, 36 insertions, 8 deletions

diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index c93953714..6b5092775 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -119,8 +119,8 @@ void test_eigensolver_selfadjoint()
     // very important to test a 3x3 matrix since we provide a special path for it
     CALL_SUBTEST( selfadjointeigensolver(Matrix3f()) );
     CALL_SUBTEST( selfadjointeigensolver(Matrix4d()) );
-    CALL_SUBTEST( selfadjointeigensolver(MatrixXf(7,7)) );
-    CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(5,5)) );
+    CALL_SUBTEST( selfadjointeigensolver(MatrixXf(4,4)) );
+    CALL_SUBTEST( selfadjointeigensolver(MatrixXcd(7,7)) );
     CALL_SUBTEST( selfadjointeigensolver(MatrixXd(19,19)) );
 
     // some trivial but implementation-wise tricky cases
diff --git a/test/product_selfadjoint.cpp b/test/product_selfadjoint.cpp
index b26b7223b..297bab1a9 100644
--- a/test/product_selfadjoint.cpp
+++ b/test/product_selfadjoint.cpp
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@gmail.com>
+// 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
@@ -29,20 +29,29 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
   typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+  typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType;
 
   int rows = m.rows();
   int cols = m.cols();
 
   MatrixType m1 = MatrixType::Random(rows, cols),
-             m2 = MatrixType::Random(rows, cols);
+             m2 = MatrixType::Random(rows, cols),
+             m3;
   VectorType v1 = VectorType::Random(rows),
              v2 = VectorType::Random(rows);
+
+  RowVectorType r1 = RowVectorType::Random(rows),
+                r2 = RowVectorType::Random(rows);
+
+  Scalar s1 = ei_random<Scalar>(),
+         s2 = ei_random<Scalar>(),
+         s3 = ei_random<Scalar>();
   
   m1 = m1.adjoint()*m1;
   
   // lower
   m2.setZero();
-  m2.template part<LowerTriangular>() = m1;
+  m2.template triangularView<LowerTriangular>() = m1;
   ei_product_selfadjoint_vector<Scalar,MatrixType::Flags&RowMajorBit,LowerTriangularBit>
     (cols,m2.data(),cols, v1.data(), v2.data());
   VERIFY_IS_APPROX(v2, m1 * v1);
@@ -50,11 +59,30 @@ template<typename MatrixType> void product_selfadjoint(const MatrixType& m)
 
   // upper
   m2.setZero();
-  m2.template part<UpperTriangular>() = m1;
+  m2.template triangularView<UpperTriangular>() = m1;
   ei_product_selfadjoint_vector<Scalar,MatrixType::Flags&RowMajorBit,UpperTriangularBit>(cols,m2.data(),cols, v1.data(), v2.data());
   VERIFY_IS_APPROX(v2, m1 * v1);
   VERIFY_IS_APPROX((m2.template selfadjointView<UpperTriangular>() * v1).eval(), m1 * v1);
 
+  // rank2 update
+  m2 = m1.template triangularView<LowerTriangular>();
+  m2.template selfadjointView<LowerTriangular>().rank2update(v1,v2);
+  VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<LowerTriangular>().toDense());
+  
+  m2 = m1.template triangularView<UpperTriangular>();
+  m2.template selfadjointView<UpperTriangular>().rank2update(-v1,s2*v2,s3);
+  VERIFY_IS_APPROX(m2, (m1 + (-s2*s3) * (v1 * v2.adjoint()+ v2 * v1.adjoint())).template triangularView<UpperTriangular>().toDense());
+
+  m2 = m1.template triangularView<UpperTriangular>();
+  m2.template selfadjointView<UpperTriangular>().rank2update(-r1.adjoint(),r2.adjoint()*s3,s1);
+  VERIFY_IS_APPROX(m2, (m1 + (-s3*s1) * (r1.adjoint() * r2 + r2.adjoint() * r1)).template triangularView<UpperTriangular>().toDense());
+
+  m2 = m1.template triangularView<LowerTriangular>();
+  m2.block(1,1,rows-1,cols-1).template selfadjointView<LowerTriangular>().rank2update(v1.end(rows-1),v2.start(cols-1));
+  m3 = m1;
+  m3.block(1,1,rows-1,cols-1) += v1.end(rows-1) * v2.start(cols-1).adjoint()+ v2.start(cols-1) * v1.end(rows-1).adjoint();
+  VERIFY_IS_APPROX(m2, m3.template triangularView<LowerTriangular>().toDense());
+
 }
 
 void test_product_selfadjoint()
@@ -65,8 +93,8 @@ void test_product_selfadjoint()
     CALL_SUBTEST( product_selfadjoint(Matrix3d()) );
     CALL_SUBTEST( product_selfadjoint(MatrixXcf(4, 4)) );
     CALL_SUBTEST( product_selfadjoint(MatrixXcd(21,21)) );
-    CALL_SUBTEST( product_selfadjoint(MatrixXd(17,17)) );
-    CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(18,18)) );
+    CALL_SUBTEST( product_selfadjoint(MatrixXd(4,4)) );
+    CALL_SUBTEST( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(17,17)) );
     CALL_SUBTEST( product_selfadjoint(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(19, 19)) );
   }
 }