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>
39 namespace 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
template_blas_sqrt
Treal template_blas_sqrt(Treal x)
mat::maxdiff
static Treal maxdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:44
mat::read_matrix
static void read_matrix(Tmatrix &A, char const *const matrixPath, int const size)
Definition: general.h:115
mat::frobdiff
static Treal frobdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:69
template_blas_fabs
Treal template_blas_fabs(Treal x)
mat::fileread
static void fileread(Treal *ptr, int size, FILE *file)
Definition: general.h:100
mat::Failure
Definition: Failure.h:57
mat
Definition: allocate.cc:39
A
#define A
mat::read_xyz
static void read_xyz(Treal *x, Treal *y, Treal *z, char *atomsPath, int const natoms, int const size)
Definition: general.h:175
mat::maxdiff_tri
static Treal maxdiff_tri(const Treal *f1, const Treal *f2, int size)
Definition: general.h:56
mat::read_sparse_matrix
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