Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_STABLENORM_H
00026 #define EIGEN_STABLENORM_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031 template<typename ExpressionType, typename Scalar>
00032 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
00033 {
00034 Scalar max = bl.cwiseAbs().maxCoeff();
00035 if (max>scale)
00036 {
00037 ssq = ssq * abs2(scale/max);
00038 scale = max;
00039 invScale = Scalar(1)/scale;
00040 }
00041
00042
00043 ssq += (bl*invScale).squaredNorm();
00044 }
00045 }
00046
00057 template<typename Derived>
00058 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00059 MatrixBase<Derived>::stableNorm() const
00060 {
00061 using std::min;
00062 const Index blockSize = 4096;
00063 RealScalar scale(0);
00064 RealScalar invScale(1);
00065 RealScalar ssq(0);
00066 enum {
00067 Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? 1 : 0
00068 };
00069 Index n = size();
00070 Index bi = internal::first_aligned(derived());
00071 if (bi>0)
00072 internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
00073 for (; bi<n; bi+=blockSize)
00074 internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
00075 return scale * internal::sqrt(ssq);
00076 }
00077
00087 template<typename Derived>
00088 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00089 MatrixBase<Derived>::blueNorm() const
00090 {
00091 using std::pow;
00092 using std::min;
00093 using std::max;
00094 static Index nmax = -1;
00095 static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
00096 if(nmax <= 0)
00097 {
00098 int nbig, ibeta, it, iemin, iemax, iexp;
00099 RealScalar abig, eps;
00100
00101
00102
00103
00104
00105
00106
00107
00108 nbig = (std::numeric_limits<Index>::max)();
00109 ibeta = std::numeric_limits<RealScalar>::radix;
00110 it = std::numeric_limits<RealScalar>::digits;
00111 iemin = std::numeric_limits<RealScalar>::min_exponent;
00112 iemax = std::numeric_limits<RealScalar>::max_exponent;
00113 rbig = (std::numeric_limits<RealScalar>::max)();
00114
00115 iexp = -((1-iemin)/2);
00116 b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
00117 iexp = (iemax + 1 - it)/2;
00118 b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
00119
00120 iexp = (2-iemin)/2;
00121 s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
00122 iexp = - ((iemax+it)/2);
00123 s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));
00124
00125 overfl = rbig*s2m;
00126 eps = RealScalar(pow(double(ibeta), 1-it));
00127 relerr = internal::sqrt(eps);
00128 abig = RealScalar(1.0/eps - 1.0);
00129 if (RealScalar(nbig)>abig) nmax = int(abig);
00130 else nmax = nbig;
00131 }
00132 Index n = size();
00133 RealScalar ab2 = b2 / RealScalar(n);
00134 RealScalar asml = RealScalar(0);
00135 RealScalar amed = RealScalar(0);
00136 RealScalar abig = RealScalar(0);
00137 for(Index j=0; j<n; ++j)
00138 {
00139 RealScalar ax = internal::abs(coeff(j));
00140 if(ax > ab2) abig += internal::abs2(ax*s2m);
00141 else if(ax < b1) asml += internal::abs2(ax*s1m);
00142 else amed += internal::abs2(ax);
00143 }
00144 if(abig > RealScalar(0))
00145 {
00146 abig = internal::sqrt(abig);
00147 if(abig > overfl)
00148 {
00149 eigen_assert(false && "overflow");
00150 return rbig;
00151 }
00152 if(amed > RealScalar(0))
00153 {
00154 abig = abig/s2m;
00155 amed = internal::sqrt(amed);
00156 }
00157 else
00158 return abig/s2m;
00159 }
00160 else if(asml > RealScalar(0))
00161 {
00162 if (amed > RealScalar(0))
00163 {
00164 abig = internal::sqrt(amed);
00165 amed = internal::sqrt(asml) / s1m;
00166 }
00167 else
00168 return internal::sqrt(asml)/s1m;
00169 }
00170 else
00171 return internal::sqrt(amed);
00172 asml = (min)(abig, amed);
00173 abig = (max)(abig, amed);
00174 if(asml <= abig*relerr)
00175 return abig;
00176 else
00177 return abig * internal::sqrt(RealScalar(1) + internal::abs2(asml/abig));
00178 }
00179
00185 template<typename Derived>
00186 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
00187 MatrixBase<Derived>::hypotNorm() const
00188 {
00189 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>());
00190 }
00191
00192 }
00193
00194 #endif // EIGEN_STABLENORM_H