ergo
general.h
Go to the documentation of this file.
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
36#ifndef MAT_GENERAL
37#define MAT_GENERAL
38#include <cassert>
39namespace mat {
40
41
42
43 template<class Treal>
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>
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>
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>
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>
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>
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>
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
Definition: Failure.h:57
#define A
Definition: allocate.cc:39
static Treal maxdiff_tri(const Treal *f1, const Treal *f2, int size)
Definition: general.h:56
static void read_xyz(Treal *x, Treal *y, Treal *z, char *atomsPath, int const natoms, int const size)
Definition: general.h:175
static Treal maxdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:44
static void fileread(Treal *ptr, int size, FILE *file)
Definition: general.h:100
static void read_sparse_matrix(Tmatrix &A, char const *const rowPath, char const *const colPath, char const *const valPath, int const nval)
Definition: general.h:131
static Treal frobdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:69
static void read_matrix(Tmatrix &A, char const *const matrixPath, int const size)
Definition: general.h:115
Treal template_blas_sqrt(Treal x)
Treal template_blas_fabs(Treal x)