Erik's patch to improve parallelisation via OpenMP.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TwoPunctures/trunk@88 b2a53a04-0f4f-0410-87ed-f9f25ced00cf

9 files changed, 171 insertions, 127 deletions

diff --git a/README b/README
index f417643..6016a54 100644
--- a/README
+++ b/README
@@ -1,11 +1,9 @@
-CVS info   : $Header$
-
 Cactus Code Thorn TwoPunctures
 Thorn Author(s)     : Marcus Ansorg <marcus.ansorg@aei.mpg.de>
-                    : Erik Schnetter <schnetter@aei.mpg.de>
-                    : Frank Loeffler <frank.loeffler@aei.mpg.de>
+                    : Erik Schnetter <schnetter@cct.lsu.edu>
+                    : Frank Löffler <knarf@cct.lsu.edu>
 Thorn Maintainer(s) : Marcus Ansorg <marcus.ansorg@aei.mpg.de>
-                    : Erik Schnetter <schnetter@aei.mpg.de>
+                    : Erik Schnetter <schnetter@cct.lsu.edu>
 --------------------------------------------------------------------------
 
 Purpose of the thorn:
@@ -13,14 +11,14 @@ Purpose of the thorn:
 Create initial for two puncture black holes using a single domain
 spectral method.  This method is described in
 
-Marcus Ansorg, Bernd Brügmann, Wolfgang Tichy: A single-domain
-spectral method for black hole puncture data, gr-qc/0404056
-
-This thorn is currently extended to include source terms to be able
-to handle matter with the same restrictions as the BHs before plus
-zero momentum/ang. momentum (K_{ij} = 0)
-This is still work in progess and might work or might not. Is should
-_not_ influence the ability of the thorn without matter. If it does,
-it is a bug and should be reported to the BTS or Frank Loeffler
-(frank.loeffler@aei.mpg.de).
+Marcus Ansorg, Bernd Brügmann, Wolfgang Tichy,
+"A single-domain spectral method for black hole puncture data",
+PRD 70, 064011 (2004),
+arXiv:gr-qc/0404056.
 
+NOTE: This thorn is currently being extended by Frank Löffler to
+include source terms to be able to handle matter, with the same
+restrictions as the black holes before, plus zero momentum/angular
+momentum (K_{ij} = 0).  This is still work in progess and might work
+or might not.  It should _not_ influence the ability of the thorn
+without matter.
diff --git a/doc/documentation.tex b/doc/documentation.tex
index 60ff50f..95f2976 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -2,7 +2,6 @@
 %  Cactus Thorn template for ThornGuide documentation
 %  Author: Ian Kelley
 %  Date: Sun Jun 02, 2002
-%  $Header$
 %
 %  Thorn documentation in the latex file doc/documentation.tex
 %  will be included in ThornGuides built with the Cactus make system.
@@ -65,9 +64,6 @@
 %
 % *======================================================================*
 
-% If you are using CVS use this line to give version information
-% $Header$
-
 \documentclass{article}
 
 % Use the Cactus ThornGuide style file
@@ -79,14 +75,12 @@
 \begin{document}
 
 % The author of the documentation
-\author{Marcus Ansorg \textless marcus.ansorg@aei.mpg.de\textgreater \\ Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+\author{Marcus Ansorg \textless marcus.ansorg@aei.mpg.de\textgreater \\ Erik Schnetter \textless schnetter@cct.lsu.edu\textgreater}
 
 % The title of the document (not necessarily the name of the Thorn)
 \title{}
 
-% the date your document was last changed, if your document is in CVS,
-% please use:
-%    \date{$ $Date$ $}
+% the date your document was last changed:
 \date{May 18 2004}
 
 \maketitle
diff --git a/interface.ccl b/interface.ccl
index 4afe81d..d022ac1 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,5 +1,4 @@
 # Interface definition for thorn TwoPunctures
-# $Header$
 
 IMPLEMENTS: TwoPunctures
 
diff --git a/param.ccl b/param.ccl
index 47e0df4..1b1f504 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,5 +1,4 @@
 # Parameter definitions for thorn TwoPunctures
-# $Header$
 
 SHARES: ADMBase
 
@@ -97,22 +96,22 @@ REAL par_b "x coordinate of the m+ puncture"
 
 REAL par_m_plus "mass of the m+ puncture" STEERABLE = ALWAYS
 {
-  0.0:* :: ""
+  0.0:*) :: ""
 } 1.0
 
 REAL par_m_minus "mass of the m- puncture" STEERABLE = ALWAYS
 {
-  0.0:* :: ""
+  0.0:*) :: ""
 } 1.0
 
 REAL target_M_plus "target ADM mass for m+"
 {
-  (0.0:*) :: ""
+  0.0:*) :: ""
 } 0.5
 
 REAL target_M_minus "target ADM mass for m-"
 {
-  (0.0:*) :: ""
+  0.0:*) :: ""
 } 0.5
 
 REAL par_P_plus[3] "momentum of the m+ puncture"
@@ -173,11 +172,11 @@ BOOLEAN do_initial_debug_output "Output debug information about initial guess"
 {
 } "no"
 
-BOOLEAN multiply_old_lapse "Multiply the old lapse to the new one"
+BOOLEAN multiply_old_lapse "Multiply the old lapse with the new one"
 {
 } "no"
 
-BOOLEAN schedule_in_ADMBase_InitialData "Schedule in or after ADMBase_InitialData"
+BOOLEAN schedule_in_ADMBase_InitialData "Schedule in (instead of after) ADMBase_InitialData"
 {
 } "yes"
 
diff --git a/schedule.ccl b/schedule.ccl
index 468f589..dab7041 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,5 +1,4 @@
 # Schedule definitions for thorn TwoPunctures
-# $Header$
 
 if (CCTK_Equals(initial_data, "twopunctures"))
 {
diff --git a/src/Equations.c b/src/Equations.c
index 1fd8876..e6355c9 100644
--- a/src/Equations.c
+++ b/src/Equations.c
@@ -90,10 +90,10 @@ BY_Aijofxyz (CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL Aij[3][3])
   r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
   r2_plus = sqrt (pow (r2_plus, 2) + pow (TP_epsilon, 4));
   r2_minus = sqrt (pow (r2_minus, 2) + pow (TP_epsilon, 4));
-  if (r2_plus < TP_Tiny)
-    r2_plus = TP_Tiny;
-  if (r2_minus < TP_Tiny)
-    r2_minus = TP_Tiny;
+  if (r2_plus < pow(TP_Tiny,2))
+    r2_plus = pow(TP_Tiny,2);
+  if (r2_minus < pow(TP_Tiny,2))
+    r2_minus = pow(TP_Tiny,2);
   r_plus = sqrt (r2_plus);
   r_minus = sqrt (r2_minus);
   r3_plus = r_plus * r2_plus;
diff --git a/src/Newton.c b/src/Newton.c
index 643404e..77be43f 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -10,33 +10,79 @@
 #include "TwoPunctures.h"
 
 static int
-bicgstab (CCTK_POINTER_TO_CONST cctkGH,
-          int nvar, int n1, int n2, int n3, derivs v,
-          derivs dv, int output, int itmax, CCTK_REAL tol, CCTK_REAL *normres);
+bicgstab (CCTK_POINTER_TO_CONST const cctkGH,
+          int const nvar, int const n1, int const n2, int const n3,
+          derivs v, derivs dv,
+          int const output, int const itmax, CCTK_REAL const tol,
+          CCTK_REAL * restrict const normres);
+static CCTK_REAL
+norm_inf (CCTK_REAL const * restrict const F,
+          int const ntotal);
 static void
-relax (CCTK_REAL *dv, int nvar, int n1, int n2, int n3,
-       CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
+relax (CCTK_REAL * restrict const dv,
+       int const nvar, int const n1, int const n2, int const n3,
+       CCTK_REAL const * restrict const rhs,
+       int const * restrict const ncols,
+       int const * restrict const * restrict const cols,
+       CCTK_REAL const * restrict const * restrict const JFD);
 static void
-resid (CCTK_REAL *res, int ntotal, CCTK_REAL *dv, CCTK_REAL *rhs,
-       int *ncols, int **cols, CCTK_REAL **JFD);
+resid (CCTK_REAL * restrict const res,
+       int const ntotal,
+       CCTK_REAL const * restrict const dv,
+       CCTK_REAL const * restrict const rhs,
+       int const * restrict const ncols,
+       int const * restrict const * restrict const cols,
+       CCTK_REAL const * restrict const * restrict const JFD);
 static void
-LineRelax_al (CCTK_REAL *dv, int j, int k, int nvar, int n1, int n2, int n3,
-	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
+LineRelax_al (CCTK_REAL * restrict const dv,
+              int const j, int const k, int const nvar,
+              int const n1, int const n2, int const n3,
+	      CCTK_REAL const * restrict const rhs,
+              int const * restrict const ncols,
+              int const * restrict const * restrict const cols,
+              CCTK_REAL const * restrict const * restrict const JFD);
 static void
-LineRelax_be (CCTK_REAL *dv, int i, int k, int nvar, int n1, int n2, int n3,
-	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD);
+LineRelax_be (CCTK_REAL * restrict const dv,
+              int const i, int const k, int const nvar,
+              int const n1, int const n2, int const n3,
+	      CCTK_REAL const * restrict const rhs,
+              int const * restrict const ncols,
+              int const * restrict const * restrict const cols,
+              CCTK_REAL const * restrict const * restrict const JFD);
+/* --------------------------------------------------------------------------*/
+static CCTK_REAL
+norm_inf (CCTK_REAL const * restrict const F,
+          int const ntotal)
+{
+  CCTK_REAL dmax = -1;
+#pragma omp parallel
+  {
+    CCTK_REAL dmax1 = -1;
+#pragma omp for
+    for (int j = 0; j < ntotal; j++)
+      if (fabs (F[j]) > dmax1)
+        dmax1 = fabs (F[j]);
+#pragma omp critical
+    if (dmax1 > dmax)
+      dmax = dmax1;
+  }
+  return dmax;
+}
 /* --------------------------------------------------------------------------*/
 static void
-resid (CCTK_REAL *res, int ntotal, CCTK_REAL *dv, CCTK_REAL *rhs,
-       int *ncols, int **cols, CCTK_REAL **JFD)
+resid (CCTK_REAL * restrict const res,
+       int const ntotal,
+       CCTK_REAL const * restrict const dv,
+       CCTK_REAL const * restrict const rhs,
+       int const * restrict const ncols,
+       int const * restrict const * restrict const cols,
+       CCTK_REAL const * restrict const * restrict const JFD)
 {
-  int i, m;
-  CCTK_REAL JFDdv_i;
-
-  for (i = 0; i < ntotal; i++)
+#pragma omp parallel for
+  for (int i = 0; i < ntotal; i++)
   {
-    JFDdv_i = 0;
-    for (m = 0; m < ncols[i]; m++)
+    CCTK_REAL JFDdv_i = 0;
+    for (int m = 0; m < ncols[i]; m++)
       JFDdv_i += JFD[i][m] * dv[cols[i][m]];
     res[i] = rhs[i] - JFDdv_i;
   }
@@ -44,8 +90,13 @@ resid (CCTK_REAL *res, int ntotal, CCTK_REAL *dv, CCTK_REAL *rhs,
 
 /* -------------------------------------------------------------------------*/
 static void
-LineRelax_al (CCTK_REAL *dv, int j, int k, int nvar, int n1, int n2, int n3,
-	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
+LineRelax_al (CCTK_REAL * restrict const dv,
+              int const j, int const k, int const nvar,
+              int const n1, int const n2, int const n3,
+	      CCTK_REAL const * restrict const rhs,
+              int const * restrict const ncols,
+              int const * restrict const * restrict const cols,
+              CCTK_REAL const * restrict const * restrict const JFD)
 {
   int i, m, Ic, Ip, Im, col, ivar;
   CCTK_REAL *a, *b, *c, *r, *u;
@@ -98,8 +149,13 @@ LineRelax_al (CCTK_REAL *dv, int j, int k, int nvar, int n1, int n2, int n3,
 
 /* --------------------------------------------------------------------------*/
 static void
-LineRelax_be (CCTK_REAL *dv, int i, int k, int nvar, int n1, int n2, int n3,
-	      CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
+LineRelax_be (CCTK_REAL * restrict const dv,
+              int const i, int const k, int const nvar,
+              int const n1, int const n2, int const n3,
+	      CCTK_REAL const * restrict const rhs,
+              int const * restrict const ncols,
+              int const * restrict const * restrict const cols,
+              CCTK_REAL const * restrict const * restrict const JFD)
 {
   int j, m, Ic, Ip, Im, col, ivar;
   CCTK_REAL *a, *b, *c, *r, *u;
@@ -152,8 +208,12 @@ LineRelax_be (CCTK_REAL *dv, int i, int k, int nvar, int n1, int n2, int n3,
 
 /* --------------------------------------------------------------------------*/
 static void
-relax (CCTK_REAL *dv, int nvar, int n1, int n2, int n3,
-       CCTK_REAL *rhs, int *ncols, int **cols, CCTK_REAL **JFD)
+relax (CCTK_REAL * restrict const dv,
+       int const nvar, int const n1, int const n2, int const n3,
+       CCTK_REAL const * restrict const rhs,
+       int const * restrict const ncols,
+       int const * restrict const * restrict const cols,
+       CCTK_REAL const * restrict const * restrict const JFD)
 {
   int i, j, k, n;
 
@@ -242,12 +302,14 @@ TestRelax (CCTK_POINTER_TO_CONST cctkGH,
 
 /* --------------------------------------------------------------------------*/
 static int
-bicgstab (CCTK_POINTER_TO_CONST cctkGH,
-          int nvar, int n1, int n2, int n3, derivs v,
-	        derivs dv, int output, int itmax, CCTK_REAL tol, CCTK_REAL *normres)
+bicgstab (CCTK_POINTER_TO_CONST const cctkGH,
+          int const nvar, int const n1, int const n2, int const n3,
+          derivs v, derivs dv,
+          int const output, int const itmax, CCTK_REAL const tol,
+          CCTK_REAL * restrict const normres)
 {
   DECLARE_CCTK_PARAMETERS;
-  int ntotal = n1 * n2 * n3 * nvar, ii, j;
+  int ntotal = n1 * n2 * n3 * nvar, ii;
   CCTK_REAL alpha = 0, beta = 0;
   CCTK_REAL rho = 0, rho1 = 1, rhotol = 1e-50;
   CCTK_REAL omega = 0, omegatol = 1e-50;
@@ -287,7 +349,8 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 
   /* compute initial residual rt = r = F - J*dv */
   J_times_dv (nvar, n1, n2, n3, dv, r, u);
-  for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+  for (int j = 0; j < ntotal; j++)
     rt[j] = r[j] = F[j] - r[j];
 
   *normres = norm2 (r, ntotal);
@@ -308,31 +371,38 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 
     /* compute direction vector p */
     if (ii == 0)
-      for (j = 0; j < ntotal; j++)
+    {
+#pragma omp parallel for
+      for (int j = 0; j < ntotal; j++)
 	p[j] = r[j];
+    }
     else
     {
       beta = (rho / rho1) * (alpha / omega);
-      for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+      for (int j = 0; j < ntotal; j++)
 	p[j] = r[j] + beta * (p[j] - omega * vv[j]);
     }
 
     /* compute direction adjusting vector ph and scalar alpha */
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
       ph.d0[j] = 0;
-    for (j = 0; j < NRELAX; j++)	/* solves JFD*ph = p by relaxation*/
+    for (int j = 0; j < NRELAX; j++)	/* solves JFD*ph = p by relaxation*/
       relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
 
     J_times_dv (nvar, n1, n2, n3, ph, vv, u);	/* vv=J*ph*/
     alpha = rho / scalarproduct (rt, vv, ntotal);
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
       s[j] = r[j] - alpha * vv[j];
 
     /* early check of tolerance */
     *normres = norm2 (s, ntotal);
     if (*normres <= tol)
     {
-      for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+      for (int j = 0; j < ntotal; j++)
 	dv.d0[j] += alpha * ph.d0[j];
       if (output == 1) {
 	printf ("bicgstab: %5d  %10.3e  %10.3e  %10.3e  %10.3e\n",
@@ -343,16 +413,18 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
     }
 
     /* compute stabilizer vector sh and scalar omega */
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
       sh.d0[j] = 0;
-    for (j = 0; j < NRELAX; j++)	/* solves JFD*sh = s by relaxation*/
+    for (int j = 0; j < NRELAX; j++)	/* solves JFD*sh = s by relaxation*/
       relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
 
     J_times_dv (nvar, n1, n2, n3, sh, t, u);	/* t=J*sh*/
     omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
 
     /* compute new solution approximation */
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
     {
       dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
       r[j] = s[j] - omega * t[j];
@@ -407,13 +479,14 @@ bicgstab (CCTK_POINTER_TO_CONST cctkGH,
 
 /* -------------------------------------------------------------------*/
 void
-Newton (CCTK_POINTER_TO_CONST cctkGH,
-  int nvar, int n1, int n2, int n3,
-	derivs v, CCTK_REAL tol, int itmax)
+Newton (CCTK_POINTER_TO_CONST const cctkGH,
+        int const nvar, int const n1, int const n2, int const n3,
+	derivs v,
+        CCTK_REAL const tol, int const itmax)
 {
   DECLARE_CCTK_PARAMETERS;
 
-  int ntotal = n1 * n2 * n3 * nvar, ii, j, it;
+  int ntotal = n1 * n2 * n3 * nvar, ii, it;
   CCTK_REAL *F, dmax, normres;
   derivs u, dv;
 
@@ -429,12 +502,10 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
     if (it == 0)
     {
       F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
-      dmax = -1;
-      for (j = 0; j < ntotal; j++)
-        if (fabs (F[j]) > dmax)
-          dmax = fabs (F[j]);
+      dmax = norm_inf (F, ntotal);
     }
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
       dv.d0[j] = 0;
     printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
     printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
@@ -442,26 +513,17 @@ Newton (CCTK_POINTER_TO_CONST cctkGH,
     ii =
       bicgstab (cctkGH,
                 nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres);
-    for (j = 0; j < ntotal; j++)
+#pragma omp parallel for
+    for (int j = 0; j < ntotal; j++)
       v.d0[j] -= dv.d0[j];
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
-    dmax = -1;
-    for (j = 0; j < ntotal; j++)
-    {
-      if (fabs (F[j]) > dmax)
-      {
-        dmax = fabs (F[j]);
-      }
-    }
+    dmax = norm_inf (F, ntotal);
     it += 1;
   }
   if (itmax==0)
   {
       F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
-      dmax = -1;
-      for (j = 0; j < ntotal; j++)
-        if (fabs (F[j]) > dmax)
-          dmax = fabs (F[j]);
+      dmax = norm_inf (F, ntotal);
   }
   printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
   fflush(stdout);
diff --git a/src/TwoPunctures.c b/src/TwoPunctures.c
index b2d687c..d8002ae 100644
--- a/src/TwoPunctures.c
+++ b/src/TwoPunctures.c
@@ -16,7 +16,7 @@
 
 /* Swap two variables */
 static inline
-void swap (CCTK_REAL * const a, CCTK_REAL * const b)
+void swap (CCTK_REAL * restrict const a, CCTK_REAL * restrict const b)
 {
   CCTK_REAL const t = *a; *a=*b; *b=t;
 }
@@ -25,6 +25,7 @@ void swap (CCTK_REAL * const a, CCTK_REAL * const b)
 
 
 
+static
 void set_initial_guess(CCTK_POINTER_TO_CONST cctkGH,
                        derivs v)
 {
@@ -205,15 +206,16 @@ TwoPunctures (CCTK_ARGUMENTS)
 
   int antisymmetric_lapse, averaged_lapse, pmn_lapse, brownsville_lapse;
 
-  int nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
+  int const nvar = 1, n1 = npoints_A, n2 = npoints_B, n3 = npoints_phi;
 
-  int imin[3], imax[3], d;
-  int i, j, k, l, n, ntotal = n1 * n2 * n3 * nvar;
+  int imin[3], imax[3];
+  int const ntotal = n1 * n2 * n3 * nvar;
+#if 0
   int percent10 = 0;
+#endif
   static CCTK_REAL *F = NULL;
   static derivs u, v;
   CCTK_REAL admMass, M_p, M_m, tmp_p, tmp_m, um, up, new_mass, old_mass;
-  CCTK_REAL old_alp;
 
   if (! F) {
     /* Solve only when called for the first time */
@@ -228,7 +230,7 @@ TwoPunctures (CCTK_ARGUMENTS)
     }
     
     /* initialise to 0 */
-    for (j = 0; j < ntotal; j++)
+    for (int j = 0; j < ntotal; j++)
     {
       v.d0[j] = 0.0;
       v.d1[j] = 0.0;
@@ -267,7 +269,7 @@ TwoPunctures (CCTK_ARGUMENTS)
 	um  = PunctIntPolAtArbitPosition 
 	  (0, nvar, n1, n2, n3, v, -par_b, 0., 0.);
 
-	for(l=0; l<32;l++) {
+	for(int l=0; l<32;l++) {
 	  tmp_p=(4*par_b*um+ *mp +4*par_b);
 	  tmp_m=(4*par_b*up+ *mm +4*par_b);
 	  *mp = *mp - ((*mp*(up+1+M_m/tmp_p)-M_p)/
@@ -291,18 +293,15 @@ TwoPunctures (CCTK_ARGUMENTS)
     Newton (cctkGH, nvar, n1, n2, n3, v, Newton_tol, Newton_maxit);
  
     F_of_v (cctkGH, nvar, n1, n2, n3, v, F, u);
-    CCTK_VInfo (CCTK_THORNSTRING,
-                "The two puncture masses are %g and %g",
-                (double) *mm, (double)* mp);
 
     CCTK_VInfo (CCTK_THORNSTRING,
-		  "The puncture masses are %g and %g",
+		  "The two puncture masses are %g and %g",
                 (double) par_m_minus, (double) par_m_plus);
 
     /* print out ADM mass, eq.: \Delta M_ADM=2*r*u=4*b*V for A=1,B=0,phi=0 */
     admMass = (*mp + *mm
                - 4*par_b*PunctEvalAtArbitPosition(v.d0, 1, 0, 0, n1, n2, n3));
-    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g\n", (double)admMass);
+    CCTK_VInfo (CCTK_THORNSTRING, "ADM mass is %g", (double)admMass);
   }
 
   if (CCTK_EQUALS(grid_setup_method, "Taylor expansion"))
@@ -344,7 +343,7 @@ TwoPunctures (CCTK_ARGUMENTS)
     *conformal_state = 0;
   }
 
-  for (d = 0; d < 3; ++ d)
+  for (int d = 0; d < 3; ++ d)
   {
     /*
     imin[d] = 0           + (cctk_bbox[2*d  ] ? 0 : cctk_nghostzones[d]);
@@ -354,22 +353,15 @@ TwoPunctures (CCTK_ARGUMENTS)
     imax[d] = cctk_lsh[d];
   }
 
-  for (k = imin[2]; k < imax[2]; ++k)
-    for (j = imin[1]; j < imax[1]; ++j)
-      for (i = imin[0]; i < imax[0]; ++i)
-        {
-          int ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
-          x[ijk] -= center_offset[0];
-          y[ijk] -= center_offset[1];
-          z[ijk] -= center_offset[2];
-        }
-
-  for (k = imin[2]; k < imax[2]; ++k)
+#pragma omp parallel for
+  for (int k = imin[2]; k < imax[2]; ++k)
   {
-    for (j = imin[1]; j < imax[1]; ++j)
+    for (int j = imin[1]; j < imax[1]; ++j)
     {
-      for (i = imin[0]; i < imax[0]; ++i)
+      for (int i = imin[0]; i < imax[0]; ++i)
       {
+#if 0
+        /* We can't output this when running in parallel */
         if (percent10 != 10*(i+j*cctk_lsh[0]+k*cctk_lsh[0]*cctk_lsh[1]) /
                             (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]))
         {
@@ -377,13 +369,14 @@ TwoPunctures (CCTK_ARGUMENTS)
                            (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]);
             CCTK_VInfo(CCTK_THORNSTRING, "%3d%% done", percent10*10);
         }
+#endif
 
         const int ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
 
         CCTK_REAL x1, y1, z1;
-        x1 = x[ind];
-        y1 = y[ind];
-        z1 = z[ind];
+        x1 = x[ind] - center_offset[0];
+        y1 = y[ind] - center_offset[1];
+        z1 = z[ind] - center_offset[2];
 
         /* We implement swapping the x and z coordinates as follows.
            The bulk of the code that performs the actual calculations
@@ -443,6 +436,7 @@ TwoPunctures (CCTK_ARGUMENTS)
         CCTK_REAL Aij[3][3];
         BY_Aijofxyz (x1, y1, z1, Aij);
 
+        CCTK_REAL old_alp;
         if (multiply_old_lapse)
             old_alp = alp[ind];
 
diff --git a/src/make.code.defn b/src/make.code.defn
index 512024e..b0ebf81 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,5 +1,4 @@
 # Main make.code.defn file for thorn TwoPunctures
-# $Header$
 
 # Source files in this directory
 SRCS = CoordTransf.c Equations.c FuncAndJacobian.c Newton.c TwoPunctures.c TP_utilities.c ParamCheck.c