ergo
truncation.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
43#ifndef MAT_TRUNCATION
44#define MAT_TRUNCATION
45#include <limits>
46#include <stdexcept>
47#include <cmath>
48namespace mat { /* Matrix namespace */
49
50 // Stuff for Euclidean norm based truncation
51
52 template<typename Tmatrix, typename Treal>
54 public:
55 explicit EuclTruncationBase( Tmatrix & A_ );
56 Treal run(Treal const threshold);
58 protected:
59 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
60 Treal const threshold ) = 0;
61 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms ) = 0;
62 virtual void frobThreshLowLevel( Treal const threshold ) = 0;
63 virtual Interval<Treal> euclIfSmall( Treal const absTol,
64 Treal const threshold ) = 0;
65 Tmatrix & A; // Matrix to be truncated
66 Tmatrix E; // Error matrix
67 };
68
69 template<typename Tmatrix, typename Treal>
71 : A(A_) {
74 A.getRows(rows);
75 A.getCols(cols);
76 E.resetSizesAndBlocks(rows, cols);
77 }
78
79 template<typename Tmatrix, typename Treal>
80 Treal EuclTruncationBase<Tmatrix, Treal>::run( Treal const threshold ) {
81 assert(threshold >= (Treal)0.0);
82 if (threshold == (Treal)0.0)
83 return (Treal)0;
84 std::vector<Treal> frobsq_norms;
85 this->getFrobSqNorms( frobsq_norms ); /*=======*/
86 std::sort(frobsq_norms.begin(), frobsq_norms.end());
87 int low = -1;
88 int high = frobsq_norms.size() - 1;
89 Treal lowFrobTrunc, highFrobTrunc;
90 this->getFrobTruncBounds( lowFrobTrunc, highFrobTrunc, threshold ); /*=======*/
91 Treal frobsqSum = 0;
92 while( low < (int)frobsq_norms.size() - 1 && frobsqSum < lowFrobTrunc ) {
93 ++low;
94 frobsqSum += frobsq_norms[low];
95 }
96 high = low; /* Removing all tom high is to much */
97 --low;
98 while( high < (int)frobsq_norms.size() - 1 && frobsqSum < highFrobTrunc ) {
99 ++high;
100 frobsqSum += frobsq_norms[high];
101 }
102 // Now we have low and high
103 int minStep = int( 0.01 * frobsq_norms.size() ); // Consider elements in chunks of at least 1 percent of all elements at a time to not get too many iterations
104 minStep = minStep > 0 ? minStep : 1; // step is at least one
105 int testIndex = high;
106 int previousTestIndex = high * 2;
107 // Now, removing everything up to and including testIndex is too much
108 Interval<Treal> euclEInt(0, threshold * 2);
109 // We go from above (too many elements in the error matrix) and stop as soon as the error matrix is small enough
110 while ( euclEInt.upp() > threshold ) {
111 // Removing everything up to and including testIndex is too much, update high:
112 high = testIndex;
113 int stepSize = (int)((high - low) * 0.01); // We can accept that only 99% of elements possible to remove are removed
114 // stepSize must be at least minStep:
115 stepSize = stepSize >= minStep ? stepSize : minStep;
116 previousTestIndex = testIndex;
117 testIndex -= stepSize;
118 // testIndex cannot be smaller than low
119 testIndex = testIndex > low ? testIndex : low;
120 /* Now we have decided the testIndex we would like to
121 use. However, we must be careful to handle the case when
122 there are several identical values in the frobsq_norms
123 list. In that case we use a modified value. */
124 while(testIndex >= 0 && frobsq_norms[testIndex] == frobsq_norms[testIndex+1])
125 testIndex--;
126 /* Note that because of the above while loop, at this point it
127 is possible that testIndex < low. */
128 if ( testIndex < 0 )
129 // testIndex == -1, we have to break
130 break;
131 assert( previousTestIndex != testIndex );
132 Treal currentFrobTrunc = frobsq_norms[testIndex];
133 frobThreshLowLevel( currentFrobTrunc ); /*=======*/
134 euclEInt = euclIfSmall( Treal(threshold * 1e-2), threshold ); /*=======*/
135 // Now we have an interval containing the Euclidean norm of E for the current testIndex
136 } // end while
137 Treal euclE;
138 if ( testIndex <= -1 ) {
139 frobThreshLowLevel( (Treal)0.0 ); /*=======*/
140 euclE = 0;
141 }
142 else {
143 euclE = euclEInt.upp();
144 }
145 return euclE;
146 } // end run
147
148
153 template<typename Tmatrix, typename Treal>
154 class EuclTruncationSymm : public EuclTruncationBase<Tmatrix, Treal> {
155 public:
156 explicit EuclTruncationSymm( Tmatrix & A_ )
157 : EuclTruncationBase<Tmatrix, Treal>(A_) {}
158 protected:
159 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
160 Treal const threshold );
161 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
162 virtual void frobThreshLowLevel( Treal const threshold );
163 virtual Interval<Treal> euclIfSmall( Treal const absTol,
164 Treal const threshold );
165 }; // end class EuclTruncationSymm
166
167 template<typename Tmatrix, typename Treal>
168 void EuclTruncationSymm<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
169 Treal const threshold ) {
170 /* Divide by 2 because of symmetry */
171 lowTrunc = (threshold * threshold) / 2;
172 highTrunc = (threshold * threshold * this->A.get_nrows()) / 2;
173 }
174
175 template<typename Tmatrix, typename Treal>
176 void EuclTruncationSymm<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
177 this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
178 }
179
180 template<typename Tmatrix, typename Treal>
182 this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
183 }
184
185 template<typename Tmatrix, typename Treal>
187 Treal const threshold ) {
188 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
189 Interval<Treal> tmpInterval = mat::euclIfSmall(this->E, absTol, relTol, threshold);
190 if ( tmpInterval.length() < 2*absTol )
191 return Interval<Treal>( tmpInterval.midPoint()-absTol,
192 tmpInterval.midPoint()+absTol );
193 return tmpInterval;
194 }
195
202 template<typename Tmatrix, typename TmatrixZ, typename Treal>
203 class EuclTruncationSymmWithZ : public EuclTruncationSymm<Tmatrix, Treal> {
204 public:
205 EuclTruncationSymmWithZ( Tmatrix & A_, TmatrixZ const & Z_ )
206 : EuclTruncationSymm<Tmatrix, Treal>(A_), Z(Z_) {}
207 protected:
208 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
209 Treal const threshold );
210 // getFrobSqNorms(...) from EuclTruncationSymm
211 // frobThreshLowLevel(...) from EuclTruncationSymm
212 virtual Interval<Treal> euclIfSmall( Treal const absTol,
213 Treal const threshold );
214 TmatrixZ const & Z;
215 }; // end class EuclTruncationSymmWithZ
216
217 template<typename Tmatrix, typename TmatrixZ, typename Treal>
219 Treal const threshold ) {
220 Treal Zfrob = Z.frob();
221 Treal thresholdTakingZIntoAccount = threshold / (Zfrob * Zfrob);
222 /* Divide by 2 because of symmetry */
223 lowTrunc = thresholdTakingZIntoAccount * thresholdTakingZIntoAccount / 2.0;
224 highTrunc = template_blas_get_num_limit_max<Treal>();
225 }
226
227 template<typename Tmatrix, typename TmatrixZ, typename Treal>
229 Treal const threshold ) {
230 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
231 mat::TripleMatrix<Tmatrix, TmatrixZ, Treal> ErrMatTriple( this->E, Z);
232 Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);
233 if ( tmpInterval.length() < 2*absTol )
234 return Interval<Treal>( tmpInterval.midPoint()-absTol,
235 tmpInterval.midPoint()+absTol );
236 return tmpInterval;
237 }
238
245 template<typename Tmatrix, typename Treal>
246 class EuclTruncationSymmElementLevel : public EuclTruncationSymm<Tmatrix, Treal> {
247 public:
248 explicit EuclTruncationSymmElementLevel( Tmatrix & A_ )
249 : EuclTruncationSymm<Tmatrix, Treal>(A_) {}
250 protected:
251 // getFrobTruncBounds(...) from EuclTruncationSymm
252 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
253 virtual void frobThreshLowLevel( Treal const threshold );
254 // Interval<Treal> euclIfSmall(...) from EuclTruncationSymm
255 }; // end class EuclTruncationSymmElementLevel
256
257 template<typename Tmatrix, typename Treal>
258 void EuclTruncationSymmElementLevel<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
259 this->A.getMatrix().getFrobSqElementLevel(frobsq_norms);
260 }
261
262 template<typename Tmatrix, typename Treal>
264 this->A.getMatrix().frobThreshElementLevel(threshold, &this->E.getMatrix() );
265 }
266
271 template<typename Tmatrix, typename Treal>
272 class EuclTruncationGeneral : public EuclTruncationBase<Tmatrix, Treal> {
273 public:
274 explicit EuclTruncationGeneral( Tmatrix & A_ )
275 : EuclTruncationBase<Tmatrix, Treal>(A_) {}
276 protected:
277 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
278 Treal const threshold );
279 virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
280 virtual void frobThreshLowLevel( Treal const threshold );
281 virtual Interval<Treal> euclIfSmall( Treal const absTol,
282 Treal const threshold );
283 }; // end class EuclTruncationGeneral
284
285 template<typename Tmatrix, typename Treal>
286 void EuclTruncationGeneral<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
287 Treal const threshold ) {
288 // Try to improve bounds based on the Frobenius norm
289 /* ||E||_F^2 <= thres^2 ->
290 * ||E||_F <= thres ->
291 * ||E||_2 <= thresh
292 */
293 lowTrunc = (threshold * threshold);
294 /* ||E||_F^2 >= thres^2 * n ->
295 * ||E||_F >= thres * sqrt(n) ->
296 * ||E||_2 >= thresh
297 */
298 highTrunc = (threshold * threshold * this->A.get_nrows());
299 }
300
301 template<typename Tmatrix, typename Treal>
302 void EuclTruncationGeneral<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
303 this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
304 }
305
306 template<typename Tmatrix, typename Treal>
308 this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
309 }
310
311 template<typename Tmatrix, typename Treal>
313 Treal const threshold ) {
314 // FIXME: this should be changed (for all trunc classes) so that
315 // some relative precision is always requested instead of the input
316 // absTol which in the current case is not used(!)
317 mat::ATAMatrix<Tmatrix, Treal> EtE(this->E);
318 Treal absTolDummy = template_blas_get_num_limit_max<Treal>(); // Treal(threshold * 1e-2)
319 Treal relTol = 100 * mat::getMachineEpsilon<Treal>();
320 Interval<Treal> tmpInterval = mat::euclIfSmall(EtE, absTolDummy, relTol, threshold);
321 tmpInterval = Interval<Treal>( template_blas_sqrt(tmpInterval.low()), template_blas_sqrt(tmpInterval.upp()) );
322 if ( tmpInterval.length() < 2*absTol )
323 return Interval<Treal>( tmpInterval.midPoint()-absTol,
324 tmpInterval.midPoint()+absTol );
325 return tmpInterval;
326 }
327
328
329
330
337 template<typename Tmatrix, typename TmatrixB, typename Treal>
339 public:
340 EuclTruncationCongrTransMeasure( Tmatrix & A_, TmatrixB const & B_ )
341 : EuclTruncationGeneral<Tmatrix, Treal>(A_), B(B_) {}
342 protected:
343 virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc,
344 Treal const threshold );
345 // getFrobSqNorms(...) from EuclTruncationGeneral
346 // frobThreshLowLevel(...) from EuclTruncationGeneral
347 virtual Interval<Treal> euclIfSmall( Treal const absTol,
348 Treal const threshold );
349 TmatrixB const & B;
350 }; // end class EuclTruncationCongrTransMeasure
351
352 template<typename Tmatrix, typename TmatrixB, typename Treal>
354 Treal & highTrunc,
355 Treal const threshold ) {
356 Treal Afrob = this->A.frob();
357 Treal Bfrob = B.frob();
358 Treal tmp = -Afrob + template_blas_sqrt( Afrob*Afrob + threshold / Bfrob );
359 lowTrunc = tmp*tmp;
360 highTrunc = template_blas_get_num_limit_max<Treal>();
361 }
362
363 template<typename Tmatrix, typename TmatrixB, typename Treal>
365 Treal const threshold ) {
366 Treal relTol = template_blas_sqrt(template_blas_sqrt(mat::getMachineEpsilon<Treal>()));
367 mat::CongrTransErrorMatrix<TmatrixB, Tmatrix, Treal> ErrMatTriple( B, this->A, this->E );
368 Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);
369 if ( tmpInterval.length() < 2*absTol ) {
370 return Interval<Treal>( tmpInterval.midPoint()-absTol,
371 tmpInterval.midPoint()+absTol );
372 }
373 return tmpInterval;
374 }
375
376
377} /* end namespace mat */
378#endif
Definition truncation.h:53
Treal run(Treal const threshold)
Definition truncation.h:80
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)=0
virtual ~EuclTruncationBase()
Definition truncation.h:57
Tmatrix E
Definition truncation.h:66
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)=0
virtual void frobThreshLowLevel(Treal const threshold)=0
EuclTruncationBase(Tmatrix &A_)
Definition truncation.h:70
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)=0
Tmatrix & A
Definition truncation.h:65
Truncation of general matrices with impact on matrix triple multiply as error measure.
Definition truncation.h:338
EuclTruncationCongrTransMeasure(Tmatrix &A_, TmatrixB const &B_)
Definition truncation.h:340
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition truncation.h:353
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition truncation.h:364
TmatrixB const & B
Definition truncation.h:349
Truncation of general matrices.
Definition truncation.h:272
virtual void frobThreshLowLevel(Treal const threshold)
Definition truncation.h:307
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition truncation.h:286
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition truncation.h:302
EuclTruncationGeneral(Tmatrix &A_)
Definition truncation.h:274
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition truncation.h:312
Truncation of symmetric matrices at the element level (used for mixed norm truncation)
Definition truncation.h:246
EuclTruncationSymmElementLevel(Tmatrix &A_)
Definition truncation.h:248
virtual void frobThreshLowLevel(Treal const threshold)
Definition truncation.h:263
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition truncation.h:258
Truncation of symmetric matrices with Z.
Definition truncation.h:203
TmatrixZ const & Z
Definition truncation.h:214
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition truncation.h:228
EuclTruncationSymmWithZ(Tmatrix &A_, TmatrixZ const &Z_)
Definition truncation.h:205
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition truncation.h:218
Truncation of symmetric matrices.
Definition truncation.h:154
virtual void frobThreshLowLevel(Treal const threshold)
Definition truncation.h:181
EuclTruncationSymm(Tmatrix &A_)
Definition truncation.h:156
virtual void getFrobTruncBounds(Treal &lowTrunc, Treal &highTrunc, Treal const threshold)
Definition truncation.h:168
virtual void getFrobSqNorms(std::vector< Treal > &frobsq_norms)
Definition truncation.h:176
virtual Interval< Treal > euclIfSmall(Treal const absTol, Treal const threshold)
Definition truncation.h:186
Definition Interval.h:46
Treal low() const
Definition Interval.h:144
Treal midPoint() const
Definition Interval.h:115
Treal length() const
Returns the length of the interval.
Definition Interval.h:109
Treal upp() const
Definition Interval.h:145
Describes dimensions of matrix and its blocks on all levels.
Definition SizesAndBlocks.h:45
mat::SizesAndBlocks rows
Definition test.cc:51
mat::SizesAndBlocks cols
Definition test.cc:52
#define B
#define A
Definition allocate.cc:39
Interval< Treal > euclIfSmall(Tmatrix const &M, Treal const requestedAbsAccuracy, Treal const requestedRelAccuracy, Treal const maxAbsVal, typename Tmatrix::VectorType *const eVecPtr=0, int maxIter=-1)
Returns interval containing the Euclidean norm of the matrix M.
Definition LanczosLargestMagnitudeEig.h:260
Definition mat_utils.h:95
Definition mat_utils.h:154
Definition mat_utils.h:123
Treal template_blas_sqrt(Treal x)