product_symm.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 Scalar, int Size, int OtherSize> void symm(int size = Size, int othersize = OtherSize)
13 {
14  typedef Matrix<Scalar, Size, Size> MatrixType;
15  typedef Matrix<Scalar, Size, OtherSize> Rhs1;
16  typedef Matrix<Scalar, OtherSize, Size> Rhs2;
17  enum { order = OtherSize==1 ? 0 : RowMajor };
18  typedef Matrix<Scalar, Size, OtherSize,order> Rhs3;
19 
20  Index rows = size;
21  Index cols = size;
22 
23  MatrixType m1 = MatrixType::Random(rows, cols),
24  m2 = MatrixType::Random(rows, cols), m3;
25 
26  m1 = (m1+m1.adjoint()).eval();
27 
28  Rhs1 rhs1 = Rhs1::Random(cols, othersize), rhs12(cols, othersize), rhs13(cols, othersize);
29  Rhs2 rhs2 = Rhs2::Random(othersize, rows), rhs22(othersize, rows), rhs23(othersize, rows);
30  Rhs3 rhs3 = Rhs3::Random(cols, othersize), rhs32(cols, othersize), rhs33(cols, othersize);
31 
32  Scalar s1 = internal::random<Scalar>(),
33  s2 = internal::random<Scalar>();
34 
35  m2 = m1.template triangularView<Lower>();
36  m3 = m2.template selfadjointView<Lower>();
37  VERIFY_IS_EQUAL(m1, m3);
38  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs1),
39  rhs13 = (s1*m1) * (s2*rhs1));
40 
41  VERIFY_IS_APPROX(rhs12 = (s1*m2).transpose().template selfadjointView<Upper>() * (s2*rhs1),
42  rhs13 = (s1*m1.transpose()) * (s2*rhs1));
43 
44  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>().transpose() * (s2*rhs1),
45  rhs13 = (s1*m1.transpose()) * (s2*rhs1));
46 
47  VERIFY_IS_APPROX(rhs12 = (s1*m2).conjugate().template selfadjointView<Lower>() * (s2*rhs1),
48  rhs13 = (s1*m1).conjugate() * (s2*rhs1));
49 
50  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>().conjugate() * (s2*rhs1),
51  rhs13 = (s1*m1).conjugate() * (s2*rhs1));
52 
53  VERIFY_IS_APPROX(rhs12 = (s1*m2).adjoint().template selfadjointView<Upper>() * (s2*rhs1),
54  rhs13 = (s1*m1).adjoint() * (s2*rhs1));
55 
56  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>().adjoint() * (s2*rhs1),
57  rhs13 = (s1*m1).adjoint() * (s2*rhs1));
58 
59  m2 = m1.template triangularView<Upper>(); rhs12.setRandom(); rhs13 = rhs12;
60  m3 = m2.template selfadjointView<Upper>();
61  VERIFY_IS_EQUAL(m1, m3);
62  VERIFY_IS_APPROX(rhs12 += (s1*m2).template selfadjointView<Upper>() * (s2*rhs1),
63  rhs13 += (s1*m1) * (s2*rhs1));
64 
65  m2 = m1.template triangularView<Lower>();
66  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Lower>() * (s2*rhs2.adjoint()),
67  rhs13 = (s1*m1) * (s2*rhs2.adjoint()));
68 
69  m2 = m1.template triangularView<Upper>();
70  VERIFY_IS_APPROX(rhs12 = (s1*m2).template selfadjointView<Upper>() * (s2*rhs2.adjoint()),
71  rhs13 = (s1*m1) * (s2*rhs2.adjoint()));
72 
73  m2 = m1.template triangularView<Upper>();
74  VERIFY_IS_APPROX(rhs12 = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs2.adjoint()),
75  rhs13 = (s1*m1.adjoint()) * (s2*rhs2.adjoint()));
76 
77  // test row major = <...>
78  m2 = m1.template triangularView<Lower>(); rhs32.setRandom(); rhs13 = rhs32;
79  VERIFY_IS_APPROX(rhs32.noalias() -= (s1*m2).template selfadjointView<Lower>() * (s2*rhs3),
80  rhs13 -= (s1*m1) * (s2 * rhs3));
81 
82  m2 = m1.template triangularView<Upper>();
83  VERIFY_IS_APPROX(rhs32.noalias() = (s1*m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate(),
84  rhs13 = (s1*m1.adjoint()) * (s2*rhs3).conjugate());
85 
86 
87  m2 = m1.template triangularView<Upper>(); rhs13 = rhs12;
88  VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<Lower>() * (s2*rhs3).conjugate()),
89  rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate());
90 
91  m2 = m1.template triangularView<Lower>();
92  VERIFY_IS_APPROX(rhs22 = (rhs2) * (m2).template selfadjointView<Lower>(), rhs23 = (rhs2) * (m1));
93  VERIFY_IS_APPROX(rhs22 = (s2*rhs2) * (s1*m2).template selfadjointView<Lower>(), rhs23 = (s2*rhs2) * (s1*m1));
94 
95  // destination with a non-default inner-stride
96  // see bug 1741
97  {
98  typedef Matrix<Scalar,Dynamic,Dynamic> MatrixX;
99  MatrixX buffer(2*cols,2*othersize);
100  Map<Rhs1,0,Stride<Dynamic,2> > map1(buffer.data(),cols,othersize,Stride<Dynamic,2>(2*rows,2));
101  buffer.setZero();
102  VERIFY_IS_APPROX( map1.noalias() = (s1*m2).template selfadjointView<Lower>() * (s2*rhs1),
103  rhs13 = (s1*m1) * (s2*rhs1));
104 
105  Map<Rhs2,0,Stride<Dynamic,2> > map2(buffer.data(),rhs22.rows(),rhs22.cols(),Stride<Dynamic,2>(2*rhs22.outerStride(),2));
106  buffer.setZero();
107  VERIFY_IS_APPROX(map2 = (rhs2) * (m2).template selfadjointView<Lower>(), rhs23 = (rhs2) * (m1));
108  }
109 }
110 
111 EIGEN_DECLARE_TEST(product_symm)
112 {
113  for(int i = 0; i < g_repeat ; i++)
114  {
115  CALL_SUBTEST_1(( symm<float,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
116  CALL_SUBTEST_2(( symm<double,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
117  CALL_SUBTEST_3(( symm<std::complex<float>,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2),internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2)) ));
118  CALL_SUBTEST_4(( symm<std::complex<double>,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2),internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2)) ));
119 
120  CALL_SUBTEST_5(( symm<float,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
121  CALL_SUBTEST_6(( symm<double,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
122  CALL_SUBTEST_7(( symm<std::complex<float>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
123  CALL_SUBTEST_8(( symm<std::complex<double>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE)) ));
124  }
125 }
symm
void symm(int size=Size, int othersize=OtherSize)
Definition: product_symm.cpp:12
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
adjoint
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
m2
MatrixType m2(n_dims)
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
Eigen::Stride
Holds strides information for Map.
Definition: Stride.h:48
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
conjugate
EIGEN_DEVICE_FUNC ConjugateReturnType conjugate() const
Definition: CommonCwiseUnaryOps.h:74
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(product_symm)
Definition: product_symm.cpp:111
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
Eigen::RowMajor
Definition: Constants.h:321
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:19