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
12template<typename MatrixType> void syrk(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::ColsAtCompileTime, Dynamic> Rhs1;
18  typedef Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> Rhs2;
19  typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic,RowMajor> Rhs3;
20
21  Index rows = m.rows();
22  Index cols = m.cols();
23
24  MatrixType m1 = MatrixType::Random(rows, cols),
25             m2 = MatrixType::Random(rows, cols);
26
27  Rhs1 rhs1 = Rhs1::Random(internal::random<int>(1,320), cols);
28  Rhs2 rhs2 = Rhs2::Random(rows, internal::random<int>(1,320));
29  Rhs3 rhs3 = Rhs3::Random(internal::random<int>(1,320), rows);
30
31  Scalar s1 = internal::random<Scalar>();
32
33  Index c = internal::random<Index>(0,cols-1);
34
35  m2.setZero();
36  VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(rhs2,s1)._expression()),
37                   ((s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
38
39  m2.setZero();
40  VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs2,s1)._expression(),
41                   (s1 * rhs2 * rhs2.adjoint()).eval().template triangularView<Upper>().toDenseMatrix());
42
43  m2.setZero();
44  VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs1.adjoint(),s1)._expression(),
45                   (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Lower>().toDenseMatrix());
46
47  m2.setZero();
48  VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs1.adjoint(),s1)._expression(),
49                   (s1 * rhs1.adjoint() * rhs1).eval().template triangularView<Upper>().toDenseMatrix());
50
51  m2.setZero();
52  VERIFY_IS_APPROX(m2.template selfadjointView<Lower>().rankUpdate(rhs3.adjoint(),s1)._expression(),
53                   (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Lower>().toDenseMatrix());
54
55  m2.setZero();
56  VERIFY_IS_APPROX(m2.template selfadjointView<Upper>().rankUpdate(rhs3.adjoint(),s1)._expression(),
57                   (s1 * rhs3.adjoint() * rhs3).eval().template triangularView<Upper>().toDenseMatrix());
58
59  m2.setZero();
60  VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c),s1)._expression()),
61                   ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
62
63  m2.setZero();
64  VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c),s1)._expression()),
65                   ((s1 * m1.col(c) * m1.col(c).adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
66
67  m2.setZero();
68  VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.col(c).conjugate(),s1)._expression()),
69                   ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
70
71  m2.setZero();
72  VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.col(c).conjugate(),s1)._expression()),
73                   ((s1 * m1.col(c).conjugate() * m1.col(c).conjugate().adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
74
75  m2.setZero();
76  VERIFY_IS_APPROX((m2.template selfadjointView<Lower>().rankUpdate(m1.row(c),s1)._expression()),
77                   ((s1 * m1.row(c).transpose() * m1.row(c).transpose().adjoint()).eval().template triangularView<Lower>().toDenseMatrix()));
78
79  m2.setZero();
80  VERIFY_IS_APPROX((m2.template selfadjointView<Upper>().rankUpdate(m1.row(c).adjoint(),s1)._expression()),
81                   ((s1 * m1.row(c).adjoint() * m1.row(c).adjoint().adjoint()).eval().template triangularView<Upper>().toDenseMatrix()));
82}
83
84void test_product_syrk()
85{
86  for(int i = 0; i < g_repeat ; i++)
87  {
88    int s;
89    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
90    CALL_SUBTEST_1( syrk(MatrixXf(s, s)) );
91    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);
92    CALL_SUBTEST_2( syrk(MatrixXd(s, s)) );
93    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
94    CALL_SUBTEST_3( syrk(MatrixXcf(s, s)) );
95    s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2);
96    CALL_SUBTEST_4( syrk(MatrixXcd(s, s)) );
97  }
98}
99