10 #ifndef EIGEN_STABLENORM_H
11 #define EIGEN_STABLENORM_H
17 template<
typename ExpressionType,
typename Scalar>
18 inline void stable_norm_kernel(
const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
20 Scalar max = bl.cwiseAbs().maxCoeff();
23 ssq = ssq * abs2(scale/max);
25 invScale = Scalar(1)/scale;
29 ssq += (bl*invScale).squaredNorm();
43 template<
typename Derived>
44 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
48 const Index blockSize = 4096;
50 RealScalar invScale(1);
56 Index bi = internal::first_aligned(derived());
58 internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
59 for (; bi<n; bi+=blockSize)
60 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).
template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
61 return scale * internal::sqrt(ssq);
73 template<
typename Derived>
80 static bool initialized =
false;
81 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
84 int ibeta, it, iemin, iemax, iexp;
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);
117 RealScalar ab2 = b2 / RealScalar(n);
118 RealScalar asml = RealScalar(0);
119 RealScalar amed = RealScalar(0);
120 RealScalar abig = RealScalar(0);
121 for(
Index j=0; j<n; ++j)
123 RealScalar ax = internal::abs(coeff(j));
124 if(ax > ab2) abig += internal::abs2(ax*s2m);
125 else if(ax < b1) asml += internal::abs2(ax*s1m);
126 else amed += internal::abs2(ax);
128 if(abig > RealScalar(0))
130 abig = internal::sqrt(abig);
135 if(amed > RealScalar(0))
138 amed = internal::sqrt(amed);
143 else if(asml > RealScalar(0))
145 if (amed > RealScalar(0))
147 abig = internal::sqrt(amed);
148 amed = internal::sqrt(asml) / s1m;
151 return internal::sqrt(asml)/s1m;
154 return internal::sqrt(amed);
155 asml = (min)(abig, amed);
156 abig = (max)(abig, amed);
157 if(asml <= abig*relerr)
160 return abig * internal::sqrt(RealScalar(1) + internal::abs2(asml/abig));
168 template<
typename Derived>
172 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
177 #endif // EIGEN_STABLENORM_H