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>
sqsumNorm(const T & v)9 EIGEN_DONT_INLINE typename T::Scalar sqsumNorm(const T& v)
10 {
11 return v.norm();
12 }
13
14 template<typename T>
hypotNorm(const T & v)15 EIGEN_DONT_INLINE typename T::Scalar hypotNorm(const T& v)
16 {
17 return v.hypotNorm();
18 }
19
20 template<typename T>
blueNorm(const T & v)21 EIGEN_DONT_INLINE typename T::Scalar blueNorm(const T& v)
22 {
23 return v.blueNorm();
24 }
25
26 template<typename T>
lapackNorm(T & v)27 EIGEN_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
49 template<typename T>
twopassNorm(T & v)50 EIGEN_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
57 template<typename T>
bl2passNorm(T & v)58 EIGEN_DONT_INLINE typename T::Scalar bl2passNorm(T& v)
59 {
60 return v.stableNorm();
61 }
62
63 template<typename T>
divacNorm(T & v)64 EIGEN_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
plt(const Packet4f & a,Packet4f & b)80 Packet4f internal::plt(const Packet4f& a, Packet4f& b) { return _mm_cmplt_ps(a,b); }
plt(const Packet2d & a,Packet2d & b)81 Packet2d internal::plt(const Packet2d& a, Packet2d& b) { return _mm_cmplt_pd(a,b); }
82
pandnot(const Packet4f & a,Packet4f & b)83 Packet4f internal::pandnot(const Packet4f& a, Packet4f& b) { return _mm_andnot_ps(a,b); }
pandnot(const Packet2d & a,Packet2d & b)84 Packet2d internal::pandnot(const Packet2d& a, Packet2d& b) { return _mm_andnot_pd(a,b); }
85 #endif
86
87 template<typename T>
pblueNorm(const T & v)88 EIGEN_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
check_accuracy(double basef,double based,int s)235 void 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
check_accuracy_var(int ef0,int ef1,int ed0,int ed1,int s)252 void 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
main(int argc,char ** argv)272 int 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