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