product_trsolve.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 #define VERIFY_TRSM(TRI,XB) { \
13  (XB).setRandom(); ref = (XB); \
14  (TRI).solveInPlace(XB); \
15  VERIFY_IS_APPROX((TRI).toDenseMatrix() * (XB), ref); \
16  (XB).setRandom(); ref = (XB); \
17  (XB) = (TRI).solve(XB); \
18  VERIFY_IS_APPROX((TRI).toDenseMatrix() * (XB), ref); \
19  }
20 
21 #define VERIFY_TRSM_ONTHERIGHT(TRI,XB) { \
22  (XB).setRandom(); ref = (XB); \
23  (TRI).transpose().template solveInPlace<OnTheRight>(XB.transpose()); \
24  VERIFY_IS_APPROX((XB).transpose() * (TRI).transpose().toDenseMatrix(), ref.transpose()); \
25  (XB).setRandom(); ref = (XB); \
26  (XB).transpose() = (TRI).transpose().template solve<OnTheRight>(XB.transpose()); \
27  VERIFY_IS_APPROX((XB).transpose() * (TRI).transpose().toDenseMatrix(), ref.transpose()); \
28  }
29 
30 template<typename Scalar,int Size, int Cols> void trsolve(int size=Size,int cols=Cols)
31 {
32  typedef typename NumTraits<Scalar>::Real RealScalar;
33 
34  Matrix<Scalar,Size,Size,ColMajor> cmLhs(size,size);
35  Matrix<Scalar,Size,Size,RowMajor> rmLhs(size,size);
36 
37  enum { colmajor = Size==1 ? RowMajor : ColMajor,
38  rowmajor = Cols==1 ? ColMajor : RowMajor };
39  Matrix<Scalar,Size,Cols,colmajor> cmRhs(size,cols);
40  Matrix<Scalar,Size,Cols,rowmajor> rmRhs(size,cols);
41  Matrix<Scalar,Dynamic,Dynamic,colmajor> ref(size,cols);
42 
43  cmLhs.setRandom(); cmLhs *= static_cast<RealScalar>(0.1); cmLhs.diagonal().array() += static_cast<RealScalar>(1);
44  rmLhs.setRandom(); rmLhs *= static_cast<RealScalar>(0.1); rmLhs.diagonal().array() += static_cast<RealScalar>(1);
45 
46  VERIFY_TRSM(cmLhs.conjugate().template triangularView<Lower>(), cmRhs);
47  VERIFY_TRSM(cmLhs.adjoint() .template triangularView<Lower>(), cmRhs);
48  VERIFY_TRSM(cmLhs .template triangularView<Upper>(), cmRhs);
49  VERIFY_TRSM(cmLhs .template triangularView<Lower>(), rmRhs);
50  VERIFY_TRSM(cmLhs.conjugate().template triangularView<Upper>(), rmRhs);
51  VERIFY_TRSM(cmLhs.adjoint() .template triangularView<Upper>(), rmRhs);
52 
53  VERIFY_TRSM(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs);
54  VERIFY_TRSM(cmLhs .template triangularView<UnitUpper>(), rmRhs);
55 
56  VERIFY_TRSM(rmLhs .template triangularView<Lower>(), cmRhs);
57  VERIFY_TRSM(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs);
58 
59 
60  VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Lower>(), cmRhs);
61  VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Upper>(), cmRhs);
62  VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<Lower>(), rmRhs);
63  VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<Upper>(), rmRhs);
64 
65  VERIFY_TRSM_ONTHERIGHT(cmLhs.conjugate().template triangularView<UnitLower>(), cmRhs);
66  VERIFY_TRSM_ONTHERIGHT(cmLhs .template triangularView<UnitUpper>(), rmRhs);
67 
68  VERIFY_TRSM_ONTHERIGHT(rmLhs .template triangularView<Lower>(), cmRhs);
69  VERIFY_TRSM_ONTHERIGHT(rmLhs.conjugate().template triangularView<UnitUpper>(), rmRhs);
70 
71  int c = internal::random<int>(0,cols-1);
72  VERIFY_TRSM(rmLhs.template triangularView<Lower>(), rmRhs.col(c));
73  VERIFY_TRSM(cmLhs.template triangularView<Lower>(), rmRhs.col(c));
74 
75  // destination with a non-default inner-stride
76  // see bug 1741
77  {
78  typedef Matrix<Scalar,Dynamic,Dynamic> MatrixX;
79  MatrixX buffer(2*cmRhs.rows(),2*cmRhs.cols());
80  Map<Matrix<Scalar,Size,Cols,colmajor>,0,Stride<Dynamic,2> > map1(buffer.data(),cmRhs.rows(),cmRhs.cols(),Stride<Dynamic,2>(2*cmRhs.outerStride(),2));
81  Map<Matrix<Scalar,Size,Cols,rowmajor>,0,Stride<Dynamic,2> > map2(buffer.data(),rmRhs.rows(),rmRhs.cols(),Stride<Dynamic,2>(2*rmRhs.outerStride(),2));
82  buffer.setZero();
83  VERIFY_TRSM(cmLhs.conjugate().template triangularView<Lower>(), map1);
84  buffer.setZero();
85  VERIFY_TRSM(cmLhs .template triangularView<Lower>(), map2);
86  }
87 
88  if(Size==Dynamic)
89  {
90  cmLhs.resize(0,0);
91  cmRhs.resize(0,cmRhs.cols());
92  Matrix<Scalar,Size,Cols,colmajor> res = cmLhs.template triangularView<Lower>().solve(cmRhs);
93  VERIFY_IS_EQUAL(res.rows(),0);
94  VERIFY_IS_EQUAL(res.cols(),cmRhs.cols());
95  res = cmRhs;
96  cmLhs.template triangularView<Lower>().solveInPlace(res);
97  VERIFY_IS_EQUAL(res.rows(),0);
98  VERIFY_IS_EQUAL(res.cols(),cmRhs.cols());
99  }
100 }
101 
102 EIGEN_DECLARE_TEST(product_trsolve)
103 {
104  for(int i = 0; i < g_repeat ; i++)
105  {
106  // matrices
107  CALL_SUBTEST_1((trsolve<float,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
108  CALL_SUBTEST_2((trsolve<double,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE),internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
109  CALL_SUBTEST_3((trsolve<std::complex<float>,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2),internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))));
110  CALL_SUBTEST_4((trsolve<std::complex<double>,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2),internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))));
111 
112  // vectors
113  CALL_SUBTEST_5((trsolve<float,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
114  CALL_SUBTEST_6((trsolve<double,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
115  CALL_SUBTEST_7((trsolve<std::complex<float>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
116  CALL_SUBTEST_8((trsolve<std::complex<double>,Dynamic,1>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))));
117 
118  // meta-unrollers
119  CALL_SUBTEST_9((trsolve<float,4,1>()));
120  CALL_SUBTEST_10((trsolve<double,4,1>()));
121  CALL_SUBTEST_11((trsolve<std::complex<float>,4,1>()));
122  CALL_SUBTEST_12((trsolve<float,1,1>()));
123  CALL_SUBTEST_13((trsolve<float,1,2>()));
124  CALL_SUBTEST_14((trsolve<float,3,1>()));
125 
126  }
127 }
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setRandom
Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setRandom(Index size)
Definition: Random.h:151
Eigen::Stride
Holds strides information for Map.
Definition: Stride.h:48
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Eigen::Matrix::outerStride
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index outerStride() const EIGEN_NOEXCEPT
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:429
VERIFY_TRSM_ONTHERIGHT
#define VERIFY_TRSM_ONTHERIGHT(TRI, XB)
Definition: product_trsolve.cpp:21
buffer
Definition: pytypes.h:2270
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
CALL_SUBTEST_11
#define CALL_SUBTEST_11(FUNC)
Definition: split_test_helper.h:64
CALL_SUBTEST_9
#define CALL_SUBTEST_9(FUNC)
Definition: split_test_helper.h:52
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
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
CALL_SUBTEST_13
#define CALL_SUBTEST_13(FUNC)
Definition: split_test_helper.h:76
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(product_trsolve)
Definition: product_trsolve.cpp:102
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
CALL_SUBTEST_10
#define CALL_SUBTEST_10(FUNC)
Definition: split_test_helper.h:58
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::rows
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: PlainObjectBase.h:143
trsolve
void trsolve(int size=Size, int cols=Cols)
Definition: product_trsolve.cpp:30
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
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
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::cols
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: PlainObjectBase.h:145
main.h
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
CALL_SUBTEST_12
#define CALL_SUBTEST_12(FUNC)
Definition: split_test_helper.h:70
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:319
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
CALL_SUBTEST_14
#define CALL_SUBTEST_14(FUNC)
Definition: split_test_helper.h:82
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
VERIFY_TRSM
#define VERIFY_TRSM(TRI, XB)
Definition: product_trsolve.cpp:12
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Cols
static const int Cols
Definition: testCallRecord.cpp:32
test_eigen.ref
ref
Definition: test_eigen.py:9

gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:34:33