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