redux.cpp (8239B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> 5 // Copyright (C) 2015 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 #define TEST_ENABLE_TEMPORARY_TRACKING 12 #define EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 8 13 // ^^ see bug 1449 14 15 #include "main.h" 16 17 template<typename MatrixType> void matrixRedux(const MatrixType& m) 18 { 19 typedef typename MatrixType::Scalar Scalar; 20 typedef typename MatrixType::RealScalar RealScalar; 21 22 Index rows = m.rows(); 23 Index cols = m.cols(); 24 25 MatrixType m1 = MatrixType::Random(rows, cols); 26 27 // The entries of m1 are uniformly distributed in [0,1], so m1.prod() is very small. This may lead to test 28 // failures if we underflow into denormals. Thus, we scale so that entries are close to 1. 29 MatrixType m1_for_prod = MatrixType::Ones(rows, cols) + RealScalar(0.2) * m1; 30 31 Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> m2(rows,rows); 32 m2.setRandom(); 33 34 VERIFY_IS_MUCH_SMALLER_THAN(MatrixType::Zero(rows, cols).sum(), Scalar(1)); 35 VERIFY_IS_APPROX(MatrixType::Ones(rows, cols).sum(), Scalar(float(rows*cols))); // the float() here to shut up excessive MSVC warning about int->complex conversion being lossy 36 Scalar s(0), p(1), minc(numext::real(m1.coeff(0))), maxc(numext::real(m1.coeff(0))); 37 for(int j = 0; j < cols; j++) 38 for(int i = 0; i < rows; i++) 39 { 40 s += m1(i,j); 41 p *= m1_for_prod(i,j); 42 minc = (std::min)(numext::real(minc), numext::real(m1(i,j))); 43 maxc = (std::max)(numext::real(maxc), numext::real(m1(i,j))); 44 } 45 const Scalar mean = s/Scalar(RealScalar(rows*cols)); 46 47 VERIFY_IS_APPROX(m1.sum(), s); 48 VERIFY_IS_APPROX(m1.mean(), mean); 49 VERIFY_IS_APPROX(m1_for_prod.prod(), p); 50 VERIFY_IS_APPROX(m1.real().minCoeff(), numext::real(minc)); 51 VERIFY_IS_APPROX(m1.real().maxCoeff(), numext::real(maxc)); 52 53 // test that partial reduction works if nested expressions is forced to evaluate early 54 VERIFY_IS_APPROX((m1.matrix() * m1.matrix().transpose()) .cwiseProduct(m2.matrix()).rowwise().sum().sum(), 55 (m1.matrix() * m1.matrix().transpose()).eval().cwiseProduct(m2.matrix()).rowwise().sum().sum()); 56 57 // test slice vectorization assuming assign is ok 58 Index r0 = internal::random<Index>(0,rows-1); 59 Index c0 = internal::random<Index>(0,cols-1); 60 Index r1 = internal::random<Index>(r0+1,rows)-r0; 61 Index c1 = internal::random<Index>(c0+1,cols)-c0; 62 VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).sum(), m1.block(r0,c0,r1,c1).eval().sum()); 63 VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).mean(), m1.block(r0,c0,r1,c1).eval().mean()); 64 VERIFY_IS_APPROX(m1_for_prod.block(r0,c0,r1,c1).prod(), m1_for_prod.block(r0,c0,r1,c1).eval().prod()); 65 VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).real().minCoeff(), m1.block(r0,c0,r1,c1).real().eval().minCoeff()); 66 VERIFY_IS_APPROX(m1.block(r0,c0,r1,c1).real().maxCoeff(), m1.block(r0,c0,r1,c1).real().eval().maxCoeff()); 67 68 // regression for bug 1090 69 const int R1 = MatrixType::RowsAtCompileTime>=2 ? MatrixType::RowsAtCompileTime/2 : 6; 70 const int C1 = MatrixType::ColsAtCompileTime>=2 ? MatrixType::ColsAtCompileTime/2 : 6; 71 if(R1<=rows-r0 && C1<=cols-c0) 72 { 73 VERIFY_IS_APPROX( (m1.template block<R1,C1>(r0,c0).sum()), m1.block(r0,c0,R1,C1).sum() ); 74 } 75 76 // test empty objects 77 VERIFY_IS_APPROX(m1.block(r0,c0,0,0).sum(), Scalar(0)); 78 VERIFY_IS_APPROX(m1.block(r0,c0,0,0).prod(), Scalar(1)); 79 80 // test nesting complex expression 81 VERIFY_EVALUATION_COUNT( (m1.matrix()*m1.matrix().transpose()).sum(), (MatrixType::IsVectorAtCompileTime && MatrixType::SizeAtCompileTime!=1 ? 0 : 1) ); 82 VERIFY_EVALUATION_COUNT( ((m1.matrix()*m1.matrix().transpose())+m2).sum(),(MatrixType::IsVectorAtCompileTime && MatrixType::SizeAtCompileTime!=1 ? 0 : 1)); 83 } 84 85 template<typename VectorType> void vectorRedux(const VectorType& w) 86 { 87 using std::abs; 88 typedef typename VectorType::Scalar Scalar; 89 typedef typename NumTraits<Scalar>::Real RealScalar; 90 Index size = w.size(); 91 92 VectorType v = VectorType::Random(size); 93 VectorType v_for_prod = VectorType::Ones(size) + Scalar(0.2) * v; // see comment above declaration of m1_for_prod 94 95 for(int i = 1; i < size; i++) 96 { 97 Scalar s(0), p(1); 98 RealScalar minc(numext::real(v.coeff(0))), maxc(numext::real(v.coeff(0))); 99 for(int j = 0; j < i; j++) 100 { 101 s += v[j]; 102 p *= v_for_prod[j]; 103 minc = (std::min)(minc, numext::real(v[j])); 104 maxc = (std::max)(maxc, numext::real(v[j])); 105 } 106 VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.head(i).sum()), Scalar(1)); 107 VERIFY_IS_APPROX(p, v_for_prod.head(i).prod()); 108 VERIFY_IS_APPROX(minc, v.real().head(i).minCoeff()); 109 VERIFY_IS_APPROX(maxc, v.real().head(i).maxCoeff()); 110 } 111 112 for(int i = 0; i < size-1; i++) 113 { 114 Scalar s(0), p(1); 115 RealScalar minc(numext::real(v.coeff(i))), maxc(numext::real(v.coeff(i))); 116 for(int j = i; j < size; j++) 117 { 118 s += v[j]; 119 p *= v_for_prod[j]; 120 minc = (std::min)(minc, numext::real(v[j])); 121 maxc = (std::max)(maxc, numext::real(v[j])); 122 } 123 VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.tail(size-i).sum()), Scalar(1)); 124 VERIFY_IS_APPROX(p, v_for_prod.tail(size-i).prod()); 125 VERIFY_IS_APPROX(minc, v.real().tail(size-i).minCoeff()); 126 VERIFY_IS_APPROX(maxc, v.real().tail(size-i).maxCoeff()); 127 } 128 129 for(int i = 0; i < size/2; i++) 130 { 131 Scalar s(0), p(1); 132 RealScalar minc(numext::real(v.coeff(i))), maxc(numext::real(v.coeff(i))); 133 for(int j = i; j < size-i; j++) 134 { 135 s += v[j]; 136 p *= v_for_prod[j]; 137 minc = (std::min)(minc, numext::real(v[j])); 138 maxc = (std::max)(maxc, numext::real(v[j])); 139 } 140 VERIFY_IS_MUCH_SMALLER_THAN(abs(s - v.segment(i, size-2*i).sum()), Scalar(1)); 141 VERIFY_IS_APPROX(p, v_for_prod.segment(i, size-2*i).prod()); 142 VERIFY_IS_APPROX(minc, v.real().segment(i, size-2*i).minCoeff()); 143 VERIFY_IS_APPROX(maxc, v.real().segment(i, size-2*i).maxCoeff()); 144 } 145 146 // test empty objects 147 VERIFY_IS_APPROX(v.head(0).sum(), Scalar(0)); 148 VERIFY_IS_APPROX(v.tail(0).prod(), Scalar(1)); 149 VERIFY_RAISES_ASSERT(v.head(0).mean()); 150 VERIFY_RAISES_ASSERT(v.head(0).minCoeff()); 151 VERIFY_RAISES_ASSERT(v.head(0).maxCoeff()); 152 } 153 154 EIGEN_DECLARE_TEST(redux) 155 { 156 // the max size cannot be too large, otherwise reduxion operations obviously generate large errors. 157 int maxsize = (std::min)(100,EIGEN_TEST_MAX_SIZE); 158 TEST_SET_BUT_UNUSED_VARIABLE(maxsize); 159 for(int i = 0; i < g_repeat; i++) { 160 CALL_SUBTEST_1( matrixRedux(Matrix<float, 1, 1>()) ); 161 CALL_SUBTEST_1( matrixRedux(Array<float, 1, 1>()) ); 162 CALL_SUBTEST_2( matrixRedux(Matrix2f()) ); 163 CALL_SUBTEST_2( matrixRedux(Array2f()) ); 164 CALL_SUBTEST_2( matrixRedux(Array22f()) ); 165 CALL_SUBTEST_3( matrixRedux(Matrix4d()) ); 166 CALL_SUBTEST_3( matrixRedux(Array4d()) ); 167 CALL_SUBTEST_3( matrixRedux(Array44d()) ); 168 CALL_SUBTEST_4( matrixRedux(MatrixXcf(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 169 CALL_SUBTEST_4( matrixRedux(ArrayXXcf(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 170 CALL_SUBTEST_5( matrixRedux(MatrixXd (internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 171 CALL_SUBTEST_5( matrixRedux(ArrayXXd (internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 172 CALL_SUBTEST_6( matrixRedux(MatrixXi (internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 173 CALL_SUBTEST_6( matrixRedux(ArrayXXi (internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) ); 174 } 175 for(int i = 0; i < g_repeat; i++) { 176 CALL_SUBTEST_7( vectorRedux(Vector4f()) ); 177 CALL_SUBTEST_7( vectorRedux(Array4f()) ); 178 CALL_SUBTEST_5( vectorRedux(VectorXd(internal::random<int>(1,maxsize))) ); 179 CALL_SUBTEST_5( vectorRedux(ArrayXd(internal::random<int>(1,maxsize))) ); 180 CALL_SUBTEST_8( vectorRedux(VectorXf(internal::random<int>(1,maxsize))) ); 181 CALL_SUBTEST_8( vectorRedux(ArrayXf(internal::random<int>(1,maxsize))) ); 182 } 183 }