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 general.h 31 32 \brief Some general utilities used by other parts of the 33 hierarchical matrix library. 34 */ 35 36 #ifndef MAT_GENERAL 37 #define MAT_GENERAL 38 #include <cassert> 39 namespace mat { 40 41 42 43 template<class Treal> maxdiff(const Treal * f1,const Treal * f2,int size)44 static Treal maxdiff(const Treal* f1,const Treal* f2,int size) { 45 Treal diff = 0; 46 Treal tmpdiff; 47 for(int i = 0; i < size * size; i++) { 48 tmpdiff = template_blas_fabs(f1[i] - f2[i]); 49 if (tmpdiff > 0) 50 diff = (diff > tmpdiff ? diff : tmpdiff); 51 } 52 return diff; 53 } 54 55 template<class Treal> maxdiff_tri(const Treal * f1,const Treal * f2,int size)56 static Treal maxdiff_tri(const Treal* f1,const Treal* f2,int size) { 57 Treal diff = 0; 58 Treal tmpdiff; 59 for (int col = 0; col < size; col++) 60 for (int row = 0; row < col + 1; row++) { 61 tmpdiff = template_blas_fabs(f1[col * size + row] - f2[col * size + row]); 62 diff = (diff > tmpdiff ? diff : tmpdiff); 63 } 64 return diff; 65 } 66 67 68 template<class Treal> frobdiff(const Treal * f1,const Treal * f2,int size)69 static Treal frobdiff(const Treal* f1,const Treal* f2,int size) { 70 Treal diff = 0; 71 Treal tmp; 72 for(int i = 0; i < size * size; i++) { 73 tmp = f1[i] - f2[i]; 74 diff += tmp * tmp; 75 } 76 return template_blas_sqrt(diff); 77 } 78 79 #if 0 80 template<class T> 81 static void fileread(T *ptr,int size,FILE*) { 82 std::cout<<"error reading file"<<std::endl; 83 } 84 template<> 85 void fileread<double>(double *ptr,int size,FILE* file) { 86 fread(ptr,sizeof(double),size*size,file); 87 } 88 template<> 89 void fileread<float>(float *ptr,int size,FILE* file) { 90 double* tmpptr=new double [size*size]; 91 fread(tmpptr,sizeof(double),size*size,file); 92 for (int i=0;i<size*size;i++) 93 { 94 ptr[i]=(float)tmpptr[i]; 95 } 96 delete[] tmpptr; 97 } 98 #else 99 template<typename Treal, typename Trealonfile> fileread(Treal * ptr,int size,FILE * file)100 static void fileread(Treal *ptr, int size, FILE* file) { 101 if (sizeof(Trealonfile) == sizeof(Treal)) 102 fread(ptr,sizeof(Treal),size,file); 103 else { 104 Trealonfile* tmpptr=new Trealonfile[size]; 105 fread(tmpptr,sizeof(Trealonfile),size,file); 106 for (int i = 0; i < size; i++) { 107 ptr[i]=(Treal)tmpptr[i]; 108 } 109 delete[] tmpptr; 110 } 111 } 112 #endif 113 114 template<typename Treal, typename Tmatrix> read_matrix(Tmatrix & A,char const * const matrixPath,int const size)115 static void read_matrix(Tmatrix& A, 116 char const * const matrixPath, 117 int const size) { 118 FILE* matrixfile=fopen(matrixPath,"rb"); 119 if (!matrixfile) { 120 throw Failure("read_matrix: Cannot open inputfile"); 121 } 122 Treal* matrixfull = new Treal [size*size]; 123 fileread<Treal, double>(matrixfull, size*size, matrixfile); 124 /* A must already have built data structure */ 125 A.assign_from_full(matrixfull, size, size); 126 delete[] matrixfull; 127 return; 128 } 129 130 template<typename Treal, typename Trealonfile, typename Tmatrix> read_sparse_matrix(Tmatrix & A,char const * const rowPath,char const * const colPath,char const * const valPath,int const nval)131 static void read_sparse_matrix(Tmatrix& A, 132 char const * const rowPath, 133 char const * const colPath, 134 char const * const valPath, 135 int const nval) { 136 FILE* rowfile=fopen(rowPath,"rb"); 137 if (!rowfile) { 138 throw Failure("read_matrix: Cannot open inputfile rowfile"); 139 } 140 FILE* colfile=fopen(colPath,"rb"); 141 if (!colfile) { 142 throw Failure("read_matrix: Cannot open inputfile colfile"); 143 } 144 FILE* valfile=fopen(valPath,"rb"); 145 if (!valfile) { 146 throw Failure("read_matrix: Cannot open inputfile valfile"); 147 } 148 int* row = new int[nval]; 149 int* col = new int[nval]; 150 Treal* val = new Treal[nval]; 151 fileread<int, int>(row, nval, rowfile); 152 fileread<int, int>(col, nval, colfile); 153 fileread<Treal, Trealonfile>(val, nval, valfile); 154 155 /* A must already have built data structure */ 156 A.assign_from_sparse(row, col, val, nval); 157 #if 0 158 Treal* compval = new Treal[nval]; 159 A.get_values(row, col, compval, nval); 160 Treal maxdiff = 0; 161 Treal diff; 162 for (int i = 0; i < nval; i++) { 163 diff = template_blas_fabs(compval[i] - val[i]); 164 maxdiff = diff > maxdiff ? diff : maxdiff; 165 } 166 std::cout<<"Maxdiff: "<<maxdiff<<std::endl; 167 #endif 168 delete[] row; 169 delete[] col; 170 delete[] val; 171 return; 172 } 173 174 template<typename Treal> read_xyz(Treal * x,Treal * y,Treal * z,char * atomsPath,int const natoms,int const size)175 static void read_xyz(Treal* x, Treal* y, Treal* z, 176 char * atomsPath, 177 int const natoms, 178 int const size) { 179 char* atomfile(atomsPath); 180 std::ifstream input(atomfile); 181 if (!input) { 182 throw Failure("read_xyz: Cannot open inputfile"); 183 } 184 input >> std::setprecision(10); 185 Treal* xtmp = new Treal[natoms]; 186 Treal* ytmp = new Treal[natoms]; 187 Treal* ztmp = new Treal[natoms]; 188 int* atomstart = new int[natoms+1]; 189 for(int i = 0 ; i < natoms ; i++) { 190 input >> x[i]; 191 input >> y[i]; 192 input >> z[i]; 193 input >> atomstart[i]; 194 } 195 atomstart[natoms] = size; 196 for (int atom = 0; atom < natoms; atom++) 197 for (int bf = atomstart[atom]; bf < atomstart[atom + 1]; bf++) { 198 x[bf] = x[atom]; 199 y[bf] = y[atom]; 200 z[bf] = z[atom]; 201 } 202 delete[] xtmp; 203 delete[] ytmp; 204 delete[] ztmp; 205 delete[] atomstart; 206 } 207 } /* end namespace mat */ 208 209 #endif 210