product_syrk.cpp
Go to the documentation of this file.
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 }
Eigen::Stride
Holds strides information for Map.
Definition: Stride.h:48
s
RealScalar s
Definition: level1_cplx_impl.h:126
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
syrk
void syrk(const MatrixType &m)
Definition: product_syrk.cpp:12
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(product_syrk)
Definition: product_syrk.cpp:131
m1
Matrix3d m1
Definition: IOFormat.cpp:2
buffer
Definition: pytypes.h:2270
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
m2
MatrixType m2(n_dims)
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
m3
static const DiscreteKey m3(M(3), 2)
main.h
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:03:42