Namespaces | Macros | Functions
bench_norm.cpp File Reference
#include <typeinfo>
#include <iostream>
#include <Eigen/Core>
#include "BenchTimer.h"

Namespaces

 Eigen
 Namespace containing all symbols from the Eigen library.
 
 Eigen::internal
 Namespace containing low-level routines from the Eigen library.
 

Macros

#define BENCH_PERF(NRM)
 

Functions

template<typename T >
EIGEN_DONT_INLINE T::Scalar sqsumNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar stableNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar hypotNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar blueNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar lapackNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar twopassNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar bl2passNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar divacNorm (T &v)
 
template<typename T >
EIGEN_DONT_INLINE T::Scalar pblueNorm (const T &v)
 
void check_accuracy (double basef, double based, int s)
 
void check_accuracy_var (int ef0, int ef1, int ed0, int ed1, int s)
 
int main (int argc, char **argv)
 

Macro Definition Documentation

◆ BENCH_PERF

#define BENCH_PERF (   NRM)
Value:
{ \
float af = 0; \
double ad = 0; \
std::complex<float> ac = 0; \
Eigen::BenchTimer tf, td, tcf; \
tf.reset(); \
td.reset(); \
tcf.reset(); \
for (int k = 0; k < tries; ++k) { \
tf.start(); \
for (int i = 0; i < iters; ++i) { \
af += NRM(vf); \
} \
tf.stop(); \
} \
for (int k = 0; k < tries; ++k) { \
td.start(); \
for (int i = 0; i < iters; ++i) { \
ad += NRM(vd); \
} \
td.stop(); \
} \
/*for (int k=0; k<std::max(1,tries/3); ++k) { \
tcf.start(); \
for (int i=0; i<iters; ++i) { ac += NRM(vcf); } \
tcf.stop(); \
} */ \
std::cout << #NRM << "\t" << tf.value() << " " << td.value() << " " << tcf.value() << "\n"; \
}
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
k
char char char int int * k
Definition: level2_impl.h:374

Function Documentation

◆ bl2passNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar bl2passNorm ( T &  v)
54  {
55  return v.stableNorm();
56 }
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1

References v.

◆ blueNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar blueNorm ( T &  v)
24  {
25  return v.blueNorm();
26 }

References v.

Referenced by check_accuracy(), and stable_norm().

◆ check_accuracy()

void check_accuracy ( double  basef,
double  based,
int  s 
)
227  {
228  double yf = basef * std::abs(internal::random<double>());
229  double yd = based * std::abs(internal::random<double>());
230  VectorXf vf = VectorXf::Ones(s) * yf;
231  VectorXd vd = VectorXd::Ones(s) * yd;
232 
233  std::cout << "reference\t" << std::sqrt(double(s)) * yf << "\t" << std::sqrt(double(s)) * yd << "\n";
234  std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n";
235  std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n";
236  std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\n";
237  std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\n";
abs
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
sqrt
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
sqsumNorm
EIGEN_DONT_INLINE T::Scalar sqsumNorm(T &v)
Definition: bench_norm.cpp:9
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
pblueNorm
EIGEN_DONT_INLINE T::Scalar pblueNorm(const T &v)
Definition: bench_norm.cpp:83
s
RealScalar s
Definition: level1_cplx_impl.h:130

References abs(), blueNorm(), hypotNorm(), lapackNorm(), pblueNorm(), s, sqrt(), and sqsumNorm().

◆ check_accuracy_var()

void check_accuracy_var ( int  ef0,
int  ef1,
int  ed0,
int  ed1,
int  s 
)
243  {
244  VectorXf vf(s);
245  VectorXd vd(s);
246  for (int i = 0; i < s; ++i) {
247  vf[i] = std::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ef0, ef1));
248  vd[i] = std::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ed0, ed1));
249  }
250 
251  // std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n";
252  std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast<long double>())
253  << "\t" << sqsumNorm(vd.cast<long double>()) << "\n";
254  std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast<long double>())
255  << "\t" << hypotNorm(vd.cast<long double>()) << "\n";
256  std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t"
257  << blueNorm(vd.cast<long double>()) << "\n";
258  std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast<long double>())
259  << "\t" << blueNorm(vd.cast<long double>()) << "\n";
260  std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>())
261  << "\t" << lapackNorm(vd.cast<long double>()) << "\n";
262  std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t"
twopassNorm
EIGEN_DONT_INLINE T::Scalar twopassNorm(T &v)
Definition: bench_norm.cpp:47
lapackNorm
EIGEN_DONT_INLINE T::Scalar lapackNorm(T &v)
Definition: bench_norm.cpp:29
Eigen::bfloat16_impl::pow
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 pow(const bfloat16 &a, const bfloat16 &b)
Definition: BFloat16.h:625

◆ divacNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar divacNorm ( T &  v)
59  {
60  int n = v.size() / 2;
61  for (int i = 0; i < n; ++i) v(i) = v(2 * i) * v(2 * i) + v(2 * i + 1) * v(2 * i + 1);
62  n = n / 2;
63  while (n > 0) {
64  for (int i = 0; i < n; ++i) v(i) = v(2 * i) + v(2 * i + 1);
65  n = n / 2;
66  }
67  return std::sqrt(v(0));
68 }
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11

References i, n, sqrt(), and v.

◆ hypotNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar hypotNorm ( T &  v)
19  {
20  return v.hypotNorm();
21 }

References v.

Referenced by check_accuracy(), and stable_norm().

◆ lapackNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar lapackNorm ( T &  v)
29  {
30  typedef typename T::Scalar Scalar;
31  int n = v.size();
32  Scalar scale = 0;
33  Scalar ssq = 1;
34  for (int i = 0; i < n; ++i) {
35  Scalar ax = std::abs(v.coeff(i));
36  if (scale >= ax) {
37  ssq += numext::abs2(ax / scale);
38  } else {
39  ssq = Scalar(1) + ssq * numext::abs2(scale / ax);
40  scale = ax;
41  }
42  }
43  return scale * std::sqrt(ssq);
44 }
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
Eigen::numext::abs2
EIGEN_DEVICE_FUNC bool abs2(bool x)
Definition: MathFunctions.h:1102
plotDoE.ax
ax
Definition: plotDoE.py:39

References abs(), Eigen::numext::abs2(), plotDoE::ax, i, n, sqrt(), and v.

Referenced by check_accuracy().

◆ main()

int main ( int argc  ,
char **  argv 
)
268  {
269  int tries = 10;
270  int iters = 100000;
271  double y = 1.1345743233455785456788e12 * internal::random<double>();
272  VectorXf v = VectorXf::Ones(1024) * y;
273 
274  // return 0;
275  int s = 10000;
276  double basef_ok = 1.1345743233455785456788e15;
277  double based_ok = 1.1345743233455785456788e95;
278 
279  double basef_under = 1.1345743233455785456788e-27;
280  double based_under = 1.1345743233455785456788e-303;
281 
282  double basef_over = 1.1345743233455785456788e+27;
283  double based_over = 1.1345743233455785456788e+302;
284 
285  std::cout.precision(20);
286 
287  std::cerr << "\nNo under/overflow:\n";
288  check_accuracy(basef_ok, based_ok, s);
289 
290  std::cerr << "\nUnderflow:\n";
291  check_accuracy(basef_under, based_under, s);
292 
293  std::cerr << "\nOverflow:\n";
294  check_accuracy(basef_over, based_over, s);
295 
296  std::cerr << "\nVarying (over):\n";
297  for (int k = 0; k < 1; ++k) {
298  check_accuracy_var(20, 27, 190, 302, s);
299  std::cout << "\n";
300  }
301 
302  std::cerr << "\nVarying (under):\n";
303  for (int k = 0; k < 1; ++k) {
304  check_accuracy_var(-27, 20, -302, -190, s);
305  std::cout << "\n";
306  }
307 
308  y = 1;
309  std::cout.precision(4);
310  int s1 = 1024 * 1024 * 32;
311  std::cerr << "Performance (out of cache, " << s1 << "):\n";
312  {
313  int iters = 1;
314  VectorXf vf = VectorXf::Random(s1) * y;
315  VectorXd vd = VectorXd::Random(s1) * y;
316  VectorXcf vcf = VectorXcf::Random(s1) * y;
317  BENCH_PERF(sqsumNorm);
318  BENCH_PERF(stableNorm);
319  BENCH_PERF(blueNorm);
320  BENCH_PERF(pblueNorm);
321  BENCH_PERF(lapackNorm);
322  BENCH_PERF(hypotNorm);
323  BENCH_PERF(twopassNorm);
324  BENCH_PERF(bl2passNorm);
325  }
326 
327  std::cerr << "\nPerformance (in cache, " << 512 << "):\n";
328  {
329  int iters = 100000;
330  VectorXf vf = VectorXf::Random(512) * y;
331  VectorXd vd = VectorXd::Random(512) * y;
332  VectorXcf vcf = VectorXcf::Random(512) * y;
333  BENCH_PERF(sqsumNorm);
334  BENCH_PERF(stableNorm);
335  BENCH_PERF(blueNorm);
336  BENCH_PERF(pblueNorm);
337  BENCH_PERF(lapackNorm);
338  BENCH_PERF(hypotNorm);
BENCH_PERF
#define BENCH_PERF(NRM)
Definition: bench_norm.cpp:196
check_accuracy
void check_accuracy(double basef, double based, int s)
Definition: bench_norm.cpp:223
check_accuracy_var
void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s)
Definition: bench_norm.cpp:239
bl2passNorm
EIGEN_DONT_INLINE T::Scalar bl2passNorm(T &v)
Definition: bench_norm.cpp:54
stableNorm
EIGEN_DONT_INLINE T::Scalar stableNorm(T &v)
Definition: bench_norm.cpp:14
y
Scalar * y
Definition: level1_cplx_impl.h:128

◆ pblueNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar pblueNorm ( const T &  v)
83  {
84 #ifndef EIGEN_VECTORIZE
85  return v.blueNorm();
86 #else
87  typedef typename T::Scalar Scalar;
88 
89  static int nmax = 0;
90  static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
91  int n;
92 
93  if (nmax <= 0) {
94  int nbig, ibeta, it, iemin, iemax, iexp;
95  Scalar abig, eps;
96 
97  nbig = NumTraits<int>::highest(); // largest integer
98  ibeta = std::numeric_limits<Scalar>::radix; // NumTraits<Scalar>::Base; // base for
99  // floating-point numbers
100  it = NumTraits<Scalar>::digits(); // NumTraits<Scalar>::Mantissa; // number of base-beta digits in
101  // mantissa
102  iemin = NumTraits<Scalar>::min_exponent(); // minimum exponent
103  iemax = NumTraits<Scalar>::max_exponent(); // maximum exponent
104  rbig = NumTraits<Scalar>::highest(); // largest floating-point number
105 
106  // Check the basic machine-dependent constants.
107  if (iemin > 1 - 2 * it || 1 + it > iemax || (it == 2 && ibeta < 5) || (it <= 4 && ibeta <= 3) || it < 2) {
108  eigen_assert(false && "the algorithm cannot be guaranteed on this computer");
109  }
110  iexp = -((1 - iemin) / 2);
111  b1 = std::pow(ibeta, iexp); // lower boundary of midrange
112  iexp = (iemax + 1 - it) / 2;
113  b2 = std::pow(ibeta, iexp); // upper boundary of midrange
114 
115  iexp = (2 - iemin) / 2;
116  s1m = std::pow(ibeta, iexp); // scaling factor for lower range
117  iexp = -((iemax + it) / 2);
118  s2m = std::pow(ibeta, iexp); // scaling factor for upper range
119 
120  overfl = rbig * s2m; // overflow boundary for abig
121  eps = std::pow(ibeta, 1 - it);
122  relerr = std::sqrt(eps); // tolerance for neglecting asml
123  abig = 1.0 / eps - 1.0;
124  if (Scalar(nbig) > abig)
125  nmax = abig; // largest safe n
126  else
127  nmax = nbig;
128  }
129 
130  typedef typename internal::packet_traits<Scalar>::type Packet;
131  const int ps = internal::packet_traits<Scalar>::size;
132  Packet pasml = internal::pset1<Packet>(Scalar(0));
133  Packet pamed = internal::pset1<Packet>(Scalar(0));
134  Packet pabig = internal::pset1<Packet>(Scalar(0));
135  Packet ps2m = internal::pset1<Packet>(s2m);
136  Packet ps1m = internal::pset1<Packet>(s1m);
137  Packet pb2 = internal::pset1<Packet>(b2);
138  Packet pb1 = internal::pset1<Packet>(b1);
139  for (int j = 0; j < v.size(); j += ps) {
140  Packet ax = internal::pabs(v.template packet<Aligned>(j));
141  Packet ax_s2m = internal::pmul(ax, ps2m);
142  Packet ax_s1m = internal::pmul(ax, ps1m);
143  Packet maskBig = internal::plt(pb2, ax);
144  Packet maskSml = internal::plt(ax, pb1);
145 
146  // Packet maskMed = internal::pand(maskSml,maskBig);
147  // Packet scale = internal::pset1(Scalar(0));
148  // scale = internal::por(scale, internal::pand(maskBig,ps2m));
149  // scale = internal::por(scale, internal::pand(maskSml,ps1m));
150  // scale = internal::por(scale, internal::pandnot(internal::pset1(Scalar(1)),maskMed));
151  // ax = internal::pmul(ax,scale);
152  // ax = internal::pmul(ax,ax);
153  // pabig = internal::padd(pabig, internal::pand(maskBig, ax));
154  // pasml = internal::padd(pasml, internal::pand(maskSml, ax));
155  // pamed = internal::padd(pamed, internal::pandnot(ax,maskMed));
156 
157  pabig = internal::padd(pabig, internal::pand(maskBig, internal::pmul(ax_s2m, ax_s2m)));
158  pasml = internal::padd(pasml, internal::pand(maskSml, internal::pmul(ax_s1m, ax_s1m)));
159  pamed = internal::padd(pamed, internal::pandnot(internal::pmul(ax, ax), internal::pand(maskSml, maskBig)));
160  }
161  Scalar abig = internal::predux(pabig);
162  Scalar asml = internal::predux(pasml);
163  Scalar amed = internal::predux(pamed);
164  if (abig > Scalar(0)) {
165  abig = std::sqrt(abig);
166  if (abig > overfl) {
167  eigen_assert(false && "overflow");
168  return rbig;
169  }
170  if (amed > Scalar(0)) {
171  abig = abig / s2m;
172  amed = std::sqrt(amed);
173  } else {
174  return abig / s2m;
175  }
176 
177  } else if (asml > Scalar(0)) {
178  if (amed > Scalar(0)) {
179  abig = std::sqrt(amed);
180  amed = std::sqrt(asml) / s1m;
181  } else {
182  return std::sqrt(asml) / s1m;
183  }
184  } else {
185  return std::sqrt(amed);
186  }
187  asml = std::min(abig, amed);
188  abig = std::max(abig, amed);
189  if (asml <= abig * relerr)
190  return abig;
191  else
192  return abig * std::sqrt(Scalar(1) + numext::abs2(asml / abig));
193 #endif
194 }
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:910
padd
EIGEN_ALWAYS_INLINE Packet2cf padd(Packet2cf &a, std::complex< float > &b)
Definition: MatrixVectorProduct.h:1277
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Packet
internal::packet_traits< Scalar >::type Packet
Definition: benchmark-blocking-sizes.cpp:54
min
#define min(a, b)
Definition: datatypes.h:22
max
#define max(a, b)
Definition: datatypes.h:23
ps
int RealScalar int RealScalar int RealScalar RealScalar * ps
Definition: level1_cplx_impl.h:124
relerr
Derived::RealScalar relerr(const MatrixBase< Derived > &A, const MatrixBase< OtherDerived > &B)
Definition: matrix_functions.h:54
CRBond_Bessel::eps
double eps
Definition: crbond_bessel.cc:24
Eigen::internal::pmul
EIGEN_STRONG_INLINE Packet4cf pmul(const Packet4cf &a, const Packet4cf &b)
Definition: AVX/Complex.h:88
Eigen::internal::pandnot
EIGEN_STRONG_INLINE Packet8h pandnot(const Packet8h &a, const Packet8h &b)
Definition: AVX/PacketMath.h:2323
Eigen::internal::pabs
EIGEN_STRONG_INLINE Packet4f pabs(const Packet4f &a)
Definition: AltiVec/PacketMath.h:1936
Eigen::internal::predux
EIGEN_DEVICE_FUNC unpacket_traits< Packet >::type predux(const Packet &a)
Definition: GenericPacketMath.h:1232
Eigen::internal::pand
EIGEN_STRONG_INLINE Packet8h pand(const Packet8h &a, const Packet8h &b)
Definition: AVX/PacketMath.h:2319
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 Eigen::numext::abs2(), plotDoE::ax, eigen_assert, CRBond_Bessel::eps, j, max, min, n, Eigen::internal::pabs(), Eigen::internal::padd(), Eigen::internal::pand(), Eigen::internal::pandnot(), Eigen::internal::pmul(), Eigen::bfloat16_impl::pow(), Eigen::internal::predux(), ps, relerr(), Eigen::internal::packet_traits< Scalar >::size, sqrt(), and v.

Referenced by check_accuracy().

◆ sqsumNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar sqsumNorm ( T &  v)
9  {
10  return v.norm();
11 }

References v.

Referenced by check_accuracy().

◆ stableNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar stableNorm ( T &  v)
14  {
15  return v.stableNorm();
16 }

References v.

Referenced by Eigen::internal::idrstabl(), and stable_norm().

◆ twopassNorm()

template<typename T >
EIGEN_DONT_INLINE T::Scalar twopassNorm ( T &  v)
47  {
48  typedef typename T::Scalar Scalar;
49  Scalar s = v.array().abs().maxCoeff();
50  return s * (v / s).norm();
51 }

References s, and v.