1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "main.h" 11 12 template<typename MatrixType> void product_selfadjoint(const MatrixType& m) 13 { 14 typedef typename MatrixType::Index Index; 15 typedef typename MatrixType::Scalar Scalar; 16 typedef typename NumTraits<Scalar>::Real RealScalar; 17 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 18 typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType; 19 20 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic, RowMajor> RhsMatrixType; 21 22 Index rows = m.rows(); 23 Index cols = m.cols(); 24 25 MatrixType m1 = MatrixType::Random(rows, cols), 26 m2 = MatrixType::Random(rows, cols), 27 m3; 28 VectorType v1 = VectorType::Random(rows), 29 v2 = VectorType::Random(rows), 30 v3(rows); 31 RowVectorType r1 = RowVectorType::Random(rows), 32 r2 = RowVectorType::Random(rows); 33 RhsMatrixType m4 = RhsMatrixType::Random(rows,10); 34 35 Scalar s1 = internal::random<Scalar>(), 36 s2 = internal::random<Scalar>(), 37 s3 = internal::random<Scalar>(); 38 39 m1 = (m1.adjoint() + m1).eval(); 40 41 // rank2 update 42 m2 = m1.template triangularView<Lower>(); 43 m2.template selfadjointView<Lower>().rankUpdate(v1,v2); 44 VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<Lower>().toDenseMatrix()); 45 46 m2 = m1.template triangularView<Upper>(); 47 m2.template selfadjointView<Upper>().rankUpdate(-v1,s2*v2,s3); 48 VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+internal::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView<Upper>().toDenseMatrix()); 49 50 m2 = m1.template triangularView<Upper>(); 51 m2.template selfadjointView<Upper>().rankUpdate(-s2*r1.adjoint(),r2.adjoint()*s3,s1); 52 VERIFY_IS_APPROX(m2, (m1 + s1*(-s2*r1.adjoint())*(r2.adjoint()*s3).adjoint() + internal::conj(s1)*(r2.adjoint()*s3) * (-s2*r1.adjoint()).adjoint()).template triangularView<Upper>().toDenseMatrix()); 53 54 if (rows>1) 55 { 56 m2 = m1.template triangularView<Lower>(); 57 m2.block(1,1,rows-1,cols-1).template selfadjointView<Lower>().rankUpdate(v1.tail(rows-1),v2.head(cols-1)); 58 m3 = m1; 59 m3.block(1,1,rows-1,cols-1) += v1.tail(rows-1) * v2.head(cols-1).adjoint()+ v2.head(cols-1) * v1.tail(rows-1).adjoint(); 60 VERIFY_IS_APPROX(m2, m3.template triangularView<Lower>().toDenseMatrix()); 61 } 62 } 63 64 void test_product_selfadjoint() 65 { 66 int s; 67 for(int i = 0; i < g_repeat ; i++) { 68 CALL_SUBTEST_1( product_selfadjoint(Matrix<float, 1, 1>()) ); 69 CALL_SUBTEST_2( product_selfadjoint(Matrix<float, 2, 2>()) ); 70 CALL_SUBTEST_3( product_selfadjoint(Matrix3d()) ); 71 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 72 CALL_SUBTEST_4( product_selfadjoint(MatrixXcf(s, s)) ); 73 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 74 CALL_SUBTEST_5( product_selfadjoint(MatrixXcd(s,s)) ); 75 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 76 CALL_SUBTEST_6( product_selfadjoint(MatrixXd(s,s)) ); 77 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 78 CALL_SUBTEST_7( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(s,s)) ); 79 } 80 EIGEN_UNUSED_VARIABLE(s) 81 } 82
