ergo
mat_acc_extrapolate.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
40 #ifndef ERGO_MAT_ACC_EXTRAPOLATE_HEADER
41 #define ERGO_MAT_ACC_EXTRAPOLATE_HEADER
42 
43 
44 #include <vector>
45 
46 
47 #include "matrix_utilities.h"
48 
49 
50 
51 template<class Treal, class Tworker>
53 {
54  public:
55  explicit MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_);
56  void Scan(const Tworker & worker,
57  Treal firstParam,
58  Treal stepFactor,
59  int nSteps);
60  void GetScanResult(Treal* threshList_,
61  Treal* errorList_frob_,
62  Treal* errorList_eucl_,
63  Treal* errorList_maxe_,
64  Treal* timeList_);
65  private:
68  Treal baseThresh;
69  std::vector<Treal> threshList;
70  std::vector<Treal> errorList_frob; // Frobenius norm
71  std::vector<Treal> errorList_eucl; // Euclidean norm
72  std::vector<Treal> errorList_maxe; // Max element norm
73  std::vector<Treal> timeList;
74 };
75 
76 
77 template<class Treal, class Tworker>
79  : matrix_size_block_info(matrix_size_block_info_)
80 {}
81 
82 
83 template<class Treal, class Tworker>
85  ::Scan(const Tworker & worker,
86  Treal firstParam,
87  Treal stepFactor,
88  int nSteps)
89 {
90  nScanSteps = nSteps;
91  baseThresh = firstParam;
92  threshList.resize(nSteps);
93  errorList_frob.resize(nSteps);
94  errorList_eucl.resize(nSteps);
95  errorList_maxe.resize(nSteps);
96  timeList.resize(nSteps);
97 
98  // Prepare matrix objects
99  symmMatrix accurateMatrix;
100  accurateMatrix.resetSizesAndBlocks(matrix_size_block_info,
101  matrix_size_block_info);
102  symmMatrix otherMatrix;
103  otherMatrix.resetSizesAndBlocks(matrix_size_block_info,
104  matrix_size_block_info);
105  symmMatrix errorMatrix;
106  errorMatrix.resetSizesAndBlocks(matrix_size_block_info,
107  matrix_size_block_info);
108 
109  // Compute "accurate" matrix
110  worker.ComputeMatrix(firstParam, accurateMatrix);
111  // Compute other matrices and compare them to "accurate" matrix
112  Treal currParam = firstParam;
113  for(int i = 0; i < nSteps; i++)
114  {
115  currParam *= stepFactor;
116  time_t startTime, endTime;
117  time(&startTime);
118  worker.ComputeMatrix(currParam, otherMatrix);
119  time(&endTime);
120  timeList[i] = endTime - startTime;
121  threshList[i] = currParam;
122  // Compute error matrix
123  errorMatrix = otherMatrix;
124  errorMatrix += (ergo_real)(-1) * accurateMatrix;
125  // Compute different norms of error matrix
126  // Frobenius norm
127  errorList_frob[i] = errorMatrix.frob();
128  // Euclidean norm
129  Treal euclAcc = 1e-11;
130  errorList_eucl[i] = errorMatrix.eucl(euclAcc);
131  // Max element norm
132  errorList_maxe[i] = compute_maxabs_sparse(errorMatrix);
133  }
134 
135 }
136 
137 
138 template<class Treal, class Tworker>
140  ::GetScanResult(Treal* threshList_,
141  Treal* errorList_frob_,
142  Treal* errorList_eucl_,
143  Treal* errorList_maxe_,
144  Treal* timeList_)
145 {
146  for(int i = 0; i < nScanSteps; i++)
147  {
148  threshList_[i] = threshList[i];
149  errorList_frob_[i] = errorList_frob[i];
150  errorList_eucl_[i] = errorList_eucl[i];
151  errorList_maxe_[i] = errorList_maxe[i];
152  timeList_ [i] = timeList [i];
153  }
154 }
155 
156 
157 
158 
159 #endif
double ergo_real
Definition: realtype.h:69
std::vector< Treal > timeList
Definition: mat_acc_extrapolate.h:73
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
Definition: mat_acc_extrapolate.h:52
ergo_real compute_maxabs_sparse(const Tmatrix &M)
Definition: matrix_utilities.h:97
MatAccInvestigator(mat::SizesAndBlocks const &matrix_size_block_info_)
Definition: mat_acc_extrapolate.h:78
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
std::vector< Treal > errorList_eucl
Definition: mat_acc_extrapolate.h:71
std::vector< Treal > threshList
Definition: mat_acc_extrapolate.h:69
std::vector< Treal > errorList_frob
Definition: mat_acc_extrapolate.h:70
void Scan(const Tworker &worker, Treal firstParam, Treal stepFactor, int nSteps)
Definition: mat_acc_extrapolate.h:85
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:667
Treal frob() const
Definition: MatrixSymmetric.h:354
int nScanSteps
Definition: mat_acc_extrapolate.h:67
Symmetric matrix.
Definition: MatrixBase.h:51
void GetScanResult(Treal *threshList_, Treal *errorList_frob_, Treal *errorList_eucl_, Treal *errorList_maxe_, Treal *timeList_)
Definition: mat_acc_extrapolate.h:140
mat::SizesAndBlocks matrix_size_block_info
Definition: mat_acc_extrapolate.h:66
Treal baseThresh
Definition: mat_acc_extrapolate.h:68
std::vector< Treal > errorList_maxe
Definition: mat_acc_extrapolate.h:72