ergo
mat_acc_extrapolate.h
Go to the documentation of this file.
1/* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 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
51template<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:
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
77template<class Treal, class Tworker>
79 : matrix_size_block_info(matrix_size_block_info_)
80{}
81
82
83template<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
138template<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
Definition: mat_acc_extrapolate.h:53
void Scan(const Tworker &worker, Treal firstParam, Treal stepFactor, int nSteps)
Definition: mat_acc_extrapolate.h:85
Treal baseThresh
Definition: mat_acc_extrapolate.h:68
std::vector< Treal > errorList_frob
Definition: mat_acc_extrapolate.h:70
std::vector< Treal > threshList
Definition: mat_acc_extrapolate.h:69
std::vector< Treal > errorList_eucl
Definition: mat_acc_extrapolate.h:71
mat::SizesAndBlocks matrix_size_block_info
Definition: mat_acc_extrapolate.h:66
std::vector< Treal > timeList
Definition: mat_acc_extrapolate.h:73
int nScanSteps
Definition: mat_acc_extrapolate.h:67
MatAccInvestigator(mat::SizesAndBlocks const &matrix_size_block_info_)
Definition: mat_acc_extrapolate.h:78
std::vector< Treal > errorList_maxe
Definition: mat_acc_extrapolate.h:72
void GetScanResult(Treal *threshList_, Treal *errorList_frob_, Treal *errorList_eucl_, Treal *errorList_maxe_, Treal *timeList_)
Definition: mat_acc_extrapolate.h:140
void resetSizesAndBlocks(SizesAndBlocks const &newRows, SizesAndBlocks const &newCols)
Definition: MatrixBase.h:76
Symmetric matrix.
Definition: MatrixSymmetric.h:68
Treal eucl(Treal const requestedAccuracy, int maxIter=-1) const
Definition: MatrixSymmetric.h:673
Treal frob() const
Definition: MatrixSymmetric.h:360
Describes dimensions of matrix and its blocks on all levels.
Definition: SizesAndBlocks.h:45
Utilities related to the hierarchical matrix library (HML), including functions for setting up permut...
ergo_real compute_maxabs_sparse(const Tmatrix &M)
Definition: matrix_utilities.h:97
double ergo_real
Definition: realtype.h:69