144 const bool orthorhombic,
const int border_mask,
const double ra[3],
145 const double rab[3],
const double radius,
const grid_basis_set *basis_set_a,
147 const int la_max,
const int lb_max,
const int cycles,
148 const int cycles_per_block,
const int npts_global[][3],
149 const int npts_local[][3],
const int shift_local[][3],
150 const int border_width[][3],
const double dh[][3][3],
153 const int ntasks = cycles;
154 const int nlevels = 1;
155 const int natoms = 2;
156 const int nkinds = 2;
157 int nblocks = cycles / cycles_per_block + 1;
164 int block_offsets[nblocks];
165 memset(block_offsets, 0, nblocks *
sizeof(
int));
166 const double atom_positions[2][3] = {
167 {ra[0], ra[1], ra[2]}, {rab[0] + ra[0], rab[1] + ra[1], rab[2] + ra[2]}};
168 const int atom_kinds[2] = {1, 2};
170 const int ipgf = o1 / ncoset(la_max) + 1;
171 const int jpgf = o2 / ncoset(lb_max) + 1;
172 assert(o1 == (ipgf - 1) * ncoset(la_max));
173 assert(o2 == (jpgf - 1) * ncoset(lb_max));
175 int level_list[ntasks], iatom_list[ntasks], jatom_list[ntasks];
176 int iset_list[ntasks], jset_list[ntasks], ipgf_list[ntasks],
178 int border_mask_list[ntasks], block_num_list[ntasks];
179 double radius_list[ntasks], rab_list_mutable[ntasks][3];
180 for (
int i = 0;
i < cycles;
i++) {
188 border_mask_list[
i] = border_mask;
189 block_num_list[
i] =
i / cycles_per_block + 1;
190 radius_list[
i] = radius;
191 rab_list_mutable[
i][0] = rab[0];
192 rab_list_mutable[
i][1] = rab[1];
193 rab_list_mutable[
i][2] = rab[2];
195 const double(*rab_list)[3] = (
const double(*)[3])rab_list_mutable;
197 grid_create_task_list(
198 orthorhombic, ntasks, nlevels, natoms, nkinds, nblocks, block_offsets,
199 atom_positions, atom_kinds, basis_sets, level_list, iatom_list,
200 jatom_list, iset_list, jset_list, ipgf_list, jpgf_list, border_mask_list,
201 block_num_list, radius_list, rab_list, npts_global,
npts_local,
202 shift_local, border_width, dh, dh_inv, task_list);
210bool grid_replay(
const char *filename,
const int cycles,
const bool collocate,
211 const bool batch,
const int cycles_per_block,
212 const double tolerance) {
215 fprintf(stderr,
"Error: Cycles have to be greater than zero.\n");
219 if (cycles_per_block < 1 || cycles_per_block > cycles) {
221 "Error: Cycles per block has to be between 1 and cycles.\n");
225 FILE *fp = fopen(filename,
"r");
227 fprintf(stderr,
"Could not open task file: %s\n", filename);
231 char header_line[100];
233 if (strcmp(header_line,
"#Grid task v10\n") != 0) {
234 fprintf(stderr,
"Error: Wrong file header.\n");
238 const bool orthorhombic =
parse_int(
"orthorhombic", fp);
239 const int border_mask =
parse_int(
"border_mask", fp);
242 const int la_max =
parse_int(
"la_max", fp);
243 const int la_min =
parse_int(
"la_min", fp);
244 const int lb_max =
parse_int(
"lb_max", fp);
245 const int lb_min =
parse_int(
"lb_min", fp);
250 double dh_mutable[3][3], dh_inv_mutable[3][3], ra[3], rab[3];
255 const double(*dh)[3] = (
const double(*)[3])dh_mutable;
256 const double(*dh_inv)[3] = (
const double(*)[3])dh_inv_mutable;
258 int npts_global[3],
npts_local[3], shift_local[3], border_width[3];
270 double pab_mutable[n2][n1];
271 char line[100], format[100];
272 for (
int i = 0;
i < n2;
i++) {
273 for (
int j = 0; j < n1; j++) {
275 snprintf(format,
sizeof(format),
"pab %i %i %%le",
i, j);
276 assert(sscanf(line, format, &pab_mutable[
i][j]) == 1);
279 const double(*pab)[n1] = (
const double(*)[n1])pab_mutable;
283 offload_create_buffer(npts_local_total, &grid_ref);
284 memset(grid_ref->
host_buffer, 0, npts_local_total *
sizeof(
double));
286 const int ngrid_nonzero =
parse_int(
"ngrid_nonzero", fp);
287 for (
int n = 0; n < ngrid_nonzero; n++) {
291 assert(sscanf(line,
"grid %i %i %i %le", &
i, &j, &k, &value) == 4);
296 double hab_ref[n2][n1];
297 memset(hab_ref, 0, n2 * n1 *
sizeof(
double));
298 for (
int i = o2;
i < ncoset(lb_max) + o2;
i++) {
299 for (
int j = o1; j < ncoset(la_max) + o1; j++) {
301 snprintf(format,
sizeof(format),
"hab %i %i %%le",
i, j);
302 assert(sscanf(line, format, &hab_ref[
i][j]) == 1);
306 double forces_ref[2][3];
310 double virial_ref[3][3];
313 char footer_line[100];
315 if (strcmp(footer_line,
"#THE_END\n") != 0) {
316 fprintf(stderr,
"Error: Wrong footer line.\n");
320 if (fclose(fp) != 0) {
321 fprintf(stderr,
"Could not close task file: %s\n", filename);
326 offload_create_buffer(npts_local_total, &grid_test);
327 double hab_test[n2][n1];
328 double forces_test[2][3];
329 double virial_test[3][3];
330 double start_time, end_time;
338 orthorhombic, border_mask, ra, rab, radius, basisa, basisb, o1, o2,
339 la_max, lb_max, cycles, cycles_per_block, (
const int(*)[3])npts_global,
340 (
const int(*)[3])
npts_local, (
const int(*)[3])shift_local,
341 (
const int(*)[3])border_width, (
const double(*)[3][3])dh,
342 (
const double(*)[3][3])dh_inv, &task_list);
344 offload_create_buffer(n1 * n2, &pab_blocks);
345 offload_create_buffer(n1 * n2, &hab_blocks);
346 const double f = (collocate) ? rscale : 1.0;
347 for (
int i = 0;
i < n1;
i++) {
348 for (
int j = 0; j < n2; j++) {
352 start_time = omp_get_wtime();
353 const int nlevels = 1;
354 const int natoms = 2;
358 grid_collocate_task_list(task_list, func, nlevels,
359 (
const int(*)[3])
npts_local, pab_blocks, grids);
363 grid_integrate_task_list(task_list, compute_tau, natoms, nlevels,
364 (
const int(*)[3])
npts_local, pab_blocks, grids,
365 hab_blocks, forces_test, virial_test);
366 for (
int i = 0;
i < n2;
i++) {
367 for (
int j = 0; j < n1; j++) {
368 hab_test[
i][j] = hab_blocks->host_buffer[
i * n1 + j];
372 end_time = omp_get_wtime();
375 grid_free_task_list(task_list);
376 offload_free_buffer(pab_blocks);
377 offload_free_buffer(hab_blocks);
379 start_time = omp_get_wtime();
382 memset(grid_test->
host_buffer, 0, npts_local_total *
sizeof(
double));
383 for (
int i = 0;
i < cycles;
i++) {
385 orthorhombic, border_mask, func, la_max, la_min, lb_max, lb_min,
386 zeta, zetb, rscale, dh, dh_inv, ra, rab, npts_global,
npts_local,
387 shift_local, border_width, radius, o1, o2, n1, n2, pab,
392 memset(hab_test, 0, n2 * n1 *
sizeof(
double));
393 memset(forces_test, 0, 2 * 3 *
sizeof(
double));
394 double virials_test[2][3][3] = {0};
395 for (
int i = 0;
i < cycles;
i++) {
397 orthorhombic, compute_tau, border_mask, la_max, la_min, lb_max,
398 lb_min, zeta, zetb, dh, dh_inv, ra, rab, npts_global,
npts_local,
399 shift_local, border_width, radius, o1, o2, n1, n2,
400 grid_ref->
host_buffer, hab_test, pab, forces_test, virials_test,
403 for (
int i = 0;
i < 3;
i++) {
404 for (
int j = 0; j < 3; j++) {
405 virial_test[
i][j] = virials_test[0][
i][j] + virials_test[1][
i][j];
409 end_time = omp_get_wtime();
412 double max_value = 0.0;
413 double max_rel_diff = 0.0;
414 const double derivatives_precision = 1e-4;
418 for (
int i = 0;
i < npts_local_total;
i++) {
419 const double ref_value = cycles * grid_ref->
host_buffer[
i];
421 const double diff = fabs(test_value - ref_value);
422 const double rel_diff = diff / fmax(1.0, fabs(ref_value));
423 max_rel_diff = fmax(max_rel_diff, rel_diff);
424 max_value = fmax(max_value, fabs(test_value));
429 for (
int i = 0;
i < n2;
i++) {
430 for (
int j = 0; j < n1; j++) {
431 const double ref_value = cycles * hab_ref[
i][j];
432 const double test_value = hab_test[
i][j];
433 const double diff = fabs(test_value - ref_value);
434 const double rel_diff = diff / fmax(1.0, fabs(ref_value));
435 max_rel_diff = fmax(max_rel_diff, rel_diff);
436 max_value = fmax(max_value, fabs(test_value));
437 if (rel_diff > tolerance) {
438 printf(
"hab[%i, %i] ref: %le test: %le diff:%le rel_diff: %le\n",
i,
439 j, ref_value, test_value, diff, rel_diff);
444 for (
int i = 0;
i < 2;
i++) {
445 for (
int j = 0; j < 3; j++) {
446 const double ref_value = cycles * forces_ref[
i][j];
447 const double test_value = forces_test[
i][j];
448 const double diff = fabs(test_value - ref_value);
449 const double rel_diff = diff / fmax(1.0, fabs(ref_value));
450 max_rel_diff = fmax(max_rel_diff, rel_diff * derivatives_precision);
451 if (rel_diff * derivatives_precision > tolerance) {
452 printf(
"forces[%i, %i] ref: %le test: %le diff:%le rel_diff: %le\n",
453 i, j, ref_value, test_value, diff, rel_diff);
458 for (
int i = 0;
i < 3;
i++) {
459 for (
int j = 0; j < 3; j++) {
460 const double ref_value = cycles * virial_ref[
i][j];
461 const double test_value = virial_test[
i][j];
462 const double diff = fabs(test_value - ref_value);
463 const double rel_diff = diff / fmax(1.0, fabs(ref_value));
464 max_rel_diff = fmax(max_rel_diff, rel_diff * derivatives_precision);
465 if (rel_diff * derivatives_precision > tolerance) {
466 printf(
"virial[ %i, %i] ref: %le test: %le diff:%le rel_diff: %le\n",
467 i, j, ref_value, test_value, diff, rel_diff);
472 printf(
"Task: %-55s %9s %-7s Cycles: %e Max value: %le "
473 "Max rel diff: %le Time: %le sec\n",
474 filename, collocate ?
"Collocate" :
"Integrate",
475 batch ?
"Batched" :
"PGF-CPU", (float)cycles, max_value, max_rel_diff,
476 end_time - start_time);
478 offload_free_buffer(grid_ref);
479 offload_free_buffer(grid_test);
482 if (fetestexcept(FE_DIVBYZERO) != 0) {
483 fprintf(stderr,
"Error: Floating point exception FE_DIVBYZERO.\n");
486 if (fetestexcept(FE_OVERFLOW) != 0) {
487 fprintf(stderr,
"Error: Floating point exception FE_OVERFLOW.\n");
491 return max_rel_diff < tolerance;
void grid_cpu_collocate_pgf_product(const bool orthorhombic, const int border_mask, const enum grid_func func, const int la_max, const int la_min, const int lb_max, const int lb_min, const double zeta, const double zetb, const double rscale, const double dh[3][3], const double dh_inv[3][3], const double ra[3], const double rab[3], const int npts_global[3], const int npts_local[3], const int shift_local[3], const int border_width[3], const double radius, const int o1, const int o2, const int n1, const int n2, const double pab[n2][n1], double *grid)
Public entry point. A thin wrapper with the only purpose of calling write_task_file when DUMP_TASKS =...
void grid_cpu_integrate_pgf_product(const bool orthorhombic, const bool compute_tau, const int border_mask, const int la_max, const int la_min, const int lb_max, const int lb_min, const double zeta, const double zetb, const double dh[3][3], const double dh_inv[3][3], const double ra[3], const double rab[3], const int npts_global[3], const int npts_local[3], const int shift_local[3], const int border_width[3], const double radius, const int o1, const int o2, const int n1, const int n2, const double *grid, double hab[n2][n1], const double pab[n2][n1], double forces[2][3], double virials[2][3][3], double hdab[n2][n1][3], double hadb[n2][n1][3], double a_hdab[n2][n1][3][3])
Integrates a single task. See grid_cpu_integrate.h for details.
static void create_dummy_task_list(const bool orthorhombic, const int border_mask, const double ra[3], const double rab[3], const double radius, const grid_basis_set *basis_set_a, const grid_basis_set *basis_set_b, const int o1, const int o2, const int la_max, const int lb_max, const int cycles, const int cycles_per_block, const int npts_global[][3], const int npts_local[][3], const int shift_local[][3], const int border_width[][3], const double dh[][3][3], const double dh_inv[][3][3], grid_task_list **task_list)
Creates mock task list with one task per cycle.