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