![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_STABLENORM_H 00011 #define EIGEN_STABLENORM_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 00017 template<typename ExpressionType, typename Scalar> 00018 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) 00019 { 00020 Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); 00021 00022 if(maxCoeff>scale) 00023 { 00024 ssq = ssq * numext::abs2(scale/maxCoeff); 00025 Scalar tmp = Scalar(1)/maxCoeff; 00026 if(tmp > NumTraits<Scalar>::highest()) 00027 { 00028 invScale = NumTraits<Scalar>::highest(); 00029 scale = Scalar(1)/invScale; 00030 } 00031 else if(maxCoeff>NumTraits<Scalar>::highest()) // we got a INF 00032 { 00033 invScale = Scalar(1); 00034 scale = maxCoeff; 00035 } 00036 else 00037 { 00038 scale = maxCoeff; 00039 invScale = tmp; 00040 } 00041 } 00042 else if(maxCoeff!=maxCoeff) // we got a NaN 00043 { 00044 scale = maxCoeff; 00045 } 00046 00047 // TODO if the maxCoeff is much much smaller than the current scale, 00048 // then we can neglect this sub vector 00049 if(scale>Scalar(0)) // if scale==0, then bl is 0 00050 ssq += (bl*invScale).squaredNorm(); 00051 } 00052 00053 template<typename Derived> 00054 inline typename NumTraits<typename traits<Derived>::Scalar>::Real 00055 blueNorm_impl(const EigenBase<Derived>& _vec) 00056 { 00057 typedef typename Derived::RealScalar RealScalar; 00058 using std::pow; 00059 using std::sqrt; 00060 using std::abs; 00061 const Derived& vec(_vec.derived()); 00062 static bool initialized = false; 00063 static RealScalar b1, b2, s1m, s2m, rbig, relerr; 00064 if(!initialized) 00065 { 00066 int ibeta, it, iemin, iemax, iexp; 00067 RealScalar eps; 00068 // This program calculates the machine-dependent constants 00069 // bl, b2, slm, s2m, relerr overfl 00070 // from the "basic" machine-dependent numbers 00071 // nbig, ibeta, it, iemin, iemax, rbig. 00072 // The following define the basic machine-dependent constants. 00073 // For portability, the PORT subprograms "ilmaeh" and "rlmach" 00074 // are used. For any specific computer, each of the assignment 00075 // statements can be replaced 00076 ibeta = std::numeric_limits<RealScalar>::radix; // base for floating-point numbers 00077 it = std::numeric_limits<RealScalar>::digits; // number of base-beta digits in mantissa 00078 iemin = std::numeric_limits<RealScalar>::min_exponent; // minimum exponent 00079 iemax = std::numeric_limits<RealScalar>::max_exponent; // maximum exponent 00080 rbig = (std::numeric_limits<RealScalar>::max)(); // largest floating-point number 00081 00082 iexp = -((1-iemin)/2); 00083 b1 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // lower boundary of midrange 00084 iexp = (iemax + 1 - it)/2; 00085 b2 = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // upper boundary of midrange 00086 00087 iexp = (2-iemin)/2; 00088 s1m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for lower range 00089 iexp = - ((iemax+it)/2); 00090 s2m = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp))); // scaling factor for upper range 00091 00092 eps = RealScalar(pow(double(ibeta), 1-it)); 00093 relerr = sqrt(eps); // tolerance for neglecting asml 00094 initialized = true; 00095 } 00096 Index n = vec.size(); 00097 RealScalar ab2 = b2 / RealScalar(n); 00098 RealScalar asml = RealScalar(0); 00099 RealScalar amed = RealScalar(0); 00100 RealScalar abig = RealScalar(0); 00101 for(typename Derived::InnerIterator it(vec, 0); it; ++it) 00102 { 00103 RealScalar ax = abs(it.value()); 00104 if(ax > ab2) abig += numext::abs2(ax*s2m); 00105 else if(ax < b1) asml += numext::abs2(ax*s1m); 00106 else amed += numext::abs2(ax); 00107 } 00108 if(amed!=amed) 00109 return amed; // we got a NaN 00110 if(abig > RealScalar(0)) 00111 { 00112 abig = sqrt(abig); 00113 if(abig > rbig) // overflow, or *this contains INF values 00114 return abig; // return INF 00115 if(amed > RealScalar(0)) 00116 { 00117 abig = abig/s2m; 00118 amed = sqrt(amed); 00119 } 00120 else 00121 return abig/s2m; 00122 } 00123 else if(asml > RealScalar(0)) 00124 { 00125 if (amed > RealScalar(0)) 00126 { 00127 abig = sqrt(amed); 00128 amed = sqrt(asml) / s1m; 00129 } 00130 else 00131 return sqrt(asml)/s1m; 00132 } 00133 else 00134 return sqrt(amed); 00135 asml = numext::mini(abig, amed); 00136 abig = numext::maxi(abig, amed); 00137 if(asml <= abig*relerr) 00138 return abig; 00139 else 00140 return abig * sqrt(RealScalar(1) + numext::abs2(asml/abig)); 00141 } 00142 00143 } // end namespace internal 00144 00155 template<typename Derived> 00156 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 00157 MatrixBase<Derived>::stableNorm() const 00158 { 00159 using std::sqrt; 00160 using std::abs; 00161 const Index blockSize = 4096; 00162 RealScalar scale(0); 00163 RealScalar invScale(1); 00164 RealScalar ssq(0); // sum of square 00165 00166 typedef typename internal::nested_eval<Derived,2>::type DerivedCopy; 00167 typedef typename internal::remove_all<DerivedCopy>::type DerivedCopyClean; 00168 DerivedCopy copy(derived()); 00169 00170 enum { 00171 CanAlign = ( (int(DerivedCopyClean::Flags)&DirectAccessBit) 00172 || (int(internal::evaluator<DerivedCopyClean>::Alignment)>0) // FIXME Alignment)>0 might not be enough 00173 ) && (blockSize*sizeof(Scalar)*2<EIGEN_STACK_ALLOCATION_LIMIT) // ifwe cannot allocate on the stack, then let's not bother about this optimization 00174 }; 00175 typedef typename internal::conditional<CanAlign, Ref<const Matrix<Scalar,Dynamic,1,0,blockSize,1>, internal::evaluator<DerivedCopyClean>::Alignment>, 00176 typename DerivedCopyClean::ConstSegmentReturnType>::type SegmentWrapper; 00177 Index n = size(); 00178 00179 if(n==1) 00180 return abs(this->coeff(0)); 00181 00182 Index bi = internal::first_default_aligned(copy); 00183 if (bi>0) 00184 internal::stable_norm_kernel(copy.head(bi), ssq, scale, invScale); 00185 for (; bi<n; bi+=blockSize) 00186 internal::stable_norm_kernel(SegmentWrapper(copy.segment(bi,numext::mini(blockSize, n - bi))), ssq, scale, invScale); 00187 return scale * sqrt(ssq); 00188 } 00189 00199 template<typename Derived> 00200 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 00201 MatrixBase<Derived>::blueNorm() const 00202 { 00203 return internal::blueNorm_impl(*this); 00204 } 00205 00211 template<typename Derived> 00212 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real 00213 MatrixBase<Derived>::hypotNorm() const 00214 { 00215 return this->cwiseAbs().redux(internal::scalar_hypot_op<RealScalar>()); 00216 } 00217 00218 } // end namespace Eigen 00219 00220 #endif // EIGEN_STABLENORM_H