grid: Convert calculate_rho_elec to C

src/grid/Makefile

@@ -1,17 +1,17 @@
-CFLAGS   = -g -O3 -std=c99 -march=native -Wall -Wextra
+CFLAGS   = -fopenmp -g -O3 -std=c99 -march=native -Wall -Wextra
 
 all: grid_collocate_miniapp.x grid_collocate_unittest.x
 
 grid_collocate_cpu.o: grid_prepare_pab.o
-grid_collocate_replay.o: grid_collocate_cpu.o
+grid_collocate_replay.o: grid_task_list.o grid_collocate_cpu.o
 
 %.o: %.c
 	$(CC) -c $(CFLAGS) $<
 
-grid_collocate_miniapp.x: grid_collocate_miniapp.o grid_collocate_replay.o grid_collocate_cpu.o grid_prepare_pab.o
+grid_collocate_miniapp.x: grid_collocate_miniapp.o grid_collocate_replay.o grid_task_list.o grid_collocate_cpu.o grid_prepare_pab.o
 	$(CC) $(CFLAGS) -o $@ $^ -lm
 
-grid_collocate_unittest.x: grid_collocate_unittest.o grid_collocate_replay.o grid_collocate_cpu.o grid_prepare_pab.o
+grid_collocate_unittest.x: grid_collocate_unittest.o grid_collocate_replay.o grid_task_list.o grid_collocate_cpu.o grid_prepare_pab.o
 	$(CC) $(CFLAGS) -o $@ $^ -lm
 
 clean:

src/grid/grid_collocate_cpu.c

@@ -5,7 +5,6 @@
 
 #include <stdio.h>
 #include <math.h>
-#include <stdbool.h>
 #include <limits.h>
 #include <assert.h>
 #include <string.h>
@@ -776,6 +775,13 @@ static void grid_collocate_internal(const bool orthorhombic,
                                     const double pab[n2][n1],
                                     double* grid){
 
+    // Check if radius is too small to be mapped onto grid of given resolution.
+    double dh_max = 0.0;
+    for (int i=0; i<3; i++)
+        for (int j=0; j<3; j++)
+           dh_max = max(dh_max, fabs(dh[i][j]));
+    if (2.0 * radius < dh_max) return;
+
     const double zetp = zeta + zetb;
     const double f = zetb / zetp;
     const double rab2 = rab[0] * rab[0] + rab[1] * rab[1] + rab[2] * rab[2];
@@ -909,7 +915,7 @@ void grid_collocate_pgf_product_cpu(const bool orthorhombic,
     const bool DUMP_TASKS = false;
 
     double* grid_before = NULL;
-    const int npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
+    const size_t npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
 
     if (DUMP_TASKS) {
         const size_t sizeof_grid = sizeof(double) * npts_local_total;
@@ -973,7 +979,7 @@ void grid_collocate_pgf_product_cpu(const bool orthorhombic,
                               pab,
                               grid);
 
-        for (int i=0; i < npts_local_total; i++) {
+        for (size_t i=0; i < npts_local_total; i++) {
             grid[i] += grid_before[i];
         }
         free(grid_before);

src/grid/grid_collocate_cpu.h

@@ -5,6 +5,8 @@
 #ifndef GRID_COLLOCATE_CPU_H
 #define GRID_COLLOCATE_CPU_H
 
+#include <stdbool.h>
+
 //******************************************************************************
 // \brief Collocates a single task. A task consists of a pair of atoms each
 //        with a position, Gaussian exponent, and a range of angular momentum.
@@ -36,6 +38,8 @@
 // \param pab           The atom-pair's density matrix block P_{ab}.
 //
 // \param grid          The output grid array to collocate into.
+//
+// \author Ole Schuett
 //******************************************************************************
 void grid_collocate_pgf_product_cpu(const bool orthorhombic,
                                     const bool use_subpatch,

src/grid/grid_collocate_miniapp.c

@@ -5,17 +5,28 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <string.h>
 
 #include "grid_collocate_replay.h"
 
 int main(int argc, char *argv[]){
-    if (argc != 3) {
-        printf("Usage: grid_base_ref_miniapp.x <cycles> <task-file>\n");
+    // Parsing of optional args.
+    int iarg = 1;
+
+    bool batch = false;
+    if (iarg < argc && strcmp(argv[iarg], "--batch") == 0) {
+        iarg++;
+        batch = true;
+    }
+
+    // All optional args have been parsed.
+    if (argc - iarg != 2) {
+        fprintf(stderr, "Usage: grid_base_ref_miniapp.x [--batch] <cycles> <task-file>\n");
         return 1;
     }
 
     int cycles;
-    if (sscanf(argv[1], "%i", &cycles) != 1) {
+    if (sscanf(argv[iarg++], "%i", &cycles) != 1) {
         fprintf(stderr, "Error: Could not parse cycles.\n");
         return 1;
     }
@@ -24,8 +35,8 @@ int main(int argc, char *argv[]){
         return 1;
     }
 
-    const double max_diff = grid_collocate_replay(argv[2], cycles);
-    if (max_diff > 1e-14 * cycles) {
+    const double max_diff = grid_collocate_replay(argv[iarg++], cycles, batch);
+    if (max_diff > 1e-12 * cycles) {
         fprintf(stderr, "Error: Maximal difference is too large.\n");
         return 2;
     }