ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 #ifndef ERGO_MAT_ACC_EXTRAPOLATE_HEADER 00029 #define ERGO_MAT_ACC_EXTRAPOLATE_HEADER 00030 00031 00032 #include <vector> 00033 00034 00035 #include "matrix_utilities.h" 00036 00037 00038 00039 template<class Treal, class Tworker> 00040 class MatAccInvestigator 00041 { 00042 public: 00043 explicit MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_); 00044 void Scan(const Tworker & worker, 00045 Treal firstParam, 00046 Treal stepFactor, 00047 int nSteps); 00048 void GetScanResult(Treal* threshList_, 00049 Treal* errorList_frob_, 00050 Treal* errorList_eucl_, 00051 Treal* errorList_maxe_, 00052 Treal* timeList_); 00053 private: 00054 mat::SizesAndBlocks matrix_size_block_info; 00055 int nScanSteps; 00056 Treal baseThresh; 00057 std::vector<Treal> threshList; 00058 std::vector<Treal> errorList_frob; // Frobenius norm 00059 std::vector<Treal> errorList_eucl; // Euclidean norm 00060 std::vector<Treal> errorList_maxe; // Max element norm 00061 std::vector<Treal> timeList; 00062 }; 00063 00064 00065 template<class Treal, class Tworker> 00066 MatAccInvestigator<Treal, Tworker>::MatAccInvestigator(mat::SizesAndBlocks const & matrix_size_block_info_) 00067 : matrix_size_block_info(matrix_size_block_info_) 00068 {} 00069 00070 00071 template<class Treal, class Tworker> 00072 void MatAccInvestigator<Treal, Tworker> 00073 ::Scan(const Tworker & worker, 00074 Treal firstParam, 00075 Treal stepFactor, 00076 int nSteps) 00077 { 00078 nScanSteps = nSteps; 00079 baseThresh = firstParam; 00080 threshList.resize(nSteps); 00081 errorList_frob.resize(nSteps); 00082 errorList_eucl.resize(nSteps); 00083 errorList_maxe.resize(nSteps); 00084 timeList.resize(nSteps); 00085 00086 // Prepare matrix objects 00087 symmMatrix accurateMatrix; 00088 accurateMatrix.resetSizesAndBlocks(matrix_size_block_info, 00089 matrix_size_block_info); 00090 symmMatrix otherMatrix; 00091 otherMatrix.resetSizesAndBlocks(matrix_size_block_info, 00092 matrix_size_block_info); 00093 symmMatrix errorMatrix; 00094 errorMatrix.resetSizesAndBlocks(matrix_size_block_info, 00095 matrix_size_block_info); 00096 00097 // Compute "accurate" matrix 00098 worker.ComputeMatrix(firstParam, accurateMatrix); 00099 // Compute other matrices and compare them to "accurate" matrix 00100 Treal currParam = firstParam; 00101 for(int i = 0; i < nSteps; i++) 00102 { 00103 currParam *= stepFactor; 00104 time_t startTime, endTime; 00105 time(&startTime); 00106 worker.ComputeMatrix(currParam, otherMatrix); 00107 time(&endTime); 00108 timeList[i] = endTime - startTime; 00109 threshList[i] = currParam; 00110 // Compute error matrix 00111 errorMatrix = otherMatrix; 00112 errorMatrix += (ergo_real)(-1) * accurateMatrix; 00113 // Compute different norms of error matrix 00114 // Frobenius norm 00115 errorList_frob[i] = errorMatrix.frob(); 00116 // Euclidean norm 00117 Treal euclAcc = 1e-11; 00118 errorList_eucl[i] = errorMatrix.eucl(euclAcc); 00119 // Max element norm 00120 errorList_maxe[i] = compute_maxabs_sparse(errorMatrix); 00121 } 00122 00123 } 00124 00125 00126 template<class Treal, class Tworker> 00127 void MatAccInvestigator<Treal, Tworker> 00128 ::GetScanResult(Treal* threshList_, 00129 Treal* errorList_frob_, 00130 Treal* errorList_eucl_, 00131 Treal* errorList_maxe_, 00132 Treal* timeList_) 00133 { 00134 for(int i = 0; i < nScanSteps; i++) 00135 { 00136 threshList_[i] = threshList[i]; 00137 errorList_frob_[i] = errorList_frob[i]; 00138 errorList_eucl_[i] = errorList_eucl[i]; 00139 errorList_maxe_[i] = errorList_maxe[i]; 00140 timeList_ [i] = timeList [i]; 00141 } 00142 } 00143 00144 00145 00146 00147 #endif