1#include <typeinfo> 2#include <iostream> 3#include <Eigen/Core> 4#include "BenchTimer.h" 5using namespace Eigen; 6using namespace std; 7 8template<typename T> 9EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v) 10{ 11 return v.norm(); 12} 13 14template<typename T> 15EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v) 16{ 17 return v.hypotNorm(); 18} 19 20template<typename T> 21EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v) 22{ 23 return v.blueNorm(); 24} 25 26template<typename T> 27EIGEN_DONT_INLINE typename T::Scalar lapackNorm(T& v) 28{ 29 typedef typename T::Scalar Scalar; 30 int n = v.size(); 31 Scalar scale = 0; 32 Scalar ssq = 1; 33 for (int i=0;i<n;++i) 34 { 35 Scalar ax = internal::abs(v.coeff(i)); 36 if (scale >= ax) 37 { 38 ssq += internal::abs2(ax/scale); 39 } 40 else 41 { 42 ssq = Scalar(1) + ssq * internal::abs2(scale/ax); 43 scale = ax; 44 } 45 } 46 return scale * internal::sqrt(ssq); 47} 48 49template<typename T> 50EIGEN_DONT_INLINE typename T::Scalar twopassNorm(T& v) 51{ 52 typedef typename T::Scalar Scalar; 53 Scalar s = v.cwise().abs().maxCoeff(); 54 return s*(v/s).norm(); 55} 56 57template<typename T> 58EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v) 59{ 60 return v.stableNorm(); 61} 62 63template<typename T> 64EIGEN_DONT_INLINE typename T::Scalar divacNorm(T& v) 65{ 66 int n =v.size() / 2; 67 for (int i=0;i<n;++i) 68 v(i) = v(2*i)*v(2*i) + v(2*i+1)*v(2*i+1); 69 n = n/2; 70 while (n>0) 71 { 72 for (int i=0;i<n;++i) 73 v(i) = v(2*i) + v(2*i+1); 74 n = n/2; 75 } 76 return internal::sqrt(v(0)); 77} 78 79#ifdef EIGEN_VECTORIZE 80Packet4f internal::plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); } 81Packet2d internal::plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); } 82 83Packet4f internal::pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); } 84Packet2d internal::pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); } 85#endif 86 87template<typename T> 88EIGEN_DONT_INLINE typename T::Scalar pblueNorm(const T& v) 89{ 90 #ifndef EIGEN_VECTORIZE 91 return v.blueNorm(); 92 #else 93 typedef typename T::Scalar Scalar; 94 95 static int nmax = 0; 96 static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; 97 int n; 98 99 if(nmax <= 0) 100 { 101 int nbig, ibeta, it, iemin, iemax, iexp; 102 Scalar abig, eps; 103 104 nbig = std::numeric_limits<int>::max(); // largest integer 105 ibeta = std::numeric_limits<Scalar>::radix; //NumTraits<Scalar>::Base; // base for floating-point numbers 106 it = std::numeric_limits<Scalar>::digits; //NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa 107 iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent 108 iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent 109 rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number 110 111 // Check the basic machine-dependent constants. 112 if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) 113 || (it<=4 && ibeta <= 3 ) || it<2) 114 { 115 eigen_assert(false && "the algorithm cannot be guaranteed on this computer"); 116 } 117 iexp = -((1-iemin)/2); 118 b1 = std::pow(ibeta, iexp); // lower boundary of midrange 119 iexp = (iemax + 1 - it)/2; 120 b2 = std::pow(ibeta,iexp); // upper boundary of midrange 121 122 iexp = (2-iemin)/2; 123 s1m = std::pow(ibeta,iexp); // scaling factor for lower range 124 iexp = - ((iemax+it)/2); 125 s2m = std::pow(ibeta,iexp); // scaling factor for upper range 126 127 overfl = rbig*s2m; // overfow boundary for abig 128 eps = std::pow(ibeta, 1-it); 129 relerr = internal::sqrt(eps); // tolerance for neglecting asml 130 abig = 1.0/eps - 1.0; 131 if (Scalar(nbig)>abig) nmax = abig; // largest safe n 132 else nmax = nbig; 133 } 134 135 typedef typename internal::packet_traits<Scalar>::type Packet; 136 const int ps = internal::packet_traits<Scalar>::size; 137 Packet pasml = internal::pset1(Scalar(0)); 138 Packet pamed = internal::pset1(Scalar(0)); 139 Packet pabig = internal::pset1(Scalar(0)); 140 Packet ps2m = internal::pset1(s2m); 141 Packet ps1m = internal::pset1(s1m); 142 Packet pb2 = internal::pset1(b2); 143 Packet pb1 = internal::pset1(b1); 144 for(int j=0; j<v.size(); j+=ps) 145 { 146 Packet ax = internal::pabs(v.template packet<Aligned>(j)); 147 Packet ax_s2m = internal::pmul(ax,ps2m); 148 Packet ax_s1m = internal::pmul(ax,ps1m); 149 Packet maskBig = internal::plt(pb2,ax); 150 Packet maskSml = internal::plt(ax,pb1); 151 152// Packet maskMed = internal::pand(maskSml,maskBig); 153// Packet scale = internal::pset1(Scalar(0)); 154// scale = internal::por(scale, internal::pand(maskBig,ps2m)); 155// scale = internal::por(scale, internal::pand(maskSml,ps1m)); 156// scale = internal::por(scale, internal::pandnot(internal::pset1(Scalar(1)),maskMed)); 157// ax = internal::pmul(ax,scale); 158// ax = internal::pmul(ax,ax); 159// pabig = internal::padd(pabig, internal::pand(maskBig, ax)); 160// pasml = internal::padd(pasml, internal::pand(maskSml, ax)); 161// pamed = internal::padd(pamed, internal::pandnot(ax,maskMed)); 162 163 164 pabig = internal::padd(pabig, internal::pand(maskBig, internal::pmul(ax_s2m,ax_s2m))); 165 pasml = internal::padd(pasml, internal::pand(maskSml, internal::pmul(ax_s1m,ax_s1m))); 166 pamed = internal::padd(pamed, internal::pandnot(internal::pmul(ax,ax),internal::pand(maskSml,maskBig))); 167 } 168 Scalar abig = internal::predux(pabig); 169 Scalar asml = internal::predux(pasml); 170 Scalar amed = internal::predux(pamed); 171 if(abig > Scalar(0)) 172 { 173 abig = internal::sqrt(abig); 174 if(abig > overfl) 175 { 176 eigen_assert(false && "overflow"); 177 return rbig; 178 } 179 if(amed > Scalar(0)) 180 { 181 abig = abig/s2m; 182 amed = internal::sqrt(amed); 183 } 184 else 185 { 186 return abig/s2m; 187 } 188 189 } 190 else if(asml > Scalar(0)) 191 { 192 if (amed > Scalar(0)) 193 { 194 abig = internal::sqrt(amed); 195 amed = internal::sqrt(asml) / s1m; 196 } 197 else 198 { 199 return internal::sqrt(asml)/s1m; 200 } 201 } 202 else 203 { 204 return internal::sqrt(amed); 205 } 206 asml = std::min(abig, amed); 207 abig = std::max(abig, amed); 208 if(asml <= abig*relerr) 209 return abig; 210 else 211 return abig * internal::sqrt(Scalar(1) + internal::abs2(asml/abig)); 212 #endif 213} 214 215#define BENCH_PERF(NRM) { \ 216 Eigen::BenchTimer tf, td, tcf; tf.reset(); td.reset(); tcf.reset();\ 217 for (int k=0; k<tries; ++k) { \ 218 tf.start(); \ 219 for (int i=0; i<iters; ++i) NRM(vf); \ 220 tf.stop(); \ 221 } \ 222 for (int k=0; k<tries; ++k) { \ 223 td.start(); \ 224 for (int i=0; i<iters; ++i) NRM(vd); \ 225 td.stop(); \ 226 } \ 227 for (int k=0; k<std::max(1,tries/3); ++k) { \ 228 tcf.start(); \ 229 for (int i=0; i<iters; ++i) NRM(vcf); \ 230 tcf.stop(); \ 231 } \ 232 std::cout << #NRM << "\t" << tf.value() << " " << td.value() << " " << tcf.value() << "\n"; \ 233} 234 235void check_accuracy(double basef, double based, int s) 236{ 237 double yf = basef * internal::abs(internal::random<double>()); 238 double yd = based * internal::abs(internal::random<double>()); 239 VectorXf vf = VectorXf::Ones(s) * yf; 240 VectorXd vd = VectorXd::Ones(s) * yd; 241 242 std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n"; 243 std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\n"; 244 std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\n"; 245 std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\n"; 246 std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\n"; 247 std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\n"; 248 std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\n"; 249 std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\n"; 250} 251 252void check_accuracy_var(int ef0, int ef1, int ed0, int ed1, int s) 253{ 254 VectorXf vf(s); 255 VectorXd vd(s); 256 for (int i=0; i<s; ++i) 257 { 258 vf[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ef0,ef1)); 259 vd[i] = internal::abs(internal::random<double>()) * std::pow(double(10), internal::random<int>(ed0,ed1)); 260 } 261 262 //std::cout << "reference\t" << internal::sqrt(double(s))*yf << "\t" << internal::sqrt(double(s))*yd << "\n"; 263 std::cout << "sqsumNorm\t" << sqsumNorm(vf) << "\t" << sqsumNorm(vd) << "\t" << sqsumNorm(vf.cast<long double>()) << "\t" << sqsumNorm(vd.cast<long double>()) << "\n"; 264 std::cout << "hypotNorm\t" << hypotNorm(vf) << "\t" << hypotNorm(vd) << "\t" << hypotNorm(vf.cast<long double>()) << "\t" << hypotNorm(vd.cast<long double>()) << "\n"; 265 std::cout << "blueNorm\t" << blueNorm(vf) << "\t" << blueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n"; 266 std::cout << "pblueNorm\t" << pblueNorm(vf) << "\t" << pblueNorm(vd) << "\t" << blueNorm(vf.cast<long double>()) << "\t" << blueNorm(vd.cast<long double>()) << "\n"; 267 std::cout << "lapackNorm\t" << lapackNorm(vf) << "\t" << lapackNorm(vd) << "\t" << lapackNorm(vf.cast<long double>()) << "\t" << lapackNorm(vd.cast<long double>()) << "\n"; 268 std::cout << "twopassNorm\t" << twopassNorm(vf) << "\t" << twopassNorm(vd) << "\t" << twopassNorm(vf.cast<long double>()) << "\t" << twopassNorm(vd.cast<long double>()) << "\n"; 269// std::cout << "bl2passNorm\t" << bl2passNorm(vf) << "\t" << bl2passNorm(vd) << "\t" << bl2passNorm(vf.cast<long double>()) << "\t" << bl2passNorm(vd.cast<long double>()) << "\n"; 270} 271 272int main(int argc, char** argv) 273{ 274 int tries = 10; 275 int iters = 100000; 276 double y = 1.1345743233455785456788e12 * internal::random<double>(); 277 VectorXf v = VectorXf::Ones(1024) * y; 278 279// return 0; 280 int s = 10000; 281 double basef_ok = 1.1345743233455785456788e15; 282 double based_ok = 1.1345743233455785456788e95; 283 284 double basef_under = 1.1345743233455785456788e-27; 285 double based_under = 1.1345743233455785456788e-303; 286 287 double basef_over = 1.1345743233455785456788e+27; 288 double based_over = 1.1345743233455785456788e+302; 289 290 std::cout.precision(20); 291 292 std::cerr << "\nNo under/overflow:\n"; 293 check_accuracy(basef_ok, based_ok, s); 294 295 std::cerr << "\nUnderflow:\n"; 296 check_accuracy(basef_under, based_under, s); 297 298 std::cerr << "\nOverflow:\n"; 299 check_accuracy(basef_over, based_over, s); 300 301 std::cerr << "\nVarying (over):\n"; 302 for (int k=0; k<1; ++k) 303 { 304 check_accuracy_var(20,27,190,302,s); 305 std::cout << "\n"; 306 } 307 308 std::cerr << "\nVarying (under):\n"; 309 for (int k=0; k<1; ++k) 310 { 311 check_accuracy_var(-27,20,-302,-190,s); 312 std::cout << "\n"; 313 } 314 315 std::cout.precision(4); 316 std::cerr << "Performance (out of cache):\n"; 317 { 318 int iters = 1; 319 VectorXf vf = VectorXf::Random(1024*1024*32) * y; 320 VectorXd vd = VectorXd::Random(1024*1024*32) * y; 321 VectorXcf vcf = VectorXcf::Random(1024*1024*32) * y; 322 BENCH_PERF(sqsumNorm); 323 BENCH_PERF(blueNorm); 324// BENCH_PERF(pblueNorm); 325// BENCH_PERF(lapackNorm); 326// BENCH_PERF(hypotNorm); 327// BENCH_PERF(twopassNorm); 328 BENCH_PERF(bl2passNorm); 329 } 330 331 std::cerr << "\nPerformance (in cache):\n"; 332 { 333 int iters = 100000; 334 VectorXf vf = VectorXf::Random(512) * y; 335 VectorXd vd = VectorXd::Random(512) * y; 336 VectorXcf vcf = VectorXcf::Random(512) * y; 337 BENCH_PERF(sqsumNorm); 338 BENCH_PERF(blueNorm); 339// BENCH_PERF(pblueNorm); 340// BENCH_PERF(lapackNorm); 341// BENCH_PERF(hypotNorm); 342// BENCH_PERF(twopassNorm); 343 BENCH_PERF(bl2passNorm); 344 } 345} 346