ergo
mat_acc_extrapolate.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 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 
28 #ifndef ERGO_MAT_ACC_EXTRAPOLATE_HEADER
29 #define ERGO_MAT_ACC_EXTRAPOLATE_HEADER
30 
31 
32 #include <vector>
33 
34 
35 #include "matrix_utilities.h"
36 
37 
38 
39 template<class Treal, class Tworker>
41 {
42  public:
43  explicit MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_);
44  void Scan(const Tworker & worker,
45  Treal firstParam,
46  Treal stepFactor,
47  int nSteps);
48  void GetScanResult(Treal* threshList_,
49  Treal* errorList_frob_,
50  Treal* errorList_eucl_,
51  Treal* errorList_maxe_,
52  Treal* timeList_);
53  private:
56  Treal baseThresh;
57  std::vector<Treal> threshList;
58  std::vector<Treal> errorList_frob; // Frobenius norm
59  std::vector<Treal> errorList_eucl; // Euclidean norm
60  std::vector<Treal> errorList_maxe; // Max element norm
61  std::vector<Treal> timeList;
62 };
63 
64 
65 template<class Treal, class Tworker>
67  : matrix_size_block_info(matrix_size_block_info_)
68 {}
69 
70 
71 template<class Treal, class Tworker>
73  ::Scan(const Tworker & worker,
74  Treal firstParam,
75  Treal stepFactor,
76  int nSteps)
77 {
78  nScanSteps = nSteps;
79  baseThresh = firstParam;
80  threshList.resize(nSteps);
81  errorList_frob.resize(nSteps);
82  errorList_eucl.resize(nSteps);
83  errorList_maxe.resize(nSteps);
84  timeList.resize(nSteps);
85 
86  // Prepare matrix objects
87  symmMatrix accurateMatrix;
88  accurateMatrix.resetSizesAndBlocks(matrix_size_block_info,
89  matrix_size_block_info);
90  symmMatrix otherMatrix;
91  otherMatrix.resetSizesAndBlocks(matrix_size_block_info,
92  matrix_size_block_info);
93  symmMatrix errorMatrix;
94  errorMatrix.resetSizesAndBlocks(matrix_size_block_info,
95  matrix_size_block_info);
96 
97  // Compute "accurate" matrix
98  worker.ComputeMatrix(firstParam, accurateMatrix);
99  // Compute other matrices and compare them to "accurate" matrix
100  Treal currParam = firstParam;
101  for(int i = 0; i < nSteps; i++)
102  {
103  currParam *= stepFactor;
104  time_t startTime, endTime;
105  time(&startTime);
106  worker.ComputeMatrix(currParam, otherMatrix);
107  time(&endTime);
108  timeList[i] = endTime - startTime;
109  threshList[i] = currParam;
110  // Compute error matrix
111  errorMatrix = otherMatrix;
112  errorMatrix += (ergo_real)(-1) * accurateMatrix;
113  // Compute different norms of error matrix
114  // Frobenius norm
115  errorList_frob[i] = errorMatrix.frob();
116  // Euclidean norm
117  Treal euclAcc = 1e-11;
118  errorList_eucl[i] = errorMatrix.eucl(euclAcc);
119  // Max element norm
120  errorList_maxe[i] = compute_maxabs_sparse(errorMatrix);
121  }
122 
123 }
124 
125 
126 template<class Treal, class Tworker>
128  ::GetScanResult(Treal* threshList_,
129  Treal* errorList_frob_,
130  Treal* errorList_eucl_,
131  Treal* errorList_maxe_,
132  Treal* timeList_)
133 {
134  for(int i = 0; i < nScanSteps; i++)
135  {
136  threshList_[i] = threshList[i];
137  errorList_frob_[i] = errorList_frob[i];
138  errorList_eucl_[i] = errorList_eucl[i];
139  errorList_maxe_[i] = errorList_maxe[i];
140  timeList_ [i] = timeList [i];
141  }
142 }
143 
144 
145 
146 
147 #endif