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 * numext::abs2(scale/max);
25 invScale = Scalar(1)/scale;
29 ssq += (bl*invScale).squaredNorm();
32 template<
typename Derived>
33 inline typename NumTraits<typename traits<Derived>::Scalar>::Real
34 blueNorm_impl(
const EigenBase<Derived>& _vec)
36 typedef typename Derived::RealScalar RealScalar;
37 typedef typename Derived::Index Index;
43 const Derived& vec(_vec.derived());
44 static bool initialized =
false;
45 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
48 int ibeta, it, iemin, iemax, iexp;
58 ibeta = std::numeric_limits<RealScalar>::radix;
59 it = std::numeric_limits<RealScalar>::digits;
60 iemin = std::numeric_limits<RealScalar>::min_exponent;
61 iemax = std::numeric_limits<RealScalar>::max_exponent;
62 rbig = (std::numeric_limits<RealScalar>::max)();
64 iexp = -((1-iemin)/2);
65 b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
66 iexp = (iemax + 1 - it)/2;
67 b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
70 s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
71 iexp = - ((iemax+it)/2);
72 s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
75 eps = RealScalar(pow(
double(ibeta), 1-it));
80 RealScalar ab2 = b2 / RealScalar(n);
81 RealScalar asml = RealScalar(0);
82 RealScalar amed = RealScalar(0);
83 RealScalar abig = RealScalar(0);
84 for(
typename Derived::InnerIterator it(vec, 0); it; ++it)
86 RealScalar ax = abs(it.value());
87 if(ax > ab2) abig += numext::abs2(ax*s2m);
88 else if(ax < b1) asml += numext::abs2(ax*s1m);
89 else amed += numext::abs2(ax);
91 if(abig > RealScalar(0))
98 if(amed > RealScalar(0))
106 else if(asml > RealScalar(0))
108 if (amed > RealScalar(0))
111 amed = sqrt(asml) / s1m;
114 return sqrt(asml)/s1m;
118 asml = (min)(abig, amed);
119 abig = (max)(abig, amed);
120 if(asml <= abig*relerr)
123 return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig));
138 template<
typename Derived>
139 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
144 const Index blockSize = 4096;
146 RealScalar invScale(1);
152 Index bi = internal::first_aligned(derived());
154 internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
155 for (; bi<n; bi+=blockSize)
156 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).
template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
157 return scale * sqrt(ssq);
169 template<
typename Derived>
173 return internal::blueNorm_impl(*
this);
181 template<
typename Derived>
185 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
190 #endif // EIGEN_STABLENORM_H