mat_acc_extrapolate.h

Go to the documentation of this file.
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

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1