1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename T> bool isNotNaN(const T& x) 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return x==x; 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// workaround aggressive optimization in ICC 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename T> bool isFinite(const T& x) 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return isNotNaN(sub(x,x)); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename T> EIGEN_DONT_INLINE T copy(const T& x) 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return x; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void stable_norm(const MatrixType& m) 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* this test covers the following files: 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath StableNorm.h 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check the basic machine-dependent constants. 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int ibeta, it, iemin, iemax; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY( (!(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) || (it<=4 && ibeta <= 3 ) || it<2)) 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "the stable norm algorithm cannot be guaranteed on this computer"); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4)); 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4)); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType vzero = MatrixType::Zero(rows, cols), 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vrand = MatrixType::Random(rows, cols), 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vbig(rows, cols), 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vsmall(rows,cols); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vbig.fill(big); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vsmall.fill(small); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm()); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm()); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm()); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar size = static_cast<RealScalar>(m.size()); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test isFinite 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(!isFinite( std::numeric_limits<RealScalar>::infinity())); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(!isFinite(internal::sqrt(-internal::abs(big)))); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test overflow 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(isFinite(internal::sqrt(size)*internal::abs(big))); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_NOT_APPROX(internal::sqrt(copy(vbig.squaredNorm())), internal::abs(internal::sqrt(size)*big)); // here the default norm must fail 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vbig.stableNorm(), internal::sqrt(size)*internal::abs(big)); 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vbig.blueNorm(), internal::sqrt(size)*internal::abs(big)); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vbig.hypotNorm(), internal::sqrt(size)*internal::abs(big)); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test underflow 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(isFinite(internal::sqrt(size)*internal::abs(small))); 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_NOT_APPROX(internal::sqrt(copy(vsmall.squaredNorm())), internal::abs(internal::sqrt(size)*small)); // here the default norm must fail 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vsmall.stableNorm(), internal::sqrt(size)*internal::abs(small)); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vsmall.blueNorm(), internal::sqrt(size)*internal::abs(small)); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vsmall.hypotNorm(), internal::sqrt(size)*internal::abs(small)); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Test compilation of cwise() version 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm()); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm()); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm()); 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm()); 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm()); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm()); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_stable_norm() 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1( stable_norm(Matrix<float, 1, 1>()) ); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( stable_norm(Vector4d()) ); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3( stable_norm(VectorXd(internal::random<int>(10,2000))) ); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4( stable_norm(VectorXf(internal::random<int>(10,2000))) ); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_5( stable_norm(VectorXcd(internal::random<int>(10,2000))) ); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 111