product_selfadjoint.cpp (3510B)
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@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::Scalar Scalar; 15 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; 16 typedef Matrix<Scalar, 1, MatrixType::RowsAtCompileTime> RowVectorType; 17 18 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic, RowMajor> RhsMatrixType; 19 20 Index rows = m.rows(); 21 Index cols = m.cols(); 22 23 MatrixType m1 = MatrixType::Random(rows, cols), 24 m2 = MatrixType::Random(rows, cols), 25 m3; 26 VectorType v1 = VectorType::Random(rows), 27 v2 = VectorType::Random(rows), 28 v3(rows); 29 RowVectorType r1 = RowVectorType::Random(rows), 30 r2 = RowVectorType::Random(rows); 31 RhsMatrixType m4 = RhsMatrixType::Random(rows,10); 32 33 Scalar s1 = internal::random<Scalar>(), 34 s2 = internal::random<Scalar>(), 35 s3 = internal::random<Scalar>(); 36 37 m1 = (m1.adjoint() + m1).eval(); 38 39 // rank2 update 40 m2 = m1.template triangularView<Lower>(); 41 m2.template selfadjointView<Lower>().rankUpdate(v1,v2); 42 VERIFY_IS_APPROX(m2, (m1 + v1 * v2.adjoint()+ v2 * v1.adjoint()).template triangularView<Lower>().toDenseMatrix()); 43 44 m2 = m1.template triangularView<Upper>(); 45 m2.template selfadjointView<Upper>().rankUpdate(-v1,s2*v2,s3); 46 VERIFY_IS_APPROX(m2, (m1 + (s3*(-v1)*(s2*v2).adjoint()+numext::conj(s3)*(s2*v2)*(-v1).adjoint())).template triangularView<Upper>().toDenseMatrix()); 47 48 m2 = m1.template triangularView<Upper>(); 49 m2.template selfadjointView<Upper>().rankUpdate(-s2*r1.adjoint(),r2.adjoint()*s3,s1); 50 VERIFY_IS_APPROX(m2, (m1 + s1*(-s2*r1.adjoint())*(r2.adjoint()*s3).adjoint() + numext::conj(s1)*(r2.adjoint()*s3) * (-s2*r1.adjoint()).adjoint()).template triangularView<Upper>().toDenseMatrix()); 51 52 if (rows>1) 53 { 54 m2 = m1.template triangularView<Lower>(); 55 m2.block(1,1,rows-1,cols-1).template selfadjointView<Lower>().rankUpdate(v1.tail(rows-1),v2.head(cols-1)); 56 m3 = m1; 57 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(); 58 VERIFY_IS_APPROX(m2, m3.template triangularView<Lower>().toDenseMatrix()); 59 } 60 } 61 62 EIGEN_DECLARE_TEST(product_selfadjoint) 63 { 64 int s = 0; 65 for(int i = 0; i < g_repeat ; i++) { 66 CALL_SUBTEST_1( product_selfadjoint(Matrix<float, 1, 1>()) ); 67 CALL_SUBTEST_2( product_selfadjoint(Matrix<float, 2, 2>()) ); 68 CALL_SUBTEST_3( product_selfadjoint(Matrix3d()) ); 69 70 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 71 CALL_SUBTEST_4( product_selfadjoint(MatrixXcf(s, s)) ); 72 TEST_SET_BUT_UNUSED_VARIABLE(s) 73 74 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 75 CALL_SUBTEST_5( product_selfadjoint(MatrixXcd(s,s)) ); 76 TEST_SET_BUT_UNUSED_VARIABLE(s) 77 78 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 79 CALL_SUBTEST_6( product_selfadjoint(MatrixXd(s,s)) ); 80 TEST_SET_BUT_UNUSED_VARIABLE(s) 81 82 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 83 CALL_SUBTEST_7( product_selfadjoint(Matrix<float,Dynamic,Dynamic,RowMajor>(s,s)) ); 84 TEST_SET_BUT_UNUSED_VARIABLE(s) 85 } 86 }