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 grid_reader.cc
31
32 @brief Grid Generator interface. Functions for opening grid file,
33 reading chunks from it, and closing the file, are provided.
34
35 @author: Pawel Salek <em>responsible</em>
36 */
37
38 #include <stdio.h>
39 #include <stdlib.h>
40 #include <string.h>
41 #include <unistd.h>
42 #include <pthread.h>
43
44 #include "dft_common.h"
45 #include "grid_reader.h"
46 #include "grid_stream.h"
47 #include "grid_hicu.h"
48 #include "grid_matrix.h"
49
50 class FullMatrixWrapper : public Dft::Matrix {
51 const Dft::FullMatrix& matrix;
52 public:
FullMatrixWrapper(const Dft::FullMatrix & m)53 explicit FullMatrixWrapper(const Dft::FullMatrix& m) : matrix(m) {}
at(int row,int col) const54 virtual ergo_real at(int row, int col) const
55 { return matrix.mat[row + col*matrix.nbast]; }
isSparse() const56 virtual bool isSparse() const { return false; }
asSparse() const57 virtual const Dft::SparseMatrix* asSparse() const
58 { throw std::runtime_error("Cannot converse full to sparse"); }
asFull() const59 virtual const ergo_real* asFull() const
60 { return matrix.mat; }
61 };
62
63 class SparseMatrixWrapper : public Dft::Matrix {
64 const Dft::SparseMatrix& matrix;
65 public:
SparseMatrixWrapper(const Dft::SparseMatrix & m)66 explicit SparseMatrixWrapper(const Dft::SparseMatrix& m) : matrix(m) {}
at(int row,int col) const67 virtual ergo_real at(int row, int col) const
68 { return matrix.at(row, col); }
isSparse() const69 virtual bool isSparse() const { return true; }
asSparse() const70 virtual const Dft::SparseMatrix* asSparse() const
71 { return & matrix; }
asFull() const72 virtual const ergo_real* asFull() const
73 { throw std::runtime_error("cannot convert sparse to full"); }
74 };
75
76
77 Dft::Matrix*
createGridMatrix(const Dft::FullMatrix & mat)78 createGridMatrix(const Dft::FullMatrix& mat)
79 {
80 return new FullMatrixWrapper(mat);
81 }
82
83 Dft::Matrix*
createGridMatrix(const Dft::SparseMatrix & mat)84 createGridMatrix(const Dft::SparseMatrix& mat)
85 {
86 return new SparseMatrixWrapper(mat);
87 }
88
89
90 /* ===================================================================
91 GENERAL GRID MANIPULATION ROUTINES.
92 =================================================================== */
93 /* grid_get_fname, grid_set_tmpdir: helper routines for creating
94 * unique grid file name and setting the temporary directory used for
95 * storing the grid file. */
96
97 std::string grid_tmpdir;
98 static char*
grid_get_fname(const char * base,int filenum)99 grid_get_fname(const char *base, int filenum)
100 {
101 char *res;
102 if(grid_tmpdir.empty()) {
103 char *tmpdir = getenv("TMPDIR");
104 if(tmpdir)
105 grid_tmpdir = tmpdir;
106 else
107 grid_tmpdir = ".";
108 }
109 res = (char*)malloc(strlen(base) + grid_tmpdir.length() + 15);
110 sprintf(res, "%s/%s.%06u.%05d", grid_tmpdir.c_str(),
111 base, (unsigned)getpid(), filenum);
112 return res;
113 }
114
115 void
grid_set_tmpdir(const char * tmpdir)116 grid_set_tmpdir(const char *tmpdir)
117 {
118 grid_tmpdir = tmpdir ? tmpdir : ".";
119 }
120
121 #define GRID_BASE_NAME "ERGO-grid"
122 #define GRID_PATT_NAME "ERGO-patt"
123
124 /* shared grid file */
125 static FILE *grid_file = NULL;
126 static int grid_file_open_count = 0;
127 static char *grid_file_name = NULL;
128 static char *patt_file_name = NULL;
129 static pthread_mutex_t grid_mutex = PTHREAD_MUTEX_INITIALIZER;
130 static pthread_mutex_t grdone_mutex = PTHREAD_MUTEX_INITIALIZER;
131
132
133 /** Frees all the cached data if any. */
134 void
grid_free_files()135 grid_free_files()
136 {
137 if(grid_file_name) {
138 unlink(grid_file_name);
139 free(grid_file_name);
140 grid_file_name = NULL;
141 }
142 if(patt_file_name) {
143 unlink(patt_file_name);
144 free(patt_file_name);
145 patt_file_name = NULL;
146 }
147 }
148
149
150 static int grid_atexit_registered = 0;
151 static void
grid_atexit(void)152 grid_atexit(void)
153 {
154 grid_free_files();
155 }
156
157
158
159 bool
grid_is_ready()160 grid_is_ready() {
161 if (grid_file_name)
162 return true;
163 return false;
164 }
165
166 struct DftGridReader {
167 FILE *f;
168 };
169
170 static const int MY_MPI_NUM = 0; /* no MPI for now. */
171 static void
grid_open_stream(const class GridGenMolInfo & molInfo,const Dft::GridParams & gss,Dft::SparsePattern * pattern,DftGridReader * reader)172 grid_open_stream(const class GridGenMolInfo& molInfo,
173 const Dft::GridParams& gss,
174 Dft::SparsePattern *pattern,
175 DftGridReader *reader)
176 {
177 pthread_mutex_lock(&grdone_mutex);
178 if (!grid_file_name) {
179 grid_file_name = grid_get_fname(GRID_BASE_NAME, MY_MPI_NUM);
180 ErgoGridStream *egStream = grid_stream_new(gss, molInfo);
181 if(pattern)
182 grid_stream_set_sparse_pattern(egStream, pattern);
183 grid_stream_generate(egStream, grid_file_name,
184 dft_get_num_threads());
185 grid_stream_free(egStream);
186 if(pattern) {
187 patt_file_name = grid_get_fname(GRID_PATT_NAME, MY_MPI_NUM);
188 FILE *f;
189 if( (f = fopen(patt_file_name, "wb")) == NULL)
190 throw "Cannot open pattern file for writing";
191 pattern->save(f);
192 if(fclose(f) != 0)
193 throw "Cannot close the sparse pattern file";
194 }
195 }
196
197 /* Grid generated if needed, setup time now! */
198 if(pattern) {
199 FILE *f;
200 if(!patt_file_name)
201 throw "Dft::SparsePattern requested but unavailable.";
202 if( (f = fopen(patt_file_name, "rb")) == NULL)
203 throw "Cannot open pattern file for reading";
204 pattern->load(f);
205 if(fclose(f) != 0)
206 throw "Cannot close the sparse pattern file";
207 }
208 if(!grid_file)
209 grid_file = fopen(grid_file_name, "rb");
210 grid_file_open_count++;
211 pthread_mutex_unlock(&grdone_mutex);
212 reader->f= grid_file;
213 if(grid_file == NULL) {
214 perror("grid_open_standard: DFT quadrature grid file " GRID_BASE_NAME
215 " not found");
216 free(reader);
217 abort();
218 }
219 }
220
221 static void
grid_open_cartesian(const BasisInfoStruct & bis,const Dft::GridParams & gss,const Dft::Matrix * dmat,Dft::SparsePattern * pattern,DftGridReader * reader)222 grid_open_cartesian(const BasisInfoStruct& bis,
223 const Dft::GridParams& gss,
224 const Dft::Matrix* dmat,
225 Dft::SparsePattern *pattern, DftGridReader *reader)
226 {
227 pthread_mutex_lock(&grdone_mutex);
228 if (!grid_file_name) {
229 grid_file_name = grid_get_fname(GRID_BASE_NAME, MY_MPI_NUM);
230 hicu_grid_generate(grid_file_name, bis,
231 gss.hicuParams.maxError,
232 gss.hicuParams.box_size,
233 gss.hicuParams.start_box_size_debug,
234 gss.hicuParams.use_error_per_volume,
235 gss.hicuParams.do_double_checking,
236 gss.hicuParams.compare_to_refined,
237 gss.hicuParams.use_energy_criterion,
238 gss.hicuParams.use_energy_criterion_only,
239 gss.hicuParams.do_variation_checking,
240 dmat, pattern, dft_get_num_threads(), false);
241 if(pattern) {
242 patt_file_name = grid_get_fname(GRID_PATT_NAME, MY_MPI_NUM);
243 FILE *f;
244 if( (f = fopen(patt_file_name, "wb")) == NULL)
245 throw "Cannot open pattern file for writing";
246 pattern->save(f);
247 if(fclose(f) != 0)
248 throw "Cannot close the sparse pattern file";
249 }
250 }
251 /* Grid generated if needed, setup time now! */
252 if(pattern) {
253 FILE *f;
254 if(!patt_file_name)
255 throw "Dft::SparsePattern requested but unavailable.";
256 if( (f = fopen(patt_file_name, "rb")) == NULL)
257 throw "Cannot open pattern file for reading";
258 pattern->load(f);
259 if(fclose(f) != 0)
260 throw "Cannot close the sparse pattern file";
261 }
262
263 if(!grid_file)
264 grid_file = fopen(grid_file_name, "rb");
265 grid_file_open_count++;
266 pthread_mutex_unlock(&grdone_mutex);
267 reader->f= grid_file;
268 if(grid_file == NULL) {
269 perror("grid_open_cartesian: DFT quadrature grid file " GRID_BASE_NAME
270 " not found");
271 free(reader);
272 abort();
273 }
274 }
275
276 /** Returns a handle to a grid file. Sets the sparse pattern if
277 passed. Observe that sparse pattern must be passed the first time
278 to get generated. Otherwise, subsequent calls will not be able to
279 set it. */
280
281 DftGridReader*
grid_open_full(const class GridGenMolInfo * mol_info,const Dft::GridParams & gss,Dft::SparsePattern * pattern,const Dft::Matrix * dmat,const BasisInfoStruct & bis)282 grid_open_full(const class GridGenMolInfo *mol_info,
283 const Dft::GridParams& gss,
284 Dft::SparsePattern *pattern,
285 const Dft::Matrix* dmat,
286 const BasisInfoStruct& bis)
287 {
288 DftGridReader *res = dal_new(1,DftGridReader);
289 if(!grid_atexit_registered) {
290 atexit(grid_atexit);
291 grid_atexit_registered = 1;
292 }
293
294 switch(gss.gridType) {
295 case Dft::GridParams::TYPE_STANDARD:
296 grid_open_stream(*mol_info, gss, pattern, res);
297 return res;
298 case Dft::GridParams::TYPE_HICU:
299 grid_open_cartesian(bis, gss, dmat, pattern, res);
300 return res;
301 default:
302 perror("Error in grid_open: unknown grid type\n");
303 free(res);
304 abort();
305 } /* END SWITCH */
306 }
307
308 /** grid_getchunk_blocked() reads grid data also with screening
309 information if only nblocks and shlblocks are provided.
310
311 @param rawgrid shared grid handle.
312 @param maxlen the upper limit on the grid point chunk length.
313
314 @param nBlocks will contain number of active b.f. blocks. May be
315 NULL if uninteresting.
316
317 @param shlBlocks pointer to the shell block range.
318
319 @param coor array with grid point coordinates.
320 @param weight array with grid point weights.
321
322 @return number of read grid points. -1 on end-of-file.
323 */
324 int
grid_getchunk_blocked(DftGridReader * rawgrid,int maxlen,int * nBlocks,int * shlBlocks,real (* coor)[3],real * weight)325 grid_getchunk_blocked(DftGridReader* rawgrid, int maxlen,
326 int *nBlocks, int *shlBlocks,
327 real (*coor)[3], real *weight)
328 {
329 int points = 0, rc, bl_cnt;
330 FILE *f;
331
332 pthread_mutex_lock(&grid_mutex);
333 f = rawgrid->f;
334 if(fread(&points, sizeof(int), 1, f) <1) {
335 points = -1; /* end of file */
336 goto unlock_and_quit;
337 }
338 if(points>maxlen) {
339 fprintf(stderr,
340 "grid_getchunk: too long vector length in file: %d > %d\n"
341 "Calculation will stop.\n", points, maxlen);
342 throw "grid_getchunk: too long vector length in file.";
343 points = -1; /* stop this! */
344 goto unlock_and_quit;
345 }
346
347 if(fread(&bl_cnt, sizeof(unsigned), 1, f) <1) {
348 puts("OCNT reading error."); points = -1; goto unlock_and_quit;
349 }
350 if(nBlocks) *nBlocks = bl_cnt;
351
352 if(shlBlocks) {
353 rc = fread(shlBlocks, sizeof(int), bl_cnt*2, f);
354 } else {
355 int buf, cnt;
356 rc = 0;
357 for(cnt=0; cnt<bl_cnt*2; cnt+=rc) {
358 rc = fread(&buf, sizeof(int), 1, f);
359 if(rc < 1)
360 break;
361 }
362 }
363 if(rc<1) {
364 fprintf(stderr,
365 "IBLOCKS reading error: failed to read %d blocks.\n",
366 *nBlocks);
367 points = -1; goto unlock_and_quit;
368 }
369
370 if((signed)fread(coor, sizeof(real), 3*points, f) < 3*points) {
371 puts("reading grid point coords failed"); points = -1;
372 goto unlock_and_quit;
373 }
374 if((signed)fread(weight, sizeof(real), points, f) < points) {
375 puts("reading grid weights failed."); points = -1;
376 }
377 unlock_and_quit:
378 pthread_mutex_unlock(&grid_mutex);
379 return points;
380 }
381
382
383 /** Closes the shared grid handle that is specifed as the argument. */
384
385 void
grid_close(DftGridReader * rawgrid)386 grid_close(DftGridReader *rawgrid)
387 {
388 pthread_mutex_lock(&grid_mutex);
389 if(--grid_file_open_count == 0) {
390 fclose(grid_file);
391 grid_file = NULL;
392 }
393 pthread_mutex_unlock(&grid_mutex);
394 free(rawgrid);
395 }
396