10 #ifndef EIGEN_STABLENORM_H
11 #define EIGEN_STABLENORM_H
16 template<
typename ExpressionType,
typename Scalar>
17 inline void stable_norm_kernel(
const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
19 Scalar max = bl.cwiseAbs().maxCoeff();
22 ssq = ssq * abs2(scale/max);
24 invScale = Scalar(1)/scale;
28 ssq += (bl*invScale).squaredNorm();
42 template<
typename Derived>
43 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
47 const Index blockSize = 4096;
49 RealScalar invScale(1);
55 Index bi = internal::first_aligned(derived());
57 internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
58 for (; bi<n; bi+=blockSize)
59 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).
template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
60 return scale * internal::sqrt(ssq);
72 template<
typename Derived>
79 static Index nmax = -1;
80 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
83 int nbig, ibeta, it, iemin, iemax, iexp;
93 nbig = (std::numeric_limits<Index>::max)();
94 ibeta = std::numeric_limits<RealScalar>::radix;
95 it = std::numeric_limits<RealScalar>::digits;
96 iemin = std::numeric_limits<RealScalar>::min_exponent;
97 iemax = std::numeric_limits<RealScalar>::max_exponent;
98 rbig = (std::numeric_limits<RealScalar>::max)();
100 iexp = -((1-iemin)/2);
101 b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
102 iexp = (iemax + 1 - it)/2;
103 b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
106 s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
107 iexp = - ((iemax+it)/2);
108 s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
111 eps = RealScalar(pow(
double(ibeta), 1-it));
112 relerr = internal::sqrt(eps);
113 abig = RealScalar(1.0/eps - 1.0);
114 if (RealScalar(nbig)>abig) nmax = int(abig);
118 RealScalar ab2 = b2 / RealScalar(n);
119 RealScalar asml = RealScalar(0);
120 RealScalar amed = RealScalar(0);
121 RealScalar abig = RealScalar(0);
122 for(
Index j=0; j<n; ++j)
124 RealScalar ax = internal::abs(coeff(j));
125 if(ax > ab2) abig += internal::abs2(ax*s2m);
126 else if(ax < b1) asml += internal::abs2(ax*s1m);
127 else amed += internal::abs2(ax);
129 if(abig > RealScalar(0))
131 abig = internal::sqrt(abig);
136 if(amed > RealScalar(0))
139 amed = internal::sqrt(amed);
144 else if(asml > RealScalar(0))
146 if (amed > RealScalar(0))
148 abig = internal::sqrt(amed);
149 amed = internal::sqrt(asml) / s1m;
152 return internal::sqrt(asml)/s1m;
155 return internal::sqrt(amed);
156 asml = (min)(abig, amed);
157 abig = (max)(abig, amed);
158 if(asml <= abig*relerr)
161 return abig * internal::sqrt(RealScalar(1) + internal::abs2(asml/abig));
169 template<
typename Derived>
173 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
178 #endif // EIGEN_STABLENORM_H