use openmp in itd.c and move log.h into jmm/include

examples/itd/itd.c

@@ -1,9 +1,12 @@
+#include <omp.h>
 #include <stdio.h>
 #include <string.h>
 
 #include <jmm/bmesh.h>
 #include <jmm/eik3.h>
+#include <jmm/log.h>
 #include <jmm/mesh3.h>
+#include <jmm/util.h>
 
 static void
 get_evaluation_grid(size_t num_az, size_t num_el, dbl **phi, dbl **theta) {
@@ -21,8 +24,9 @@ get_evaluation_grid(size_t num_az, size_t num_el, dbl **phi, dbl **theta) {
 }
 
 int main(int argc, char const *argv[]) {
-  if (argc != 4) {
-    printf("usage: %s <n> <off_path>\n", argv[0]);
+  if (argc < 4) {
+    printf("usage: %s <n> <off_path> [nthreads]\n", argv[0]);
+    exit(EXIT_FAILURE);
   }
 
   dbl eps = 1e-5;
@@ -36,6 +40,13 @@ int main(int argc, char const *argv[]) {
 
   char const *off_path = argv[2];
 
+  size_t nthreads = omp_get_num_threads();
+  if (argc >= 4)
+    nthreads = atoi(argv[3]);
+  omp_set_num_threads(nthreads);
+
+  printf("running ITD simulation w/ %lu threads\n", nthreads);
+
   /* NOTE: units in mm! */
 
   /* For meshes in `./old_meshes`: */
@@ -52,6 +63,8 @@ int main(int argc, char const *argv[]) {
   dbl *phi_grid, *theta_grid;
   get_evaluation_grid(num_az, num_el, &phi_grid, &theta_grid);
 
+  toc();
+
   mesh3_data_s data;
   mesh3_data_init_from_off_file(&data, off_path, maxvol, verbose);
 
@@ -68,6 +81,8 @@ int main(int argc, char const *argv[]) {
   mesh3_dump_verts(mesh, "verts.bin");
   mesh3_dump_cells(mesh, "cells.bin");
 
+  printf("set up tetrahedron mesh [%.2fs]\n", toc());
+
   /* Solve the eikonal equation for the left ear */
   eik3_s *eik_L;
   eik3_alloc(&eik_L);
@@ -76,6 +91,8 @@ int main(int argc, char const *argv[]) {
   eik3_solve(eik_L);
   eik3_dump_jet(eik_L, "jet_L.bin");
 
+  printf("computed eikonal for left ear [%.2fs]\n", toc());
+
   /* Solve the eikonal equation for the right ear */
   eik3_s *eik_R;
   eik3_alloc(&eik_R);
@@ -84,6 +101,8 @@ int main(int argc, char const *argv[]) {
   eik3_solve(eik_R);
   eik3_dump_jet(eik_R, "jet_R.bin");
 
+  printf("computed eikonal for right ear [%.2fs]\n", toc());
+
   jet31t const *jet_L = eik3_get_jet_ptr(eik_L);
   jet31t const *jet_R = eik3_get_jet_ptr(eik_R);
 
@@ -97,23 +116,36 @@ int main(int argc, char const *argv[]) {
   bmesh33_alloc(&tau_R);
   bmesh33_init_from_mesh3_and_jets(tau_R, mesh, jet_R);
 
+  printf("set up tetrahedral splines for interpolating L/R eikonals [%.2fs]\n", toc());
+
   FILE *fp = NULL;
 
+  dbl *itd = malloc(num_el*num_az*sizeof(dbl));
+
+  { size_t k = 0;
+#pragma omp parallel for
+    for (size_t i = 0; i < num_el; ++i) {
+      dbl el = theta_grid[i];
+      for (size_t j = 0; j < num_az; ++j) {
+        dbl az = phi_grid[j];
+        dbl3 point = {cos(az)*sin(el), sin(az)*sin(el), cos(el)};
+        dbl3_dbl_mul_inplace(point, r_grid);
+        dbl T_L = bmesh33_f(tau_L, point)/c;
+        dbl T_R = bmesh33_f(tau_R, point)/c;
+        itd[k++] = T_R - T_L;
+      }
+    } }
+
+  printf("interpolated ITD to grid [%.2fs]\n", toc());
+
   fp = fopen("itd_grid.bin", "wb");
-  for (size_t i = 0; i < num_el; ++i) {
-    dbl el = theta_grid[i];
-    for (size_t j = 0; j < num_az; ++j) {
-      dbl az = phi_grid[j];
-      dbl3 point = {cos(az)*sin(el), sin(az)*sin(el), cos(el)};
-      dbl3_dbl_mul_inplace(point, r_grid);
-      dbl T_L = bmesh33_f(tau_L, point)/c;
-      dbl T_R = bmesh33_f(tau_R, point)/c;
-      dbl itd = T_R - T_L;
-      fwrite(&itd, sizeof(dbl), 1, fp);
-    }
-  }
+  fwrite(itd, sizeof(dbl), num_el*num_az, fp);
   fclose(fp);
 
+  printf("wrote ITD to itd_grid.bin [%.2fs]\n", toc());
+
+  free(itd);
+
   bmesh33_deinit(tau_R);
   bmesh33_dealloc(&tau_R);
 

examples/itd/meson.build

@@ -1,4 +1,4 @@
-executable('itd', 'itd.c', dependencies : jmm_dep)
+executable('itd', 'itd.c', dependencies : [jmm_dep, openmp_dep])
 
 files = ['scratch.py']
 

meson.build

@@ -15,6 +15,7 @@ fs = import('fs')
 m_dep = meson.get_compiler('c').find_library('m', required : false)
 argp_dep = meson.get_compiler('c').find_library('argp', required : false)
 gsl_dep = dependency('gsl')
+openmp_dep = dependency('openmp')
 tetgen_dep = dependency('tetgen')
 
 jmm_lib_src = [

src/eik2m1.c

@@ -5,11 +5,10 @@
 #include <stdlib.h>
 
 #include <jmm/heap.h>
+#include <jmm/log.h>
 #include <jmm/mat.h>
 #include <jmm/utri21.h>
 
-#include "log.h"
-
 struct eik2m1 {
   mesh22_s const *mesh;
   jet22t *jet;

src/eik2mp.c

@@ -5,11 +5,10 @@
 #include <stdlib.h>
 
 #include <jmm/heap.h>
+#include <jmm/log.h>
 #include <jmm/mat.h>
 #include <jmm/utri21.h>
 
-#include "log.h"
-
 struct eik2mp {
   mesh22_s const *mesh;
   jet22t *jet; // eik_vals, grads from eik_grid plus the hessian
@@ -28,7 +27,7 @@ static dbl value(eik2mp_s const *eik, int l) {
 
 static void setpos(eik2mp_s const *eik, int l, int pos) {
   // changes the index in the priority_queue (all these methods can be
-  // found in priority_queue.c 
+  // found in priority_queue.c
   eik->pos[l] = pos;
 }
 
@@ -80,7 +79,7 @@ void eik2mp_deinit(eik2mp_s *eik) {
 }
 
 size_t eik2mp_peek(eik2mp_s const *eik) {
-  // returns the value in the 0th place in the heap 
+  // returns the value in the 0th place in the heap
   return heap_front(eik->heap);
 }
 
@@ -165,7 +164,7 @@ size_t eik2mp_step(eik2mp_s *eik) {
 
   for (size_t i = 0; i < nvv; ++i) {
     size_t l = vv[i]; // neighbor of x0
-    if (eik->state[l] == FAR) { 
+    if (eik->state[l] == FAR) {
       eik->state[l] = TRIAL;
       heap_insert(eik->heap, l);
     }

src/eik3.c

@@ -13,6 +13,7 @@
 #include <jmm/edge.h>
 #include <jmm/eik3_transport.h>
 #include <jmm/heap.h>
+#include <jmm/log.h>
 #include <jmm/mat.h>
 #include <jmm/mesh1.h>
 #include <jmm/mesh2.h>
@@ -25,7 +26,6 @@
 #include <jmm/utri_cache.h>
 #include <jmm/vec.h>
 
-#include "log.h"
 #include "macros.h"
 
 static bool l_OK(size_t l) {

src/eik3_transport.c

@@ -3,11 +3,10 @@
 #include <assert.h>
 #include <math.h>
 
+#include <jmm/log.h>
 #include <jmm/mesh3.h>
 #include <jmm/util.h>
 
-#include "log.h"
-
 static void transport_dbl(eik3_s const *eik, size_t l0, dbl *values) {
   par3_s par = eik3_get_par(eik, l0);
 

src/log.c

@@ -20,7 +20,7 @@
  * IN THE SOFTWARE.
  */
 
-#include "log.h"
+#include <jmm/log.h>
 
 #define MAX_CALLBACKS 32
 

src/mesh2.c

@@ -7,9 +7,9 @@
 #include <string.h>
 
 #include <jmm/def.h>
+#include <jmm/log.h>
 #include <jmm/vec.h>
 
-#include "log.h"
 #include "macros.h"
 #include "mesh_util.h"
 

src/mesh22.c

@@ -4,10 +4,9 @@
 #include <stdlib.h>
 #include <string.h>
 
+#include <jmm/log.h>
 #include <jmm/util.h>
 
-#include "log.h"
-
 struct mesh22 {
   dbl2 *verts;
   size_t nverts;

src/rtree.c

@@ -7,12 +7,12 @@
 
 #include <jmm/bmesh.h>
 #include <jmm/def.h>
+#include <jmm/log.h>
 #include <jmm/mesh2.h>
 #include <jmm/mesh3.h>
 #include <jmm/stats.h>
 #include <jmm/vec.h>
 
-#include "log.h"
 #include "macros.h"
 #include "pool.h"
 

src/uline.c

@@ -5,10 +5,9 @@
 #include <stdlib.h>
 
 #include <jmm/eik3.h>
+#include <jmm/log.h>
 #include <jmm/mesh3.h>
 
-#include "log.h"
-
 static size_t MAX_NUM_ITER = 100;
 
 struct uline {

src/utetra.c

@@ -8,13 +8,13 @@
 #include <jmm/array.h>
 #include <jmm/bb.h>
 #include <jmm/eik3.h>
+#include <jmm/log.h>
 #include <jmm/mat.h>
 #include <jmm/mesh3.h>
 #include <jmm/opt.h>
 #include <jmm/uline.h>
 #include <jmm/util.h>
 
-#include "log.h"
 #include "macros.h"
 
 #define MAX_NITER 100

src/utri21.c

@@ -3,10 +3,9 @@
 #include <assert.h>
 
 #include <jmm/hybrid.h>
+#include <jmm/log.h>
 #include <jmm/mat.h>
 
-#include "log.h"
-
 void utri21_init(utri21_s *utri, dbl2 const xhat, dbl2 const x[2],
                  jet21t const jet[2]) {
   dbl2_copy(xhat, utri->xhat);