25 #ifndef EIGEN_STABLENORM_H
26 #define EIGEN_STABLENORM_H
31 template<
typename ExpressionType,
typename Scalar>
32 inline void stable_norm_kernel(
const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
34 Scalar max = bl.cwiseAbs().maxCoeff();
37 ssq = ssq *
abs2(scale/max);
39 invScale = Scalar(1)/scale;
43 ssq += (bl*invScale).squaredNorm();
57 template<
typename Derived>
58 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
62 const Index blockSize = 4096;
70 Index bi = internal::first_aligned(derived());
73 for (; bi<n; bi+=blockSize)
74 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).
template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
87 template<
typename Derived>
94 static Index nmax = -1;
95 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
98 int nbig, ibeta, it, iemin, iemax, iexp;
108 nbig = (std::numeric_limits<Index>::max)();
109 ibeta = std::numeric_limits<RealScalar>::radix;
110 it = std::numeric_limits<RealScalar>::digits;
111 iemin = std::numeric_limits<RealScalar>::min_exponent;
112 iemax = std::numeric_limits<RealScalar>::max_exponent;
113 rbig = (std::numeric_limits<RealScalar>::max)();
115 iexp = -((1-iemin)/2);
117 iexp = (iemax + 1 - it)/2;
122 iexp = - ((iemax+it)/2);
137 for(
Index j=0; j<n; ++j)
172 asml = (min)(abig, amed);
173 abig = (max)(abig, amed);
174 if(asml <= abig*relerr)
185 template<
typename Derived>
189 return this->
cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
194 #endif // EIGEN_STABLENORM_H