mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

bugfix in blueNorm

Browse Source
This commit is contained in:
Gael Guennebaud 2009-07-16 14:20:36 +02:00
parent 65fc70b750
commit 525da6a464
3 changed files with 89 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
Eigen/src/Core/Dot.h
View File

@ -304,29 +304,6 @@ MatrixBase<Derived>::stableNorm() const
return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>());
}
/** \internal Computes ibeta^iexp by binary expansion of iexp,
* exact if ibeta is the machine base */
template<typename T> inline T bexp(int ibeta, int iexp)
{
T tbeta = T(ibeta);
T res = 1.0;
int n = iexp;
if (n<0)
{
n = - n;
tbeta = 1.0/tbeta;
}
for(;;)
{
if ((n % 2)==0)
res = res * tbeta;
n = n/2;
if (n==0) return res;
tbeta = tbeta*tbeta;
}
return res;
}
/** \returns the \em l2 norm of \c *this using the Blue's algorithm.
* A Portable Fortran Program to Find the Euclidean Norm of a Vector,
* ACM TOMS, Vol 4, Issue 1, 1978.
@ -337,7 +314,7 @@ template<typename Derived>
inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
MatrixBase<Derived>::blueNorm() const
{
static int nmax;
static int nmax = -1;
static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
int n;
Scalar ax, abig, amed, asml;
@ -355,8 +332,8 @@ MatrixBase<Derived>::blueNorm() const
// are used. For any specific computer, each of the assignment
// statements can be replaced
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers
it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers
it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
@ -368,17 +345,17 @@ MatrixBase<Derived>::blueNorm() const
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange
b1 = std::pow(ibeta, iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange
b2 = std::pow(ibeta,iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range
s1m = std::pow(ibeta,iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range
s2m = std::pow(ibeta,iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = bexp<Scalar>(ibeta, 1-it);
eps = std::pow(ibeta, 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (Scalar(nbig)>abig) nmax = abig; // largest safe n

8
Eigen/src/Core/NumTraits.h
View File

@ -70,9 +70,7 @@ template<> struct NumTraits<float>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1,
Base = 2,
Mantissa = 23
MulCost = 1
};
};
@ -85,9 +83,7 @@ template<> struct NumTraits<double>
HasFloatingPoint = 1,
ReadCost = 1,
AddCost = 1,
MulCost = 1,
Base = 2,
Mantissa = 52
MulCost = 1
};
};

126
bench/bench_norm.cpp
View File

@ -1,3 +1,4 @@
#include <typeinfo>
#include <Eigen/Core>
#include "BenchTimer.h"
using namespace Eigen;
@ -58,18 +59,23 @@ EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v)
return ei_sqrt(v(0));
}
#ifdef EIGEN_VECTORIZE
Packet4f ei_plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); }
Packet2d ei_plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); }
Packet4f ei_pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); }
Packet2d ei_pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); }
#endif
template<typename T>
EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
{
#ifndef EIGEN_VECTORIZE
return v.blueNorm();
#else
typedef typename T::Scalar Scalar;
static int nmax;
static int nmax = 0;
static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr;
int n;
@ -79,8 +85,8 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
Scalar abig, eps;
nbig = std::numeric_limits<int>::max(); // largest integer
ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers
it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers
it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa
iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent
iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent
rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number
@ -92,23 +98,23 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
ei_assert(false && "the algorithm cannot be guaranteed on this computer");
}
iexp = -((1-iemin)/2);
b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange
b1 = std::pow(ibeta, iexp); // lower boundary of midrange
iexp = (iemax + 1 - it)/2;
b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange
b2 = std::pow(ibeta,iexp); // upper boundary of midrange
iexp = (2-iemin)/2;
s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range
s1m = std::pow(ibeta,iexp); // scaling factor for lower range
iexp = - ((iemax+it)/2);
s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range
s2m = std::pow(ibeta,iexp); // scaling factor for upper range
overfl = rbig*s2m; // overfow boundary for abig
eps = bexp<Scalar>(ibeta, 1-it);
eps = std::pow(ibeta, 1-it);
relerr = ei_sqrt(eps); // tolerance for neglecting asml
abig = 1.0/eps - 1.0;
if (Scalar(nbig)>abig) nmax = abig; // largest safe n
else nmax = nbig;
}
typedef typename ei_packet_traits<Scalar>::type Packet;
const int ps = ei_packet_traits<Scalar>::size;
Packet pasml = ei_pset1(Scalar(0));
@ -173,6 +179,7 @@ EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v)
return abig;
else
return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig));
#endif
}
#define BENCH_PERF(NRM) { \
@ -196,7 +203,7 @@ void check_accuracy(double basef, double based, int s)
double yd = based * ei_abs(ei_random<double>());
VectorXf vf = VectorXf::Ones(s) * yf;
VectorXd vd = VectorXd::Ones(s) * yd;
std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_sqrt(double(s))*yd << "\n";
std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n";
std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n";
@ -205,55 +212,80 @@ void check_accuracy(double basef, double based, int s)
std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n";
}
int main(int argc, char** argv)
void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s)
{
VectorXf vf(s);
VectorXd vd(s);
for (int i=0; i<s; ++i)
{
vf[i] = ei_abs(ei_random<double>()) * std::pow(double(10), ei_random<int>(ef0,ef1));
vd[i] = ei_abs(ei_random<double>()) * std::pow(double(10), ei_random<int>(ed0,ed1));
}
//std::cout << "reference\t" << ei_sqrt(double(s))*yf << "\t" << ei_sqrt(double(s))*yd << "\n";
std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast<long double>()) << "\t" << sqsumNorm(vd.cast<long double>()) << "\n";
std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast<long double>()) << "\t" << hypotNorm(vd.cast<long double>()) << "\n";
std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n";
std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>()) << "\t" << lapackNorm(vd.cast<long double>()) << "\n";
}
int main(int argc, char** argv)
{
int tries = 5;
int iters = 100000;
double y = 1.1345743233455785456788e12 * ei_random<double>();
VectorXf v = VectorXf::Ones(1024) * y;
// std::cerr << "Performance (out of cache):\n";
// {
// int iters = 1;
// VectorXf vf = VectorXf::Ones(1024*1024*32) * y;
// VectorXd vd = VectorXd::Ones(1024*1024*32) * y;
// BENCH_PERF(sqsumNorm);
// BENCH_PERF(blueNorm);
// BENCH_PERF(pblueNorm);
// BENCH_PERF(lapackNorm);
// BENCH_PERF(hypotNorm);
// }
//
// std::cerr << "\nPerformance (in cache):\n";
// {
// int iters = 100000;
// VectorXf vf = VectorXf::Ones(512) * y;
// VectorXd vd = VectorXd::Ones(512) * y;
// BENCH_PERF(sqsumNorm);
// BENCH_PERF(blueNorm);
// BENCH_PERF(pblueNorm);
// BENCH_PERF(lapackNorm);
// BENCH_PERF(hypotNorm);
// }
std::cerr << "Performance (out of cache):\n";
{
int iters = 1;
VectorXf vf = VectorXf::Ones(1024*1024*32) * y;
VectorXd vd = VectorXd::Ones(1024*1024*32) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(blueNorm);
BENCH_PERF(pblueNorm);
BENCH_PERF(lapackNorm);
BENCH_PERF(hypotNorm);
}
std::cerr << "\nPerformance (in cache):\n";
{
int iters = 100000;
VectorXf vf = VectorXf::Ones(512) * y;
VectorXd vd = VectorXd::Ones(512) * y;
BENCH_PERF(sqsumNorm);
BENCH_PERF(blueNorm);
BENCH_PERF(pblueNorm);
BENCH_PERF(lapackNorm);
BENCH_PERF(hypotNorm);
}
int s = 10000;
double basef_ok = 1.1345743233455785456788e12;
double based_ok = 1.1345743233455785456788e32;
double basef_under = 1.1345743233455785456788e-23;
double based_under = 1.1345743233455785456788e-180;
double basef_ok = 1.1345743233455785456788e15;
double based_ok = 1.1345743233455785456788e95;
double basef_under = 1.1345743233455785456788e-27;
double based_under = 1.1345743233455785456788e-315;
double basef_over = 1.1345743233455785456788e+27;
double based_over = 1.1345743233455785456788e+185;
double based_over = 1.1345743233455785456788e+302;
std::cout.precision(20);
std::cerr << "\nNo under/overflow:\n";
check_accuracy(basef_ok, based_ok, s);
std::cerr << "\nUnderflow:\n";
check_accuracy(basef_under, based_under, 1);
std::cerr << "\nOverflow:\n";
check_accuracy(basef_over, based_over, s);
std::cerr << "\nVarying (over):\n";
for (int k=0; k<5; ++k)
{
check_accuracy_var(20,27,190,302,s);
std::cout << "\n";
}
}