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 (at) maths.leeds.ac.uk> 5 // Copyright (C) 2016 Gael Guennebaud <gael.guennebaud (at) 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> 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> 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> 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> 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 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
