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
00026
00027
00028
00029
00043 #ifndef MAT_TRUNCATION
00044 #define MAT_TRUNCATION
00045 #include <limits>
00046 #include <stdexcept>
00047 #include <cmath>
00048 namespace mat {
00049
00050
00051
00052 template<typename Tmatrix, typename Treal>
00053 class EuclTruncationBase {
00054 public:
00055 explicit EuclTruncationBase( Tmatrix & A_ );
00056 Treal run(Treal const threshold);
00057 virtual ~EuclTruncationBase() {}
00058 protected:
00059 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00060 Treal const threshold ) = 0;
00061 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms ) = 0;
00062 virtual void frobThreshLowLevel( Treal const threshold ) = 0;
00063 virtual Interval<Treal> euclIfSmall( Treal const absTol,
00064 Treal const threshold ) = 0;
00065 Tmatrix & A;
00066 Tmatrix E;
00067 };
00068
00069 template<typename Tmatrix, typename Treal>
00070 EuclTruncationBase<Tmatrix, Treal>::EuclTruncationBase( Tmatrix & A_ )
00071 : A(A_) {
00072 SizesAndBlocks rows;
00073 SizesAndBlocks cols;
00074 A.getRows(rows);
00075 A.getCols(cols);
00076 E.resetSizesAndBlocks(rows, cols);
00077 }
00078
00079 template<typename Tmatrix, typename Treal>
00080 Treal EuclTruncationBase<Tmatrix, Treal>::run( Treal const threshold ) {
00081 assert(threshold >= (Treal)0.0);
00082 if (threshold == (Treal)0.0)
00083 return (Treal)0;
00084 std::vector<Treal> frobsq_norms;
00085 this->getFrobSqNorms( frobsq_norms );
00086 std::sort(frobsq_norms.begin(), frobsq_norms.end());
00087 int low = -1;
00088 int high = frobsq_norms.size() - 1;
00089 Treal lowFrobTrunc, highFrobTrunc;
00090 this->getFrobTruncBounds( lowFrobTrunc, highFrobTrunc, threshold );
00091 Treal frobsqSum = 0;
00092 while( low < (int)frobsq_norms.size() - 1 && frobsqSum < lowFrobTrunc ) {
00093 ++low;
00094 frobsqSum += frobsq_norms[low];
00095 }
00096 high = low;
00097 --low;
00098 while( high < (int)frobsq_norms.size() - 1 && frobsqSum < highFrobTrunc ) {
00099 ++high;
00100 frobsqSum += frobsq_norms[high];
00101 }
00102
00103 int minStep = int( 0.01 * frobsq_norms.size() );
00104 minStep = minStep > 0 ? minStep : 1;
00105 int testIndex = high;
00106 int previousTestIndex = high * 2;
00107
00108 Interval<Treal> euclEInt(0, threshold * 2);
00109
00110 while ( euclEInt.upp() > threshold ) {
00111
00112 high = testIndex;
00113 int stepSize = (int)((high - low) * 0.01);
00114
00115 stepSize = stepSize >= minStep ? stepSize : minStep;
00116 previousTestIndex = testIndex;
00117 testIndex -= stepSize;
00118
00119 testIndex = testIndex > low ? testIndex : low;
00120
00121
00122
00123
00124 while(testIndex >= 0 && frobsq_norms[testIndex] == frobsq_norms[testIndex+1])
00125 testIndex--;
00126
00127
00128 if ( testIndex < 0 )
00129
00130 break;
00131 assert( previousTestIndex != testIndex );
00132 Treal currentFrobTrunc = frobsq_norms[testIndex];
00133 frobThreshLowLevel( currentFrobTrunc );
00134 euclEInt = euclIfSmall( Treal(threshold * 1e-2), threshold );
00135
00136 }
00137 Treal euclE;
00138 if ( testIndex <= -1 ) {
00139 frobThreshLowLevel( (Treal)0.0 );
00140 euclE = 0;
00141 }
00142 else {
00143 euclE = euclEInt.upp();
00144 }
00145 return euclE;
00146 }
00147
00148
00153 template<typename Tmatrix, typename Treal>
00154 class EuclTruncationSymm : public EuclTruncationBase<Tmatrix, Treal> {
00155 public:
00156 explicit EuclTruncationSymm( Tmatrix & A_ )
00157 : EuclTruncationBase<Tmatrix, Treal>(A_) {}
00158 protected:
00159 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00160 Treal const threshold );
00161 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00162 virtual void frobThreshLowLevel( Treal const threshold );
00163 virtual Interval<Treal> euclIfSmall( Treal const absTol,
00164 Treal const threshold );
00165 };
00166
00167 template<typename Tmatrix, typename Treal>
00168 void EuclTruncationSymm<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00169 Treal const threshold ) {
00170
00171 lowTrunc = (threshold * threshold) / 2;
00172 highTrunc = (threshold * threshold * this->A.get_nrows()) / 2;
00173 }
00174
00175 template<typename Tmatrix, typename Treal>
00176 void EuclTruncationSymm<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00177 this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
00178 }
00179
00180 template<typename Tmatrix, typename Treal>
00181 void EuclTruncationSymm<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00182 this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
00183 }
00184
00185 template<typename Tmatrix, typename Treal>
00186 Interval<Treal> EuclTruncationSymm<Tmatrix, Treal>::euclIfSmall( Treal const absTol,
00187 Treal const threshold ) {
00188 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
00189 Interval<Treal> tmpInterval = mat::euclIfSmall(this->E, absTol, relTol, threshold);
00190 if ( tmpInterval.length() < 2*absTol )
00191 return Interval<Treal>( tmpInterval.midPoint()-absTol,
00192 tmpInterval.midPoint()+absTol );
00193 return tmpInterval;
00194 }
00195
00202 template<typename Tmatrix, typename TmatrixZ, typename Treal>
00203 class EuclTruncationSymmWithZ : public EuclTruncationSymm<Tmatrix, Treal> {
00204 public:
00205 EuclTruncationSymmWithZ( Tmatrix & A_, TmatrixZ const & Z_ )
00206 : EuclTruncationSymm<Tmatrix, Treal>(A_), Z(Z_) {}
00207 protected:
00208 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00209 Treal const threshold );
00210
00211
00212 virtual Interval<Treal> euclIfSmall( Treal const absTol,
00213 Treal const threshold );
00214 TmatrixZ const & Z;
00215 };
00216
00217 template<typename Tmatrix, typename TmatrixZ, typename Treal>
00218 void EuclTruncationSymmWithZ<Tmatrix, TmatrixZ, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00219 Treal const threshold ) {
00220 Treal Zfrob = Z.frob();
00221 Treal thresholdTakingZIntoAccount = threshold / (Zfrob * Zfrob);
00222
00223 lowTrunc = thresholdTakingZIntoAccount * thresholdTakingZIntoAccount / 2.0;
00224 highTrunc = template_blas_get_num_limit_max<Treal>();
00225 }
00226
00227 template<typename Tmatrix, typename TmatrixZ, typename Treal>
00228 Interval<Treal> EuclTruncationSymmWithZ<Tmatrix, TmatrixZ, Treal>::euclIfSmall( Treal const absTol,
00229 Treal const threshold ) {
00230 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
00231 mat::TripleMatrix<Tmatrix, TmatrixZ, Treal> ErrMatTriple( this->E, Z);
00232 Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);
00233 if ( tmpInterval.length() < 2*absTol )
00234 return Interval<Treal>( tmpInterval.midPoint()-absTol,
00235 tmpInterval.midPoint()+absTol );
00236 return tmpInterval;
00237 }
00238
00245 template<typename Tmatrix, typename Treal>
00246 class EuclTruncationSymmElementLevel : public EuclTruncationSymm<Tmatrix, Treal> {
00247 public:
00248 explicit EuclTruncationSymmElementLevel( Tmatrix & A_ )
00249 : EuclTruncationSymm<Tmatrix, Treal>(A_) {}
00250 protected:
00251
00252 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00253 virtual void frobThreshLowLevel( Treal const threshold );
00254
00255 };
00256
00257 template<typename Tmatrix, typename Treal>
00258 void EuclTruncationSymmElementLevel<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00259 this->A.getMatrix().getFrobSqElementLevel(frobsq_norms);
00260 }
00261
00262 template<typename Tmatrix, typename Treal>
00263 void EuclTruncationSymmElementLevel<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00264 this->A.getMatrix().frobThreshElementLevel(threshold, &this->E.getMatrix() );
00265 }
00266
00271 template<typename Tmatrix, typename Treal>
00272 class EuclTruncationGeneral : public EuclTruncationBase<Tmatrix, Treal> {
00273 public:
00274 explicit EuclTruncationGeneral( Tmatrix & A_ )
00275 : EuclTruncationBase<Tmatrix, Treal>(A_) {}
00276 protected:
00277 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00278 Treal const threshold );
00279 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00280 virtual void frobThreshLowLevel( Treal const threshold );
00281 virtual Interval<Treal> euclIfSmall( Treal const absTol,
00282 Treal const threshold );
00283 };
00284
00285 template<typename Tmatrix, typename Treal>
00286 void EuclTruncationGeneral<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00287 Treal const threshold ) {
00288
00289
00290
00291
00292
00293 lowTrunc = (threshold * threshold);
00294
00295
00296
00297
00298 highTrunc = (threshold * threshold * this->A.get_nrows());
00299 }
00300
00301 template<typename Tmatrix, typename Treal>
00302 void EuclTruncationGeneral<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00303 this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
00304 }
00305
00306 template<typename Tmatrix, typename Treal>
00307 void EuclTruncationGeneral<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00308 this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
00309 }
00310
00311 template<typename Tmatrix, typename Treal>
00312 Interval<Treal> EuclTruncationGeneral<Tmatrix, Treal>::euclIfSmall( Treal const absTol,
00313 Treal const threshold ) {
00314
00315
00316
00317 mat::ATAMatrix<Tmatrix, Treal> EtE(this->E);
00318 Treal absTolDummy = template_blas_get_num_limit_max<Treal>();
00319 Treal relTol = 100 * mat::getMachineEpsilon<Treal>();
00320 Interval<Treal> tmpInterval = mat::euclIfSmall(EtE, absTolDummy, relTol, threshold);
00321 tmpInterval = Interval<Treal>( template_blas_sqrt(tmpInterval.low()), template_blas_sqrt(tmpInterval.upp()) );
00322 if ( tmpInterval.length() < 2*absTol )
00323 return Interval<Treal>( tmpInterval.midPoint()-absTol,
00324 tmpInterval.midPoint()+absTol );
00325 return tmpInterval;
00326 }
00327
00328
00329
00330
00337 template<typename Tmatrix, typename TmatrixB, typename Treal>
00338 class EuclTruncationCongrTransMeasure : public EuclTruncationGeneral<Tmatrix, Treal> {
00339 public:
00340 EuclTruncationCongrTransMeasure( Tmatrix & A_, TmatrixB const & B_ )
00341 : EuclTruncationGeneral<Tmatrix, Treal>(A_), B(B_) {}
00342 protected:
00343 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
00344 Treal const threshold );
00345
00346
00347 virtual Interval<Treal> euclIfSmall( Treal const absTol,
00348 Treal const threshold );
00349 TmatrixB const & B;
00350 };
00351
00352 template<typename Tmatrix, typename TmatrixB, typename Treal>
00353 void EuclTruncationCongrTransMeasure<Tmatrix, TmatrixB, Treal>::getFrobTruncBounds( Treal & lowTrunc,
00354 Treal & highTrunc,
00355 Treal const threshold ) {
00356 Treal Afrob = this->A.frob();
00357 Treal Bfrob = B.frob();
00358 Treal tmp = -Afrob + template_blas_sqrt( Afrob*Afrob + threshold / Bfrob );
00359 lowTrunc = tmp*tmp;
00360 highTrunc = template_blas_get_num_limit_max<Treal>();
00361 }
00362
00363 template<typename Tmatrix, typename TmatrixB, typename Treal>
00364 Interval<Treal> EuclTruncationCongrTransMeasure<Tmatrix, TmatrixB, Treal>::euclIfSmall( Treal const absTol,
00365 Treal const threshold ) {
00366 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
00367 mat::CongrTransErrorMatrix<TmatrixB, Tmatrix, Treal> ErrMatTriple( B, this->A, this->E );
00368 Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);
00369 if ( tmpInterval.length() < 2*absTol ) {
00370 return Interval<Treal>( tmpInterval.midPoint()-absTol,
00371 tmpInterval.midPoint()+absTol );
00372 }
00373 return tmpInterval;
00374 }
00375
00376
00377 }
00378 #endif