1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 20013 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 // This unit test cannot be easily written to work with EIGEN_DEFAULT_TO_ROW_MAJOR
11 #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
12 #undef EIGEN_DEFAULT_TO_ROW_MAJOR
13 #endif
14
15 #define TEST_ENABLE_TEMPORARY_TRACKING
16
17 #include "main.h"
18
19 // test Ref.h
20
21 // Deal with i387 extended precision
22 #if EIGEN_ARCH_i386 && !(EIGEN_ARCH_x86_64)
23
24 #if EIGEN_COMP_GNUC_STRICT && EIGEN_GNUC_AT_LEAST(4,4)
25 #pragma GCC optimize ("-ffloat-store")
26 #else
27 #undef VERIFY_IS_EQUAL
28 #define VERIFY_IS_EQUAL(X,Y) VERIFY_IS_APPROX(X,Y)
29 #endif
30
31 #endif
32
ref_matrix(const MatrixType & m)33 template<typename MatrixType> void ref_matrix(const MatrixType& m)
34 {
35 typedef typename MatrixType::Index Index;
36 typedef typename MatrixType::Scalar Scalar;
37 typedef typename MatrixType::RealScalar RealScalar;
38 typedef Matrix<Scalar,Dynamic,Dynamic,MatrixType::Options> DynMatrixType;
39 typedef Matrix<RealScalar,Dynamic,Dynamic,MatrixType::Options> RealDynMatrixType;
40
41 typedef Ref<MatrixType> RefMat;
42 typedef Ref<DynMatrixType> RefDynMat;
43 typedef Ref<const DynMatrixType> ConstRefDynMat;
44 typedef Ref<RealDynMatrixType , 0, Stride<Dynamic,Dynamic> > RefRealMatWithStride;
45
46 Index rows = m.rows(), cols = m.cols();
47
48 MatrixType m1 = MatrixType::Random(rows, cols),
49 m2 = m1;
50
51 Index i = internal::random<Index>(0,rows-1);
52 Index j = internal::random<Index>(0,cols-1);
53 Index brows = internal::random<Index>(1,rows-i);
54 Index bcols = internal::random<Index>(1,cols-j);
55
56 RefMat rm0 = m1;
57 VERIFY_IS_EQUAL(rm0, m1);
58 RefDynMat rm1 = m1;
59 VERIFY_IS_EQUAL(rm1, m1);
60 RefDynMat rm2 = m1.block(i,j,brows,bcols);
61 VERIFY_IS_EQUAL(rm2, m1.block(i,j,brows,bcols));
62 rm2.setOnes();
63 m2.block(i,j,brows,bcols).setOnes();
64 VERIFY_IS_EQUAL(m1, m2);
65
66 m2.block(i,j,brows,bcols).setRandom();
67 rm2 = m2.block(i,j,brows,bcols);
68 VERIFY_IS_EQUAL(m1, m2);
69
70 ConstRefDynMat rm3 = m1.block(i,j,brows,bcols);
71 m1.block(i,j,brows,bcols) *= 2;
72 m2.block(i,j,brows,bcols) *= 2;
73 VERIFY_IS_EQUAL(rm3, m2.block(i,j,brows,bcols));
74 RefRealMatWithStride rm4 = m1.real();
75 VERIFY_IS_EQUAL(rm4, m2.real());
76 rm4.array() += 1;
77 m2.real().array() += 1;
78 VERIFY_IS_EQUAL(m1, m2);
79 }
80
ref_vector(const VectorType & m)81 template<typename VectorType> void ref_vector(const VectorType& m)
82 {
83 typedef typename VectorType::Index Index;
84 typedef typename VectorType::Scalar Scalar;
85 typedef typename VectorType::RealScalar RealScalar;
86 typedef Matrix<Scalar,Dynamic,1,VectorType::Options> DynMatrixType;
87 typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixType;
88 typedef Matrix<RealScalar,Dynamic,1,VectorType::Options> RealDynMatrixType;
89
90 typedef Ref<VectorType> RefMat;
91 typedef Ref<DynMatrixType> RefDynMat;
92 typedef Ref<const DynMatrixType> ConstRefDynMat;
93 typedef Ref<RealDynMatrixType , 0, InnerStride<> > RefRealMatWithStride;
94 typedef Ref<DynMatrixType , 0, InnerStride<> > RefMatWithStride;
95
96 Index size = m.size();
97
98 VectorType v1 = VectorType::Random(size),
99 v2 = v1;
100 MatrixType mat1 = MatrixType::Random(size,size),
101 mat2 = mat1,
102 mat3 = MatrixType::Random(size,size);
103
104 Index i = internal::random<Index>(0,size-1);
105 Index bsize = internal::random<Index>(1,size-i);
106
107 RefMat rm0 = v1;
108 VERIFY_IS_EQUAL(rm0, v1);
109 RefDynMat rv1 = v1;
110 VERIFY_IS_EQUAL(rv1, v1);
111 RefDynMat rv2 = v1.segment(i,bsize);
112 VERIFY_IS_EQUAL(rv2, v1.segment(i,bsize));
113 rv2.setOnes();
114 v2.segment(i,bsize).setOnes();
115 VERIFY_IS_EQUAL(v1, v2);
116
117 v2.segment(i,bsize).setRandom();
118 rv2 = v2.segment(i,bsize);
119 VERIFY_IS_EQUAL(v1, v2);
120
121 ConstRefDynMat rm3 = v1.segment(i,bsize);
122 v1.segment(i,bsize) *= 2;
123 v2.segment(i,bsize) *= 2;
124 VERIFY_IS_EQUAL(rm3, v2.segment(i,bsize));
125
126 RefRealMatWithStride rm4 = v1.real();
127 VERIFY_IS_EQUAL(rm4, v2.real());
128 rm4.array() += 1;
129 v2.real().array() += 1;
130 VERIFY_IS_EQUAL(v1, v2);
131
132 RefMatWithStride rm5 = mat1.row(i).transpose();
133 VERIFY_IS_EQUAL(rm5, mat1.row(i).transpose());
134 rm5.array() += 1;
135 mat2.row(i).array() += 1;
136 VERIFY_IS_EQUAL(mat1, mat2);
137 rm5.noalias() = rm4.transpose() * mat3;
138 mat2.row(i) = v2.real().transpose() * mat3;
139 VERIFY_IS_APPROX(mat1, mat2);
140 }
141
check_const_correctness(const PlainObjectType &)142 template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
143 {
144 // verify that ref-to-const don't have LvalueBit
145 typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
146 VERIFY( !(internal::traits<Ref<ConstPlainObjectType> >::Flags & LvalueBit) );
147 VERIFY( !(internal::traits<Ref<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
148 VERIFY( !(Ref<ConstPlainObjectType>::Flags & LvalueBit) );
149 VERIFY( !(Ref<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
150 }
151
152 template<typename B>
call_ref_1(Ref<VectorXf> a,const B & b)153 EIGEN_DONT_INLINE void call_ref_1(Ref<VectorXf> a, const B &b) { VERIFY_IS_EQUAL(a,b); }
154 template<typename B>
call_ref_2(const Ref<const VectorXf> & a,const B & b)155 EIGEN_DONT_INLINE void call_ref_2(const Ref<const VectorXf>& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
156 template<typename B>
call_ref_3(Ref<VectorXf,0,InnerStride<>> a,const B & b)157 EIGEN_DONT_INLINE void call_ref_3(Ref<VectorXf,0,InnerStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
158 template<typename B>
call_ref_4(const Ref<const VectorXf,0,InnerStride<>> & a,const B & b)159 EIGEN_DONT_INLINE void call_ref_4(const Ref<const VectorXf,0,InnerStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
160 template<typename B>
call_ref_5(Ref<MatrixXf,0,OuterStride<>> a,const B & b)161 EIGEN_DONT_INLINE void call_ref_5(Ref<MatrixXf,0,OuterStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
162 template<typename B>
call_ref_6(const Ref<const MatrixXf,0,OuterStride<>> & a,const B & b)163 EIGEN_DONT_INLINE void call_ref_6(const Ref<const MatrixXf,0,OuterStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
164 template<typename B>
call_ref_7(Ref<Matrix<float,Dynamic,3>> a,const B & b)165 EIGEN_DONT_INLINE void call_ref_7(Ref<Matrix<float,Dynamic,3> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
166
call_ref()167 void call_ref()
168 {
169 VectorXcf ca = VectorXcf::Random(10);
170 VectorXf a = VectorXf::Random(10);
171 RowVectorXf b = RowVectorXf::Random(10);
172 MatrixXf A = MatrixXf::Random(10,10);
173 RowVector3f c = RowVector3f::Random();
174 const VectorXf& ac(a);
175 VectorBlock<VectorXf> ab(a,0,3);
176 const VectorBlock<VectorXf> abc(a,0,3);
177
178
179 VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0);
180 VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0);
181 // call_ref_1(ac,a<c); // does not compile because ac is const
182 VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0);
183 VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0);
184 VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0);
185 VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0);
186 // call_ref_1(A.row(3),A.row(3)); // does not compile because innerstride!=1
187 VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0);
188 VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0);
189 // call_ref_1(a+a, a+a); // does not compile for obvious reason
190
191 MatrixXf tmp = A*A.col(1);
192 VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1); // evaluated into a temp
193 VERIFY_EVALUATION_COUNT( call_ref_2(ac.head(5),ac.head(5)), 0);
194 VERIFY_EVALUATION_COUNT( call_ref_2(ac,ac), 0);
195 VERIFY_EVALUATION_COUNT( call_ref_2(a,a), 0);
196 VERIFY_EVALUATION_COUNT( call_ref_2(ab,ab), 0);
197 VERIFY_EVALUATION_COUNT( call_ref_2(a.head(4),a.head(4)), 0);
198 tmp = a+a;
199 VERIFY_EVALUATION_COUNT( call_ref_2(a+a,tmp), 1); // evaluated into a temp
200 VERIFY_EVALUATION_COUNT( call_ref_2(ca.imag(),ca.imag()), 1); // evaluated into a temp
201
202 VERIFY_EVALUATION_COUNT( call_ref_4(ac.head(5),ac.head(5)), 0);
203 tmp = a+a;
204 VERIFY_EVALUATION_COUNT( call_ref_4(a+a,tmp), 1); // evaluated into a temp
205 VERIFY_EVALUATION_COUNT( call_ref_4(ca.imag(),ca.imag()), 0);
206
207 VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0);
208 VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0);
209 VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0);
210 // call_ref_5(A.transpose(),A.transpose()); // does not compile because storage order does not match
211 VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
212 VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0); // storage order do not match, but this is a degenerate case that should work
213 VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0);
214
215 VERIFY_EVALUATION_COUNT( call_ref_6(a,a), 0);
216 VERIFY_EVALUATION_COUNT( call_ref_6(a.head(3),a.head(3)), 0);
217 VERIFY_EVALUATION_COUNT( call_ref_6(A.row(3),A.row(3)), 1); // evaluated into a temp thouth it could be avoided by viewing it as a 1xn matrix
218 tmp = A+A;
219 VERIFY_EVALUATION_COUNT( call_ref_6(A+A,tmp), 1); // evaluated into a temp
220 VERIFY_EVALUATION_COUNT( call_ref_6(A,A), 0);
221 VERIFY_EVALUATION_COUNT( call_ref_6(A.transpose(),A.transpose()), 1); // evaluated into a temp because the storage orders do not match
222 VERIFY_EVALUATION_COUNT( call_ref_6(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
223
224 VERIFY_EVALUATION_COUNT( call_ref_7(c,c), 0);
225 }
226
227 typedef Matrix<double,Dynamic,Dynamic,RowMajor> RowMatrixXd;
test_ref_overload_fun1(Ref<MatrixXd>)228 int test_ref_overload_fun1(Ref<MatrixXd> ) { return 1; }
test_ref_overload_fun1(Ref<RowMatrixXd>)229 int test_ref_overload_fun1(Ref<RowMatrixXd> ) { return 2; }
test_ref_overload_fun1(Ref<MatrixXf>)230 int test_ref_overload_fun1(Ref<MatrixXf> ) { return 3; }
231
test_ref_overload_fun2(Ref<const MatrixXd>)232 int test_ref_overload_fun2(Ref<const MatrixXd> ) { return 4; }
test_ref_overload_fun2(Ref<const MatrixXf>)233 int test_ref_overload_fun2(Ref<const MatrixXf> ) { return 5; }
234
test_ref_ambiguous(const Ref<const ArrayXd> & A,Ref<ArrayXd> B)235 void test_ref_ambiguous(const Ref<const ArrayXd> &A, Ref<ArrayXd> B)
236 {
237 B = A;
238 B = A - A;
239 }
240
241 // See also bug 969
test_ref_overloads()242 void test_ref_overloads()
243 {
244 MatrixXd Ad, Bd;
245 RowMatrixXd rAd, rBd;
246 VERIFY( test_ref_overload_fun1(Ad)==1 );
247 VERIFY( test_ref_overload_fun1(rAd)==2 );
248
249 MatrixXf Af, Bf;
250 VERIFY( test_ref_overload_fun2(Ad)==4 );
251 VERIFY( test_ref_overload_fun2(Ad+Bd)==4 );
252 VERIFY( test_ref_overload_fun2(Af+Bf)==5 );
253
254 ArrayXd A, B;
255 test_ref_ambiguous(A, B);
256 }
257
test_ref()258 void test_ref()
259 {
260 for(int i = 0; i < g_repeat; i++) {
261 CALL_SUBTEST_1( ref_vector(Matrix<float, 1, 1>()) );
262 CALL_SUBTEST_1( check_const_correctness(Matrix<float, 1, 1>()) );
263 CALL_SUBTEST_2( ref_vector(Vector4d()) );
264 CALL_SUBTEST_2( check_const_correctness(Matrix4d()) );
265 CALL_SUBTEST_3( ref_vector(Vector4cf()) );
266 CALL_SUBTEST_4( ref_vector(VectorXcf(8)) );
267 CALL_SUBTEST_5( ref_vector(VectorXi(12)) );
268 CALL_SUBTEST_5( check_const_correctness(VectorXi(12)) );
269
270 CALL_SUBTEST_1( ref_matrix(Matrix<float, 1, 1>()) );
271 CALL_SUBTEST_2( ref_matrix(Matrix4d()) );
272 CALL_SUBTEST_1( ref_matrix(Matrix<float,3,5>()) );
273 CALL_SUBTEST_4( ref_matrix(MatrixXcf(internal::random<int>(1,10),internal::random<int>(1,10))) );
274 CALL_SUBTEST_4( ref_matrix(Matrix<std::complex<double>,10,15>()) );
275 CALL_SUBTEST_5( ref_matrix(MatrixXi(internal::random<int>(1,10),internal::random<int>(1,10))) );
276 CALL_SUBTEST_6( call_ref() );
277 }
278
279 CALL_SUBTEST_7( test_ref_overloads() );
280 }
281