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 30 /** @file mat_acc_extrapolate.h 31 32 @brief Functionality for performing a scan testing different 33 threshold values for some kind of matrix computation. Can be used 34 for accuracy scans of e.g. Coulomb of HF exchange matrix 35 evaluation. 36 37 @author: Elias Rudberg <em>responsible</em> 38 */ 39 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> 52 class MatAccInvestigator 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: 66 mat::SizesAndBlocks matrix_size_block_info; 67 int nScanSteps; 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> MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_)78 MatAccInvestigator<Treal, Tworker>::MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_) 79 : matrix_size_block_info(matrix_size_block_info_) 80 {} 81 82 83 template<class Treal, class Tworker> 84 void MatAccInvestigator<Treal, Tworker> Scan(const Tworker & worker,Treal firstParam,Treal stepFactor,int nSteps)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> 139 void MatAccInvestigator<Treal, Tworker> GetScanResult(Treal * threshList_,Treal * errorList_frob_,Treal * errorList_eucl_,Treal * errorList_maxe_,Treal * timeList_)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 160