41 #ifndef MAT_TRUNCATION
42 #define MAT_TRUNCATION
50 template<
typename Tmatrix,
typename Treal>
54 Treal
run(Treal
const threshold);
58 Treal
const threshold ) = 0;
59 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms ) = 0;
62 Treal
const threshold ) = 0;
67 template<
typename Tmatrix,
typename Treal>
74 E.resetSizesAndBlocks(rows, cols);
77 template<
typename Tmatrix,
typename Treal>
79 assert(threshold >= (Treal)0.0);
80 if (threshold == (Treal)0.0)
82 std::vector<Treal> frobsq_norms;
83 this->getFrobSqNorms( frobsq_norms );
84 std::sort(frobsq_norms.begin(), frobsq_norms.end());
86 int high = frobsq_norms.size() - 1;
87 Treal lowFrobTrunc, highFrobTrunc;
88 this->getFrobTruncBounds( lowFrobTrunc, highFrobTrunc, threshold );
90 while( low < (
int)frobsq_norms.size() - 1 && frobsqSum < lowFrobTrunc ) {
92 frobsqSum += frobsq_norms[low];
96 while( high < (
int)frobsq_norms.size() - 1 && frobsqSum < highFrobTrunc ) {
98 frobsqSum += frobsq_norms[high];
101 int minStep = int( 0.01 * frobsq_norms.size() );
102 minStep = minStep > 0 ? minStep : 1;
103 int testIndex = high;
104 int previousTestIndex = high * 2;
108 while ( euclEInt.
upp() > threshold ) {
111 int stepSize = (int)((high - low) * 0.01);
113 stepSize = stepSize >= minStep ? stepSize : minStep;
114 previousTestIndex = testIndex;
115 testIndex -= stepSize;
117 testIndex = testIndex > low ? testIndex : low;
122 while(testIndex >= 0 && frobsq_norms[testIndex] == frobsq_norms[testIndex+1])
129 assert( previousTestIndex != testIndex );
130 Treal currentFrobTrunc = frobsq_norms[testIndex];
131 frobThreshLowLevel( currentFrobTrunc );
132 euclEInt =
euclIfSmall( Treal(threshold * 1e-2), threshold );
136 if ( testIndex <= -1 ) {
137 frobThreshLowLevel( (Treal)0.0 );
141 euclE = euclEInt.
upp();
151 template<
typename Tmatrix,
typename Treal>
158 Treal
const threshold );
162 Treal
const threshold );
165 template<
typename Tmatrix,
typename Treal>
167 Treal
const threshold ) {
169 lowTrunc = (threshold * threshold) / 2;
170 highTrunc = (threshold * threshold * this->
A.get_nrows()) / 2;
173 template<
typename Tmatrix,
typename Treal>
175 this->
A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
178 template<
typename Tmatrix,
typename Treal>
180 this->
A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
183 template<
typename Tmatrix,
typename Treal>
185 Treal
const threshold ) {
186 Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
188 if ( tmpInterval.
length() < 2*absTol )
200 template<
typename Tmatrix,
typename TmatrixZ,
typename Treal>
207 Treal
const threshold );
211 Treal
const threshold );
215 template<
typename Tmatrix,
typename TmatrixZ,
typename Treal>
217 Treal
const threshold ) {
218 Treal Zfrob = Z.frob();
219 Treal thresholdTakingZIntoAccount = threshold / (Zfrob * Zfrob);
221 lowTrunc = thresholdTakingZIntoAccount * thresholdTakingZIntoAccount / 2.0;
225 template<
typename Tmatrix,
typename TmatrixZ,
typename Treal>
227 Treal
const threshold ) {
228 Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
231 if ( tmpInterval.
length() < 2*absTol )
243 template<
typename Tmatrix,
typename Treal>
255 template<
typename Tmatrix,
typename Treal>
257 this->
A.getMatrix().getFrobSqElementLevel(frobsq_norms);
260 template<
typename Tmatrix,
typename Treal>
262 this->
A.getMatrix().frobThreshElementLevel(threshold, &this->E.getMatrix() );
269 template<
typename Tmatrix,
typename Treal>
276 Treal
const threshold );
280 Treal
const threshold );
283 template<
typename Tmatrix,
typename Treal>
285 Treal
const threshold ) {
291 lowTrunc = (threshold * threshold);
296 highTrunc = (threshold * threshold * this->
A.get_nrows());
299 template<
typename Tmatrix,
typename Treal>
301 this->
A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
304 template<
typename Tmatrix,
typename Treal>
306 this->
A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
309 template<
typename Tmatrix,
typename Treal>
311 Treal
const threshold ) {
317 Treal relTol = 100 * std::numeric_limits<Treal>::epsilon();
320 if ( tmpInterval.
length() < 2*absTol )
335 template<
typename Tmatrix,
typename TmatrixB,
typename Treal>
342 Treal
const threshold );
346 Treal
const threshold );
350 template<
typename Tmatrix,
typename TmatrixB,
typename Treal>
353 Treal
const threshold ) {
354 Treal Afrob = this->
A.frob();
355 Treal Bfrob =
B.frob();
356 Treal tmp = -Afrob + std::sqrt( Afrob*Afrob + threshold / Bfrob );
361 template<
typename Tmatrix,
typename TmatrixB,
typename Treal>
363 Treal
const threshold ) {
364 Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
367 if ( tmpInterval.
length() < 2*absTol ) {
virtual void frobThreshLowLevel(Treal const threshold)
Definition: truncation.h:179
EuclTruncationBase(Tmatrix &A_)
Definition: truncation.h:68
Treal upp() const
Definition: Interval.h:143
Tmatrix & A
Definition: truncation.h:63
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
TmatrixZ const & Z
Definition: truncation.h:212
EuclTruncationSymm(Tmatrix &A_)
Definition: truncation.h:154
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition: truncation.h:174
Truncation of symmetric matrices at the element level (used for mixed norm truncation) ...
Definition: truncation.h:244
EuclTruncationSymmWithZ(Tmatrix &A_, TmatrixZ const &Z_)
Definition: truncation.h:203
Truncation of general matrices.
Definition: truncation.h:270
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition: truncation.h:284
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)=0
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition: truncation.h:216
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:37
Treal midPoint() const
Definition: Interval.h:113
Definition: Interval.h:44
EuclTruncationCongrTransMeasure(Tmatrix &A_, TmatrixB const &B_)
Definition: truncation.h:338
EuclTruncationSymmElementLevel(Tmatrix &A_)
Definition: truncation.h:246
EuclTruncationGeneral(Tmatrix &A_)
Definition: truncation.h:272
#define max(a, b)
Definition: integrator.cc:92
virtual void frobThreshLowLevel(Treal const threshold)
Definition: truncation.h:261
TmatrixB const & B
Definition: truncation.h:347
virtual void frobThreshLowLevel(Treal const threshold)=0
Treal low() const
Definition: Interval.h:142
Treal run(Treal const threshold)
Definition: truncation.h:78
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition: truncation.h:256
Interval< Treal > euclIfSmall(Tmatrix const &M, Treal const requestedAbsAccuracy, Treal const requestedRelAccuracy, Treal const maxAbsVal, typename Tmatrix::VectorType *const eVecPtr=0)
Returns interval containing the Euclidean norm of the matrix M.
Definition: LanczosLargestMagnitudeEig.h:257
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition: truncation.h:362
Definition: mat_utils.h:97
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition: truncation.h:351
Definition: mat_utils.h:128
Definition: mat_utils.h:69
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)=0
Definition: truncation.h:51
virtual void frobThreshLowLevel(Treal const threshold)
Definition: truncation.h:305
virtual ~EuclTruncationBase()
Definition: truncation.h:55
Tmatrix E
Definition: truncation.h:64
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)=0
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition: truncation.h:184
Truncation of symmetric matrices.
Definition: truncation.h:152
Truncation of symmetric matrices with Z.
Definition: truncation.h:201
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition: truncation.h:310
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition: truncation.h:300
Truncation of general matrices with impact on matrix triple multiply as error measure.
Definition: truncation.h:336
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition: truncation.h:226
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition: truncation.h:166