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 #if !defined(_SPARSE_MATRIX_H_) 00029 #define _SPARSE_MATRIX_H_ 1 00030 00036 #include <stdio.h> 00037 #include <vector> 00038 #include <algorithm> 00039 00040 00041 #include "realtype.h" 00042 #include "matrix_typedefs.h" 00043 #include "basisinfo.h" 00044 #include "sparse_pattern.h" 00045 00046 #if !defined(BEGIN_NAMESPACE) 00047 #define BEGIN_NAMESPACE(x) namespace x { 00048 #define END_NAMESPACE(x) } /* x */ 00049 #endif 00050 00051 BEGIN_NAMESPACE(Dft) 00052 00053 00054 class SparseMatrix { 00055 class Exception : public std::exception { 00056 const char *msg; 00057 public: 00058 explicit Exception(const char *msg_) : msg(msg_) {} 00059 virtual const char *what() const throw() { return msg; } 00060 }; 00061 00062 const SparsePattern& pattern; 00063 ergo_real **columns; 00064 int **offsets; 00065 int **his; 00066 int *cnt; 00067 int n; 00069 void createOffsets(const SparsePattern& pattern); 00070 public: 00073 explicit SparseMatrix(const SparsePattern& pattern_); 00074 SparseMatrix(const SparsePattern& pattern_, 00075 const symmMatrix& m, const int *aoMap, 00076 std::vector<int> const & permutationHML); 00077 00078 ~SparseMatrix() { 00079 for(int i=0; i<n; i++) { 00080 delete [](columns[i]); 00081 delete [](offsets[i]); 00082 delete [](his[i]); 00083 } 00084 delete []columns; 00085 delete []offsets; 00086 delete []his; 00087 delete []cnt; 00088 } 00089 00090 void print(const char *title) const; 00091 00093 void addSymmetrizedTo(symmMatrix& sMat, 00094 const int *aoMap, 00095 std::vector<int> const & permutationHML) const; 00096 00100 void add(int row, int col, ergo_real val) { 00101 ergo_real *columnData = columns[col]; 00102 const int *hi = his[col]; 00103 int idx; 00104 for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx); 00105 //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi; 00106 if(idx >= cnt[col]) 00107 throw Exception("SparseMatrix::add called with incorrect args"); 00108 int offset = offsets[col][idx]; 00109 /* Add it... */ 00110 //assert(row-offset>=0); 00111 //assert(row-offset<pattern.getColumnSize(col)); 00112 columnData[row-offset] += val; 00113 } 00114 00115 /* This operator[] syntax glue that we could in principle use can be 00116 expensive performance-wise, do it the old-fashioned way. 00117 Checking against intervals.end() is *terribly* expensive!!! 00118 */ 00119 ergo_real at(int row, int col) const { 00120 const ergo_real *columnData = columns[col]; 00121 const int *hi = his[col]; 00122 int idx; for(idx = 0; idx < cnt[col] && row >hi[idx]; ++idx); 00123 if(idx >= cnt[col]) 00124 throw Exception("SparseMatrix::at called with incorrect args"); 00125 //int idx = std::upper_bound(hi, hi+cnt[col], row)-hi; 00126 int offset = offsets[col][idx]; 00127 /* return it... */ 00128 //assert(row-offset>=0); 00129 //assert(row-offset<pattern.getColumnSize(col)); 00130 return columnData[row-offset]; 00131 } 00132 }; 00133 00134 00135 END_NAMESPACE(Dft) 00136 00137 void 00138 getrho_blocked_lda(int nbast, const Dft::SparseMatrix& dmat, 00139 const ergo_real * gao, 00140 const int* nblocks, const int (*iblocks)[2], 00141 int ldaib, ergo_real *tmp, int nvclen, ergo_real *rho); 00142 void 00143 getrho_blocked_gga(int nbast, const Dft::SparseMatrix& dmat, 00144 const ergo_real * gao, 00145 const int* nblocks, const int (*iblocks)[2], 00146 int ldaib, ergo_real *tmp, int nvclen, 00147 ergo_real *rho, ergo_real (*grad)[3]); 00148 00149 #endif /* _SPARSE_MATRIX_H_ */