product_syrk.cpp (7856B)
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 syrk(const MatrixType& m) 13 { 14 typedef typename MatrixType::Scalar Scalar; 15 typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime, RowMajor> RMatrixType; 16 typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> Rhs1; 17 typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2; 18 typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> Rhs3; 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 = MatrixType::Random(rows, cols); 26 RMatrixType rm2 = MatrixType::Random(rows, cols); 27 28 Rhs1 rhs1 = Rhs1::Random(internal::random<int>(1,320), cols); Rhs1 rhs11 = Rhs1::Random(rhs1.rows(), cols); 29 Rhs2 rhs2 = Rhs2::Random(rows, internal::random<int>(1,320)); Rhs2 rhs22 = Rhs2::Random(rows, rhs2.cols()); 30 Rhs3 rhs3 = Rhs3::Random(internal::random<int>(1,320), rows); 31 32 Scalar s1 = internal::random<Scalar>(); 33 34 Index c = internal::random<Index>(0,cols-1); 35 36 m2.setZero(); 37 VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()), 38 ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 39 m2.setZero(); 40 VERIFY_IS_APPROX(((m2.template triangularView<Lower>() += s1 * rhs2 * rhs22.adjoint()).nestedExpression()), 41 ((s1 * rhs2 * rhs22.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 42 43 44 m2.setZero(); 45 VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(), 46 (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix()); 47 m2.setZero(); 48 VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * rhs22 * rhs2.adjoint()).nestedExpression(), 49 (s1 * rhs22 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix()); 50 51 52 m2.setZero(); 53 VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(), 54 (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix()); 55 m2.setZero(); 56 VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * rhs11.adjoint() * rhs1).nestedExpression(), 57 (s1 * rhs11.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix()); 58 59 60 m2.setZero(); 61 VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(), 62 (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix()); 63 VERIFY_IS_APPROX((m2.template triangularView<Upper>() = s1 * rhs1.adjoint() * rhs11).nestedExpression(), 64 (s1 * rhs1.adjoint() * rhs11).eval().template triangularView<Upper>().toDenseMatrix()); 65 66 67 m2.setZero(); 68 VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(), 69 (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix()); 70 71 m2.setZero(); 72 VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs3.adjoint(),s1)._expression(), 73 (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Upper>().toDenseMatrix()); 74 75 m2.setZero(); 76 VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c),s1)._expression()), 77 ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 78 79 m2.setZero(); 80 VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()), 81 ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 82 rm2.setZero(); 83 VERIFY_IS_APPROX((rm2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()), 84 ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 85 m2.setZero(); 86 VERIFY_IS_APPROX((m2.template triangularView<Upper>() += s1 * m3.col(c) * m1.col(c).adjoint()).nestedExpression(), 87 ((s1 * m3.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 88 rm2.setZero(); 89 VERIFY_IS_APPROX((rm2.template triangularView<Upper>() += s1 * m1.col(c) * m3.col(c).adjoint()).nestedExpression(), 90 ((s1 * m1.col(c) * m3.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 91 92 m2.setZero(); 93 VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c).conjugate(),s1)._expression()), 94 ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 95 96 m2.setZero(); 97 VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c).conjugate(),s1)._expression()), 98 ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 99 100 101 m2.setZero(); 102 VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()), 103 ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 104 rm2.setZero(); 105 VERIFY_IS_APPROX((rm2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()), 106 ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 107 m2.setZero(); 108 VERIFY_IS_APPROX((m2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), 109 ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 110 rm2.setZero(); 111 VERIFY_IS_APPROX((rm2.template triangularView<Lower>() += s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).nestedExpression(), 112 ((s1 * m3.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 113 114 115 m2.setZero(); 116 VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.row(c).adjoint(),s1)._expression()), 117 ((s1 * m1.row(c).adjoint() * m1.row(c).adjoint().adjoint()).eval().template triangularView<Upper>().toDenseMatrix())); 118 119 // destination with a non-default inner-stride 120 // see bug 1741 121 { 122 typedef Matrix<Scalar,Dynamic,Dynamic> MatrixX; 123 MatrixX buffer(2*rows,2*cols); 124 Map<MatrixType,0,Stride<Dynamic,2> > map1(buffer.data(),rows,cols,Stride<Dynamic,2>(2*rows,2)); 125 buffer.setZero(); 126 VERIFY_IS_APPROX((map1.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()), 127 ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix())); 128 } 129 } 130 131 EIGEN_DECLARE_TEST(product_syrk) 132 { 133 for(int i = 0; i < g_repeat ; i++) 134 { 135 int s; 136 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE); 137 CALL_SUBTEST_1( syrk(MatrixXf(s, s)) ); 138 CALL_SUBTEST_2( syrk(MatrixXd(s, s)) ); 139 TEST_SET_BUT_UNUSED_VARIABLE(s) 140 141 s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2); 142 CALL_SUBTEST_3( syrk(MatrixXcf(s, s)) ); 143 CALL_SUBTEST_4( syrk(MatrixXcd(s, s)) ); 144 TEST_SET_BUT_UNUSED_VARIABLE(s) 145 } 146 }