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