Functions
stable_norm.cpp File Reference
#include "main.h"

Functions

template<typename T >
EIGEN_DONT_INLINE T copy (const T &x)
 
template<typename MatrixType >
void stable_norm (const MatrixType &m)
 
void test_empty ()
 
template<typename Scalar >
void test_hypot ()
 
 EIGEN_DECLARE_TEST (stable_norm)
 

Function Documentation

◆ copy()

template<typename T >
EIGEN_DONT_INLINE T copy ( const T &  x)
13  {
14  return x;
15 }
plotDoE.x
list x
Definition: plotDoE.py:28

References plotDoE::x.

Referenced by stable_norm().

◆ EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST ( stable_norm  )
242  {
243  CALL_SUBTEST_1(test_empty());
244 
245  for (int i = 0; i < g_repeat; i++) {
246  CALL_SUBTEST_3(test_hypot<double>());
247  CALL_SUBTEST_4(test_hypot<float>());
248  CALL_SUBTEST_5(test_hypot<std::complex<double> >());
249  CALL_SUBTEST_6(test_hypot<std::complex<float> >());
250 
251  CALL_SUBTEST_1(stable_norm(Matrix<float, 1, 1>()));
252  CALL_SUBTEST_2(stable_norm(Vector4d()));
253  CALL_SUBTEST_3(stable_norm(VectorXd(internal::random<int>(10, 2000))));
254  CALL_SUBTEST_3(stable_norm(MatrixXd(internal::random<int>(10, 200), internal::random<int>(10, 200))));
255  CALL_SUBTEST_4(stable_norm(VectorXf(internal::random<int>(10, 2000))));
256  CALL_SUBTEST_5(stable_norm(VectorXcd(internal::random<int>(10, 2000))));
257  CALL_SUBTEST_6(stable_norm(VectorXcf(internal::random<int>(10, 2000))));
258  }
259 }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: Eigen/Eigen/src/Core/Matrix.h:186
Eigen::g_repeat
static int g_repeat
Definition: main.h:191
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
test_hypot
void test_hypot()
Definition: stable_norm.cpp:218
test_empty
void test_empty()
Definition: stable_norm.cpp:212
stable_norm
void stable_norm(const MatrixType &m)
Definition: stable_norm.cpp:18

References CALL_SUBTEST_1, CALL_SUBTEST_2, CALL_SUBTEST_3, CALL_SUBTEST_4, CALL_SUBTEST_5, CALL_SUBTEST_6, Eigen::g_repeat, i, stable_norm(), test_empty(), and test_hypot().

◆ stable_norm()

template<typename MatrixType >
void stable_norm ( const MatrixType &  m)
18  {
19  /* this test covers the following files:
20  StableNorm.h
21  */
22  using std::abs;
23  using std::sqrt;
24  typedef typename MatrixType::Scalar Scalar;
25  typedef typename NumTraits<Scalar>::Real RealScalar;
26 
27  bool complex_real_product_ok = true;
28 
29  // Check the basic machine-dependent constants.
30  {
31  int ibeta, it, iemin, iemax;
32 
33  ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers
34  it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa
35  iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent
36  iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent
37 
38  VERIFY((!(iemin > 1 - 2 * it || 1 + it > iemax || (it == 2 && ibeta < 5) || (it <= 4 && ibeta <= 3) || it < 2)) &&
39  "the stable norm algorithm cannot be guaranteed on this computer");
40 
41  Scalar inf = std::numeric_limits<RealScalar>::infinity();
42  if (NumTraits<Scalar>::IsComplex && (numext::isnan)(inf * RealScalar(1))) {
43  complex_real_product_ok = false;
44  static bool first = true;
45  if (first)
46  std::cerr << "WARNING: compiler mess up complex*real product, " << inf << " * " << 1.0 << " = "
47  << inf * RealScalar(1) << std::endl;
48  first = false;
49  }
50  }
51 
52  Index rows = m.rows();
53  Index cols = m.cols();
54 
55  // get a non-zero random factor
56  Scalar factor = internal::random<Scalar>();
57  while (numext::abs2(factor) < RealScalar(1e-4)) factor = internal::random<Scalar>();
58  Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
59 
60  factor = internal::random<Scalar>();
61  while (numext::abs2(factor) < RealScalar(1e-4)) factor = internal::random<Scalar>();
62  Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
63 
64  Scalar one(1);
65 
66  MatrixType vzero = MatrixType::Zero(rows, cols), vrand = MatrixType::Random(rows, cols), vbig(rows, cols),
67  vsmall(rows, cols);
68 
69  vbig.fill(big);
70  vsmall.fill(small);
71 
72  VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1));
73  VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm());
74  VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm());
75  VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm());
76 
77  // test with expressions as input
78  VERIFY_IS_APPROX((one * vrand).stableNorm(), vrand.norm());
79  VERIFY_IS_APPROX((one * vrand).blueNorm(), vrand.norm());
80  VERIFY_IS_APPROX((one * vrand).hypotNorm(), vrand.norm());
81  VERIFY_IS_APPROX((one * vrand + one * vrand - one * vrand).stableNorm(), vrand.norm());
82  VERIFY_IS_APPROX((one * vrand + one * vrand - one * vrand).blueNorm(), vrand.norm());
83  VERIFY_IS_APPROX((one * vrand + one * vrand - one * vrand).hypotNorm(), vrand.norm());
84 
85  RealScalar size = static_cast<RealScalar>(m.size());
86 
87  // test numext::isfinite
88  VERIFY(!(numext::isfinite)(std::numeric_limits<RealScalar>::infinity()));
89  VERIFY(!(numext::isfinite)(sqrt(-abs(big))));
90 
91  // test overflow
92  VERIFY((numext::isfinite)(sqrt(size) * abs(big)));
93  VERIFY_IS_NOT_APPROX(sqrt(copy(vbig.squaredNorm())), abs(sqrt(size) * big)); // here the default norm must fail
94  VERIFY_IS_APPROX(vbig.stableNorm(), sqrt(size) * abs(big));
95  VERIFY_IS_APPROX(vbig.blueNorm(), sqrt(size) * abs(big));
96  VERIFY_IS_APPROX(vbig.hypotNorm(), sqrt(size) * abs(big));
97 
98  // test underflow
99  VERIFY((numext::isfinite)(sqrt(size) * abs(small)));
100  VERIFY_IS_NOT_APPROX(sqrt(copy(vsmall.squaredNorm())), abs(sqrt(size) * small)); // here the default norm must fail
101  VERIFY_IS_APPROX(vsmall.stableNorm(), sqrt(size) * abs(small));
102  VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size) * abs(small));
103  VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size) * abs(small));
104 
105  // Test compilation of cwise() version
106  VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm());
107  VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm());
108  VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm());
109  VERIFY_IS_APPROX(vrand.rowwise().stableNorm(), vrand.rowwise().norm());
110  VERIFY_IS_APPROX(vrand.rowwise().blueNorm(), vrand.rowwise().norm());
111  VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(), vrand.rowwise().norm());
112 
113  // test NaN, +inf, -inf
114  MatrixType v;
115  Index i = internal::random<Index>(0, rows - 1);
116  Index j = internal::random<Index>(0, cols - 1);
117 
118  // NaN
119  {
120  v = vrand;
121  v(i, j) = std::numeric_limits<RealScalar>::quiet_NaN();
122  VERIFY(!(numext::isfinite)(v.squaredNorm()));
123  VERIFY((numext::isnan)(v.squaredNorm()));
124  VERIFY(!(numext::isfinite)(v.norm()));
125  VERIFY((numext::isnan)(v.norm()));
126  VERIFY(!(numext::isfinite)(v.stableNorm()));
127  VERIFY((numext::isnan)(v.stableNorm()));
128  VERIFY(!(numext::isfinite)(v.blueNorm()));
129  VERIFY((numext::isnan)(v.blueNorm()));
130  VERIFY(!(numext::isfinite)(v.hypotNorm()));
131  VERIFY((numext::isnan)(v.hypotNorm()));
132  }
133 
134  // +inf
135  {
136  v = vrand;
137  v(i, j) = std::numeric_limits<RealScalar>::infinity();
138  VERIFY(!(numext::isfinite)(v.squaredNorm()));
139  VERIFY(isPlusInf(v.squaredNorm()));
140  VERIFY(!(numext::isfinite)(v.norm()));
141  VERIFY(isPlusInf(v.norm()));
142  VERIFY(!(numext::isfinite)(v.stableNorm()));
143  if (complex_real_product_ok) {
144  VERIFY(isPlusInf(v.stableNorm()));
145  }
146  VERIFY(!(numext::isfinite)(v.blueNorm()));
147  VERIFY(isPlusInf(v.blueNorm()));
148  VERIFY(!(numext::isfinite)(v.hypotNorm()));
149  VERIFY(isPlusInf(v.hypotNorm()));
150  }
151 
152  // -inf
153  {
154  v = vrand;
155  v(i, j) = -std::numeric_limits<RealScalar>::infinity();
156  VERIFY(!(numext::isfinite)(v.squaredNorm()));
157  VERIFY(isPlusInf(v.squaredNorm()));
158  VERIFY(!(numext::isfinite)(v.norm()));
159  VERIFY(isPlusInf(v.norm()));
160  VERIFY(!(numext::isfinite)(v.stableNorm()));
161  if (complex_real_product_ok) {
162  VERIFY(isPlusInf(v.stableNorm()));
163  }
164  VERIFY(!(numext::isfinite)(v.blueNorm()));
165  VERIFY(isPlusInf(v.blueNorm()));
166  VERIFY(!(numext::isfinite)(v.hypotNorm()));
167  VERIFY(isPlusInf(v.hypotNorm()));
168  }
169 
170  // mix
171  {
172  Index i2 = internal::random<Index>(0, rows - 1);
173  Index j2 = internal::random<Index>(0, cols - 1);
174  v = vrand;
175  v(i, j) = -std::numeric_limits<RealScalar>::infinity();
176  v(i2, j2) = std::numeric_limits<RealScalar>::quiet_NaN();
177  VERIFY(!(numext::isfinite)(v.squaredNorm()));
178  VERIFY((numext::isnan)(v.squaredNorm()));
179  VERIFY(!(numext::isfinite)(v.norm()));
180  VERIFY((numext::isnan)(v.norm()));
181  VERIFY(!(numext::isfinite)(v.stableNorm()));
182  VERIFY((numext::isnan)(v.stableNorm()));
183  VERIFY(!(numext::isfinite)(v.blueNorm()));
184  VERIFY((numext::isnan)(v.blueNorm()));
185  if (i2 != i || j2 != j) {
186  // hypot propagates inf over NaN.
187  VERIFY(!(numext::isfinite)(v.hypotNorm()));
188  VERIFY((numext::isinf)(v.hypotNorm()));
189  } else {
190  // inf is overwritten by NaN, expect norm to be NaN.
191  VERIFY(!(numext::isfinite)(v.hypotNorm()));
192  VERIFY((numext::isnan)(v.hypotNorm()));
193  }
194  }
195 
196  // stableNormalize[d]
197  {
198  VERIFY_IS_APPROX(vrand.stableNormalized(), vrand.normalized());
199  MatrixType vcopy(vrand);
200  vcopy.stableNormalize();
201  VERIFY_IS_APPROX(vcopy, vrand.normalized());
202  VERIFY_IS_APPROX((vrand.stableNormalized()).norm(), RealScalar(1));
203  VERIFY_IS_APPROX(vcopy.norm(), RealScalar(1));
204  VERIFY_IS_APPROX((vbig.stableNormalized()).norm(), RealScalar(1));
205  VERIFY_IS_APPROX((vsmall.stableNormalized()).norm(), RealScalar(1));
206  RealScalar big_scaling = ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
207  VERIFY_IS_APPROX(vbig / big_scaling, (vbig.stableNorm() * vbig.stableNormalized()).eval() / big_scaling);
208  VERIFY_IS_APPROX(vsmall, vsmall.stableNorm() * vsmall.stableNormalized());
209  }
210 }
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:46
blueNorm
EIGEN_DONT_INLINE T::Scalar blueNorm(T &v)
Definition: bench_norm.cpp:24
hypotNorm
EIGEN_DONT_INLINE T::Scalar hypotNorm(T &v)
Definition: bench_norm.cpp:19
stableNorm
EIGEN_DONT_INLINE T::Scalar stableNorm(T &v)
Definition: bench_norm.cpp:14
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:13
VERIFY_IS_NOT_APPROX
#define VERIFY_IS_NOT_APPROX(a, b)
Definition: integer_types.cpp:15
m
int * m
Definition: level2_cplx_impl.h:294
isfinite
#define isfinite(X)
Definition: main.h:111
VERIFY
#define VERIFY(a)
Definition: main.h:362
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:371
isnan
#define isnan(X)
Definition: main.h:109
isinf
#define isinf(X)
Definition: main.h:110
Eigen::numext::abs2
EIGEN_DEVICE_FUNC bool abs2(bool x)
Definition: MathFunctions.h:1102
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:83
Eigen::isPlusInf
bool isPlusInf(const T &x)
Definition: main.h:715
constants::inf
const Mdouble inf
Definition: GeneralDefine.h:23
oomph::PseudoSolidHelper::Zero
double Zero
Definition: pseudosolid_node_update_elements.cc:35
copy
EIGEN_DONT_INLINE T copy(const T &x)
Definition: stable_norm.cpp:13
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:217
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References abs(), Eigen::numext::abs2(), blueNorm(), cols, copy(), e(), hypotNorm(), i, constants::inf, isfinite, isinf, isnan, Eigen::isPlusInf(), j, m, max, min, rows, size, sqrt(), stableNorm(), v, VERIFY, VERIFY_IS_APPROX, VERIFY_IS_MUCH_SMALLER_THAN, VERIFY_IS_NOT_APPROX, and oomph::PseudoSolidHelper::Zero.

Referenced by EIGEN_DECLARE_TEST().

◆ test_empty()

void test_empty ( )
212  {
213  Eigen::VectorXf empty(0);
214  VERIFY_IS_EQUAL(empty.stableNorm(), 0.0f);
215 }
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:367

References VERIFY_IS_EQUAL.

Referenced by EIGEN_DECLARE_TEST().

◆ test_hypot()

template<typename Scalar >
void test_hypot ( )
218  {
219  typedef typename NumTraits<Scalar>::Real RealScalar;
220  Scalar factor = internal::random<Scalar>();
221  while (numext::abs2(factor) < RealScalar(1e-4)) factor = internal::random<Scalar>();
222  Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
223 
224  factor = internal::random<Scalar>();
225  while (numext::abs2(factor) < RealScalar(1e-4)) factor = internal::random<Scalar>();
226  Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
227 
228  Scalar one(1), zero(0), sqrt2(std::sqrt(2)), nan(std::numeric_limits<RealScalar>::quiet_NaN());
229 
230  Scalar a = internal::random<Scalar>(-1, 1);
231  Scalar b = internal::random<Scalar>(-1, 1);
232  VERIFY_IS_APPROX(numext::hypot(a, b), std::sqrt(numext::abs2(a) + numext::abs2(b)));
233  VERIFY_IS_EQUAL(numext::hypot(zero, zero), zero);
234  VERIFY_IS_APPROX(numext::hypot(one, one), sqrt2);
235  VERIFY_IS_APPROX(numext::hypot(big, big), sqrt2 * numext::abs(big));
236  VERIFY_IS_APPROX(numext::hypot(small, small), sqrt2 * numext::abs(small));
237  VERIFY_IS_APPROX(numext::hypot(small, big), numext::abs(big));
238  VERIFY((numext::isnan)(numext::hypot(nan, a)));
239  VERIFY((numext::isnan)(numext::hypot(a, nan)));
240 }
b
Scalar * b
Definition: benchVecAdd.cpp:17
a
const Scalar * a
Definition: level2_cplx_impl.h:32
zero
EIGEN_DONT_INLINE Scalar zero()
Definition: svd_common.h:232

References a, abs(), Eigen::numext::abs2(), b, e(), isnan, max, min, sqrt(), VERIFY, VERIFY_IS_APPROX, VERIFY_IS_EQUAL, and zero().

Referenced by EIGEN_DECLARE_TEST().