ergo
general.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28 #ifndef MAT_GENERAL
29 #define MAT_GENERAL
30 #include <cassert>
31 namespace mat {
32 
33 
34 
35  template<class Treal>
36  static Treal maxdiff(const Treal* f1,const Treal* f2,int size) {
37  Treal diff = 0;
38  Treal tmpdiff;
39  for(int i = 0; i < size * size; i++) {
40  tmpdiff = template_blas_fabs(f1[i] - f2[i]);
41  if (tmpdiff > 0)
42  diff = (diff > tmpdiff ? diff : tmpdiff);
43  }
44  return diff;
45  }
46 
47  template<class Treal>
48  static Treal maxdiff_tri(const Treal* f1,const Treal* f2,int size) {
49  Treal diff = 0;
50  Treal tmpdiff;
51  for (int col = 0; col < size; col++)
52  for (int row = 0; row < col + 1; row++) {
53  tmpdiff = template_blas_fabs(f1[col * size + row] - f2[col * size + row]);
54  diff = (diff > tmpdiff ? diff : tmpdiff);
55  }
56  return diff;
57  }
58 
59 
60  template<class Treal>
61  static Treal frobdiff(const Treal* f1,const Treal* f2,int size) {
62  Treal diff = 0;
63  Treal tmp;
64  for(int i = 0; i < size * size; i++) {
65  tmp = f1[i] - f2[i];
66  diff += tmp * tmp;
67  }
68  return template_blas_sqrt(diff);
69  }
70 
71 #if 0
72  template<class T>
73  static void fileread(T *ptr,int size,FILE*) {
74  std::cout<<"error reading file"<<std::endl;
75  }
76  template<>
77  void fileread<double>(double *ptr,int size,FILE* file) {
78  fread(ptr,sizeof(double),size*size,file);
79  }
80  template<>
81  void fileread<float>(float *ptr,int size,FILE* file) {
82  double* tmpptr=new double [size*size];
83  fread(tmpptr,sizeof(double),size*size,file);
84  for (int i=0;i<size*size;i++)
85  {
86  ptr[i]=(float)tmpptr[i];
87  }
88  delete[] tmpptr;
89  }
90 #else
91  template<typename Treal, typename Trealonfile>
92  static void fileread(Treal *ptr, int size, FILE* file) {
93  if (sizeof(Trealonfile) == sizeof(Treal))
94  fread(ptr,sizeof(Treal),size,file);
95  else {
96  Trealonfile* tmpptr=new Trealonfile[size];
97  fread(tmpptr,sizeof(Trealonfile),size,file);
98  for (int i = 0; i < size; i++) {
99  ptr[i]=(Treal)tmpptr[i];
100  }
101  delete[] tmpptr;
102  }
103  }
104 #endif
105 
106  template<typename Treal, typename Tmatrix>
107  static void read_matrix(Tmatrix& A,
108  char const * const matrixPath,
109  int const size) {
110  FILE* matrixfile=fopen(matrixPath,"rb");
111  if (!matrixfile) {
112  throw Failure("read_matrix: Cannot open inputfile");
113  }
114  Treal* matrixfull = new Treal [size*size];
115  fileread<Treal, double>(matrixfull, size*size, matrixfile);
116  /* A must already have built data structure */
117  A.assign_from_full(matrixfull, size, size);
118  delete[] matrixfull;
119  return;
120  }
121 
122  template<typename Treal, typename Trealonfile, typename Tmatrix>
123  static void read_sparse_matrix(Tmatrix& A,
124  char const * const rowPath,
125  char const * const colPath,
126  char const * const valPath,
127  int const nval) {
128  FILE* rowfile=fopen(rowPath,"rb");
129  if (!rowfile) {
130  throw Failure("read_matrix: Cannot open inputfile rowfile");
131  }
132  FILE* colfile=fopen(colPath,"rb");
133  if (!colfile) {
134  throw Failure("read_matrix: Cannot open inputfile colfile");
135  }
136  FILE* valfile=fopen(valPath,"rb");
137  if (!valfile) {
138  throw Failure("read_matrix: Cannot open inputfile valfile");
139  }
140  int* row = new int[nval];
141  int* col = new int[nval];
142  Treal* val = new Treal[nval];
143  fileread<int, int>(row, nval, rowfile);
144  fileread<int, int>(col, nval, colfile);
145  fileread<Treal, Trealonfile>(val, nval, valfile);
146 
147  /* A must already have built data structure */
148  A.assign_from_sparse(row, col, val, nval);
149 #if 0
150  Treal* compval = new Treal[nval];
151  A.get_values(row, col, compval, nval);
152  Treal maxdiff = 0;
153  Treal diff;
154  for (int i = 0; i < nval; i++) {
155  diff = template_blas_fabs(compval[i] - val[i]);
156  maxdiff = diff > maxdiff ? diff : maxdiff;
157  }
158  std::cout<<"Maxdiff: "<<maxdiff<<std::endl;
159 #endif
160  delete[] row;
161  delete[] col;
162  delete[] val;
163  return;
164  }
165 
166  template<typename Treal>
167  static void read_xyz(Treal* x, Treal* y, Treal* z,
168  char * atomsPath,
169  int const natoms,
170  int const size) {
171  char* atomfile(atomsPath);
172  std::ifstream input(atomfile);
173  if (!input) {
174  throw Failure("read_xyz: Cannot open inputfile");
175  }
176  input >> std::setprecision(10);
177  Treal* xtmp = new Treal[natoms];
178  Treal* ytmp = new Treal[natoms];
179  Treal* ztmp = new Treal[natoms];
180  int* atomstart = new int[natoms+1];
181  for(int i = 0 ; i < natoms ; i++) {
182  input >> x[i];
183  input >> y[i];
184  input >> z[i];
185  input >> atomstart[i];
186  }
187  atomstart[natoms] = size;
188  for (int atom = 0; atom < natoms; atom++)
189  for (int bf = atomstart[atom]; bf < atomstart[atom + 1]; bf++) {
190  x[bf] = x[atom];
191  y[bf] = y[atom];
192  z[bf] = z[atom];
193  }
194  delete[] xtmp;
195  delete[] ytmp;
196  delete[] ztmp;
197  delete[] atomstart;
198  }
199 } /* end namespace mat */
200 
201 #endif
#define A
Definition: allocate.cc:30
Treal template_blas_fabs(Treal x)
static Treal maxdiff_tri(const Treal *f1, const Treal *f2, int size)
Definition: general.h:48
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:123
static void read_xyz(Treal *x, Treal *y, Treal *z, char *atomsPath, int const natoms, int const size)
Definition: general.h:167
static void read_matrix(Tmatrix &A, char const *const matrixPath, int const size)
Definition: general.h:107
static Treal maxdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:36
static int input(void)
Definition: ergo_input_parser.c:1455
static Treal frobdiff(const Treal *f1, const Treal *f2, int size)
Definition: general.h:61
Definition: Failure.h:47
static void fileread(Treal *ptr, int size, FILE *file)
Definition: general.h:92
Treal template_blas_sqrt(Treal x)