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 lowTrunc = -Afrob + std::sqrt( Afrob*Afrob + threshold / Bfrob );
360 template<
typename Tmatrix,
typename TmatrixB,
typename Treal>
362 Treal
const threshold ) {
363 Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
366 if ( tmpInterval.
length() < 2*absTol ) {