1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
5 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 #include "main.h"
12
13 template<typename MatrixType>
equalsIdentity(const MatrixType & A)14 bool equalsIdentity(const MatrixType& A)
15 {
16 typedef typename MatrixType::Scalar Scalar;
17 Scalar zero = static_cast<Scalar>(0);
18
19 bool offDiagOK = true;
20 for (Index i = 0; i < A.rows(); ++i) {
21 for (Index j = i+1; j < A.cols(); ++j) {
22 offDiagOK = offDiagOK && (A(i,j) == zero);
23 }
24 }
25 for (Index i = 0; i < A.rows(); ++i) {
26 for (Index j = 0; j < (std::min)(i, A.cols()); ++j) {
27 offDiagOK = offDiagOK && (A(i,j) == zero);
28 }
29 }
30
31 bool diagOK = (A.diagonal().array() == 1).all();
32 return offDiagOK && diagOK;
33
34 }
35
36 template<typename VectorType>
check_extremity_accuracy(const VectorType & v,const typename VectorType::Scalar & low,const typename VectorType::Scalar & high)37 void check_extremity_accuracy(const VectorType &v, const typename VectorType::Scalar &low, const typename VectorType::Scalar &high)
38 {
39 typedef typename VectorType::Scalar Scalar;
40 typedef typename VectorType::RealScalar RealScalar;
41
42 RealScalar prec = internal::is_same<RealScalar,float>::value ? NumTraits<RealScalar>::dummy_precision()*10 : NumTraits<RealScalar>::dummy_precision()/10;
43 Index size = v.size();
44
45 if(size<20)
46 return;
47
48 for (int i=0; i<size; ++i)
49 {
50 if(i<5 || i>size-6)
51 {
52 Scalar ref = (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1);
53 if(std::abs(ref)>1)
54 {
55 if(!internal::isApprox(v(i), ref, prec))
56 std::cout << v(i) << " != " << ref << " ; relative error: " << std::abs((v(i)-ref)/ref) << " ; required precision: " << prec << " ; range: " << low << "," << high << " ; i: " << i << "\n";
57 VERIFY(internal::isApprox(v(i), (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1), prec));
58 }
59 }
60 }
61 }
62
63 template<typename VectorType>
testVectorType(const VectorType & base)64 void testVectorType(const VectorType& base)
65 {
66 typedef typename VectorType::Scalar Scalar;
67 typedef typename VectorType::RealScalar RealScalar;
68
69 const Index size = base.size();
70
71 Scalar high = internal::random<Scalar>(-500,500);
72 Scalar low = (size == 1 ? high : internal::random<Scalar>(-500,500));
73 if (low>high) std::swap(low,high);
74
75 // check low==high
76 if(internal::random<float>(0.f,1.f)<0.05f)
77 low = high;
78 // check abs(low) >> abs(high)
79 else if(size>2 && std::numeric_limits<RealScalar>::max_exponent10>0 && internal::random<float>(0.f,1.f)<0.1f)
80 low = -internal::random<Scalar>(1,2) * RealScalar(std::pow(RealScalar(10),std::numeric_limits<RealScalar>::max_exponent10/2));
81
82 const Scalar step = ((size == 1) ? 1 : (high-low)/(size-1));
83
84 // check whether the result yields what we expect it to do
85 VectorType m(base);
86 m.setLinSpaced(size,low,high);
87
88 if(!NumTraits<Scalar>::IsInteger)
89 {
90 VectorType n(size);
91 for (int i=0; i<size; ++i)
92 n(i) = low+i*step;
93 VERIFY_IS_APPROX(m,n);
94
95 CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
96 }
97
98 if((!NumTraits<Scalar>::IsInteger) || ((high-low)>=size && (Index(high-low)%(size-1))==0) || (Index(high-low+1)<size && (size%Index(high-low+1))==0))
99 {
100 VectorType n(size);
101 if((!NumTraits<Scalar>::IsInteger) || (high-low>=size))
102 for (int i=0; i<size; ++i)
103 n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/(size-1));
104 else
105 for (int i=0; i<size; ++i)
106 n(i) = size==1 ? low : low + Scalar((double(high-low+1)*double(i))/double(size));
107 VERIFY_IS_APPROX(m,n);
108
109 // random access version
110 m = VectorType::LinSpaced(size,low,high);
111 VERIFY_IS_APPROX(m,n);
112 VERIFY( internal::isApprox(m(m.size()-1),high) );
113 VERIFY( size==1 || internal::isApprox(m(0),low) );
114 VERIFY_IS_EQUAL(m(m.size()-1) , high);
115 if(!NumTraits<Scalar>::IsInteger)
116 CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
117 }
118
119 VERIFY( m(m.size()-1) <= high );
120 VERIFY( (m.array() <= high).all() );
121 VERIFY( (m.array() >= low).all() );
122
123
124 VERIFY( m(m.size()-1) >= low );
125 if(size>=1)
126 {
127 VERIFY( internal::isApprox(m(0),low) );
128 VERIFY_IS_EQUAL(m(0) , low);
129 }
130
131 // check whether everything works with row and col major vectors
132 Matrix<Scalar,Dynamic,1> row_vector(size);
133 Matrix<Scalar,1,Dynamic> col_vector(size);
134 row_vector.setLinSpaced(size,low,high);
135 col_vector.setLinSpaced(size,low,high);
136 // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit
137 // when computing the squared sum in isApprox, thus the 2x factor.
138 VERIFY( row_vector.isApprox(col_vector.transpose(), Scalar(2)*NumTraits<Scalar>::epsilon()));
139
140 Matrix<Scalar,Dynamic,1> size_changer(size+50);
141 size_changer.setLinSpaced(size,low,high);
142 VERIFY( size_changer.size() == size );
143
144 typedef Matrix<Scalar,1,1> ScalarMatrix;
145 ScalarMatrix scalar;
146 scalar.setLinSpaced(1,low,high);
147 VERIFY_IS_APPROX( scalar, ScalarMatrix::Constant(high) );
148 VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) );
149
150 // regression test for bug 526 (linear vectorized transversal)
151 if (size > 1 && (!NumTraits<Scalar>::IsInteger)) {
152 m.tail(size-1).setLinSpaced(low, high);
153 VERIFY_IS_APPROX(m(size-1), high);
154 }
155
156 // regression test for bug 1383 (LinSpaced with empty size/range)
157 {
158 Index n0 = VectorType::SizeAtCompileTime==Dynamic ? 0 : VectorType::SizeAtCompileTime;
159 low = internal::random<Scalar>();
160 m = VectorType::LinSpaced(n0,low,low-1);
161 VERIFY(m.size()==n0);
162
163 if(VectorType::SizeAtCompileTime==Dynamic)
164 {
165 VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,0,Scalar(n0-1)).sum(),Scalar(0));
166 VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,low,low-1).sum(),Scalar(0));
167 }
168
169 m.setLinSpaced(n0,0,Scalar(n0-1));
170 VERIFY(m.size()==n0);
171 m.setLinSpaced(n0,low,low-1);
172 VERIFY(m.size()==n0);
173
174 // empty range only:
175 VERIFY_IS_APPROX(VectorType::LinSpaced(size,low,low),VectorType::Constant(size,low));
176 m.setLinSpaced(size,low,low);
177 VERIFY_IS_APPROX(m,VectorType::Constant(size,low));
178
179 if(NumTraits<Scalar>::IsInteger)
180 {
181 VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,Scalar(low+size-1)), VectorType::LinSpaced(size,Scalar(low+size-1),low).reverse() );
182
183 if(VectorType::SizeAtCompileTime==Dynamic)
184 {
185 // Check negative multiplicator path:
186 for(Index k=1; k<5; ++k)
187 VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,Scalar(low+(size-1)*k)), VectorType::LinSpaced(size,Scalar(low+(size-1)*k),low).reverse() );
188 // Check negative divisor path:
189 for(Index k=1; k<5; ++k)
190 VERIFY_IS_APPROX( VectorType::LinSpaced(size*k,low,Scalar(low+size-1)), VectorType::LinSpaced(size*k,Scalar(low+size-1),low).reverse() );
191 }
192 }
193 }
194 }
195
196 template<typename MatrixType>
testMatrixType(const MatrixType & m)197 void testMatrixType(const MatrixType& m)
198 {
199 using std::abs;
200 const Index rows = m.rows();
201 const Index cols = m.cols();
202 typedef typename MatrixType::Scalar Scalar;
203 typedef typename MatrixType::RealScalar RealScalar;
204
205 Scalar s1;
206 do {
207 s1 = internal::random<Scalar>();
208 } while(abs(s1)<RealScalar(1e-5) && (!NumTraits<Scalar>::IsInteger));
209
210 MatrixType A;
211 A.setIdentity(rows, cols);
212 VERIFY(equalsIdentity(A));
213 VERIFY(equalsIdentity(MatrixType::Identity(rows, cols)));
214
215
216 A = MatrixType::Constant(rows,cols,s1);
217 Index i = internal::random<Index>(0,rows-1);
218 Index j = internal::random<Index>(0,cols-1);
219 VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1)(i,j), s1 );
220 VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1).coeff(i,j), s1 );
221 VERIFY_IS_APPROX( A(i,j), s1 );
222 }
223
test_nullary()224 void test_nullary()
225 {
226 CALL_SUBTEST_1( testMatrixType(Matrix2d()) );
227 CALL_SUBTEST_2( testMatrixType(MatrixXcf(internal::random<int>(1,300),internal::random<int>(1,300))) );
228 CALL_SUBTEST_3( testMatrixType(MatrixXf(internal::random<int>(1,300),internal::random<int>(1,300))) );
229
230 for(int i = 0; i < g_repeat*10; i++) {
231 CALL_SUBTEST_4( testVectorType(VectorXd(internal::random<int>(1,30000))) );
232 CALL_SUBTEST_5( testVectorType(Vector4d()) ); // regression test for bug 232
233 CALL_SUBTEST_6( testVectorType(Vector3d()) );
234 CALL_SUBTEST_7( testVectorType(VectorXf(internal::random<int>(1,30000))) );
235 CALL_SUBTEST_8( testVectorType(Vector3f()) );
236 CALL_SUBTEST_8( testVectorType(Vector4f()) );
237 CALL_SUBTEST_8( testVectorType(Matrix<float,8,1>()) );
238 CALL_SUBTEST_8( testVectorType(Matrix<float,1,1>()) );
239
240 CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(1,10))) );
241 CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(9,300))) );
242 CALL_SUBTEST_9( testVectorType(Matrix<int,1,1>()) );
243 }
244
245 #ifdef EIGEN_TEST_PART_6
246 // Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
247 VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<double>::epsilon() );
248 #endif
249
250 #ifdef EIGEN_TEST_PART_9
251 // Check possible overflow issue
252 {
253 int n = 60000;
254 ArrayXi a1(n), a2(n);
255 a1.setLinSpaced(n, 0, n-1);
256 for(int i=0; i<n; ++i)
257 a2(i) = i;
258 VERIFY_IS_APPROX(a1,a2);
259 }
260 #endif
261
262 #ifdef EIGEN_TEST_PART_10
263 // check some internal logic
264 VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<double> >::value ));
265 VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<double> >::value ));
266 VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<double> >::value ));
267 VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<double> >::ret ));
268
269 VERIFY(( !internal::has_nullary_operator<internal::scalar_identity_op<double> >::value ));
270 VERIFY(( !internal::has_unary_operator<internal::scalar_identity_op<double> >::value ));
271 VERIFY(( internal::has_binary_operator<internal::scalar_identity_op<double> >::value ));
272 VERIFY(( !internal::functor_has_linear_access<internal::scalar_identity_op<double> >::ret ));
273
274 VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<float,float> >::value ));
275 VERIFY(( internal::has_unary_operator<internal::linspaced_op<float,float> >::value ));
276 VERIFY(( !internal::has_binary_operator<internal::linspaced_op<float,float> >::value ));
277 VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<float,float> >::ret ));
278
279 // Regression unit test for a weird MSVC bug.
280 // Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details.
281 // See also traits<Ref>::match.
282 {
283 MatrixXf A = MatrixXf::Random(3,3);
284 Ref<const MatrixXf> R = 2.0*A;
285 VERIFY_IS_APPROX(R, A+A);
286
287 Ref<const MatrixXf> R1 = MatrixXf::Random(3,3)+A;
288
289 VectorXi V = VectorXi::Random(3);
290 Ref<const VectorXi> R2 = VectorXi::LinSpaced(3,1,3)+V;
291 VERIFY_IS_APPROX(R2, V+Vector3i(1,2,3));
292
293 VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<float> >::value ));
294 VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<float> >::value ));
295 VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<float> >::value ));
296 VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<float> >::ret ));
297
298 VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<int,int> >::value ));
299 VERIFY(( internal::has_unary_operator<internal::linspaced_op<int,int> >::value ));
300 VERIFY(( !internal::has_binary_operator<internal::linspaced_op<int,int> >::value ));
301 VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<int,int> >::ret ));
302 }
303 #endif
304 }
305