Fixed wrong argument types to pow(3) which caused the testsuites on DEC Alphas

(with native C compiler) to fail if the second argument to this function
was given as an integer.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAnalyticBH/trunk@84 6a3ddf76-46e1-4315-99d9-bc56cac1ef84

4 files changed, 213 insertions, 318 deletions

diff --git a/src/BrillLindquist.c b/src/BrillLindquist.c
index 0edc302..5d47e9c 100644
--- a/src/BrillLindquist.c
+++ b/src/BrillLindquist.c
@@ -1,13 +1,15 @@
  /*@@
    @file      BrillLindquist.F
-   @date      
-   @author    
-   @desc 
-      Set up initial data for Brill Lindquist black holes
-   @enddesc 
+   @date
+   @author
+   @desc
+              Set up initial data for Brill Lindquist black holes
+   @enddesc
+   @version   $Id$
  @@*/
 
 #include <math.h>
+#include <string.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -15,6 +17,10 @@
 
 #include "CactusEinstein/Einstein/src/Einstein.h"
 
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusEinstein_IDAnalyticBH_BrillLindquist_c)
+
+
 #define SQR(a) ((a)*(a))
 
 #define MAX_HOLES 4
@@ -23,19 +29,12 @@ void BrillLindquist(CCTK_ARGUMENTS);
 
  /*@@
    @routine    BrillLindquist
-   @date       
-   @author     
-   @desc 
-      Set up initial data for Brill Lindquist black holes
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
- 
-   @endhistory 
-
+   @date
+   @author
+   @desc
+               Set up initial data for Brill Lindquist black holes
+   @enddesc
 @@*/
-
 void BrillLindquist(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS
@@ -47,13 +46,10 @@ void BrillLindquist(CCTK_ARGUMENTS)
   CCTK_REAL tmp1, tmp2, tmp3;
   CCTK_REAL xval, yval, zval;
   CCTK_REAL x_2, y_2, z_2;
-  int i,j,k;
-  int nx,ny,nz;
-  int index;
+  int i, npoints;
+
 
-  nx = cctk_lsh[0];
-  ny = cctk_lsh[1];
-  nz = cctk_lsh[2];
+  npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
 
   /*     Put parameters into arrays for following calculations
    *     -----------------------------------------------------
@@ -79,84 +75,67 @@ void BrillLindquist(CCTK_ARGUMENTS)
   hole_mass[3] = bl_M_4;
 
 
-  for(k=0; k < nz; k++)
+  if (use_conformal_derivs == 1)
   {
-    for(j=0; j < ny; j++)
+    memset (psix, 0, npoints * sizeof (CCTK_REAL));
+    memset (psiy, 0, npoints * sizeof (CCTK_REAL));
+    memset (psiz, 0, npoints * sizeof (CCTK_REAL));
+    memset (psixx, 0, npoints * sizeof (CCTK_REAL));
+    memset (psixy, 0, npoints * sizeof (CCTK_REAL));
+    memset (psixz, 0, npoints * sizeof (CCTK_REAL));
+    memset (psiyy, 0, npoints * sizeof (CCTK_REAL));
+    memset (psiyz, 0, npoints * sizeof (CCTK_REAL));
+    memset (psizz, 0, npoints * sizeof (CCTK_REAL));
+  }
+
+  for (i = 0; i < npoints; i++)
+  {
+    /*     Initialize to zero and then use +=
+     *     ----------------------------------
+     */
+
+    psi[i] = 1.0;
+
+    x_2 = SQR(x[i]);
+    y_2 = SQR(y[i]);
+    z_2 = SQR(z[i]);
+
+    xval = x[i];
+    yval = y[i];
+    zval = z[i];
+
+    for(n = 0;  n < bl_nbh; n++)
     {
-      for(i=0; i < nx; i++)
+
+      /*     Maple Output
+       *     ------------
+       */
+
+      tmp1 = sqrt(x_2+2.0*xval*hole_x0[n]+SQR(hole_x0[n])
+                 +y_2+2.0*yval*hole_y0[n]+SQR(hole_y0[n])
+                 +z_2+2.0*zval*hole_z0[n]+SQR(hole_z0[n]));
+
+      psi[i] += hole_mass[n]/tmp1*0.5;
+
+      if (use_conformal_derivs == 1)
       {
-        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-
-        /*     Initialize to zero and then use +=
-         *     ----------------------------------
-         */
-
-        psi[index] = 1.0;
-
-
-        if (use_conformal_derivs == 1)
-        {
-          psix[index]  = 0.0;
-          psiy[index]  = 0.0;
-          psiz[index]  = 0.0;
-          psixx[index] = 0.0;
-          psixy[index] = 0.0;
-          psixz[index] = 0.0;
-          psiyy[index] = 0.0;
-          psiyz[index] = 0.0;
-          psizz[index] = 0.0;
-        }
-
-        x_2 = SQR(x[index]);
-        y_2 = SQR(y[index]);
-        z_2 = SQR(z[index]);
-
-        xval = x[index];
-        yval = y[index];
-        zval = z[index];
-
-        for(n = 0;  n < bl_nbh; n++)
-        {
-
-          /*     Maple Output
-           *     ------------
-           */
-
-          tmp1 = sqrt(x_2+2.0*xval*hole_x0[n]+SQR(hole_x0[n])
-                     +y_2+2.0*yval*hole_y0[n]+SQR(hole_y0[n])
-                     +z_2+2.0*zval*hole_z0[n]+SQR(hole_z0[n]));
-
-          psi[index] += hole_mass[n]/tmp1/2.0;
-
-          if (use_conformal_derivs == 1)
-          {
-            tmp2 = pow(tmp1, 3);
-
-            psix[index] +=  -hole_mass[n]/tmp2*(2.0*xval+2.0*hole_x0[n])/4.0;
-            psiy[index] +=  -hole_mass[n]/tmp2*(2.0*yval+2.0*hole_y0[n])/4.0;
-            psiz[index] +=  -hole_mass[n]/tmp2*(2.0*zval+2.0*hole_z0[n])/4.0;
-
-            tmp2 = pow(tmp1, 5.0);
-            tmp3 = pow(tmp1, 3.0);
-
-            psixx[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      SQR(2.0*xval+2.0*hole_x0[n])-hole_mass[n]/tmp3/2.0;
-        
-            psixy[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      (2.0*xval+2.0*hole_x0[n])*(2.0*yval+2.0*hole_y0[n]);
-            psixz[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      (2.0*xval+2.0*hole_x0[n])*(2.0*zval+2.0*hole_z0[n]);
-
-            psiyy[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      SQR(2.0*yval+2.0*hole_y0[n])-hole_mass[n]/tmp3/2.0;
-            psiyz[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      (2.0*yval+2.0*hole_y0[n])*(2.0*zval+2.0*hole_z0[n]);
-
-            psizz[index] += 3.0/8.0*hole_mass[n]/tmp2*
-	      SQR(2.0*zval+2.0*hole_z0[n])-hole_mass[n]/tmp3/2.0;
-
-          }
-        }
+        tmp2 = 1 / (tmp1 * tmp1 * tmp1);
+        tmp3 = 4 * (3.0 / 8.0) * hole_mass[n] * tmp2 / (tmp1 * tmp1);
+        tmp2 *= -0.5 * hole_mass[n];
+
+        psix[i] +=  tmp2 * (xval+hole_x0[n]);
+        psiy[i] +=  tmp2 * (yval+hole_y0[n]);
+        psiz[i] +=  tmp2 * (zval+hole_z0[n]);
+
+        psixx[i] += tmp3 * SQR(xval+hole_x0[n]) + tmp2;
+
+        psixy[i] += tmp3 * (xval+hole_x0[n]) * (yval+hole_y0[n]);
+        psixz[i] += tmp3 * (xval+hole_x0[n]) * (zval+hole_z0[n]);
+
+        psiyy[i] += tmp3 * SQR(yval+hole_y0[n]) + tmp2;
+        psiyz[i] += tmp3 * (yval+hole_y0[n]) * (zval+hole_z0[n]);
+
+        psizz[i] += tmp3 * SQR(zval+hole_z0[n]) + tmp2;
       }
     }
   }
@@ -167,25 +146,19 @@ void BrillLindquist(CCTK_ARGUMENTS)
 
   if (use_conformal_derivs == 1)
   {
-    for(k=0; k < nz; k++)
+    for (i = 0; i < npoints; i++)
     {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-          
-          psix[index]  /= psi[index];
-          psiy[index]  /= psi[index];
-          psiz[index]  /= psi[index];
-          psixx[index] /= psi[index];
-          psixy[index] /= psi[index];
-          psixz[index] /= psi[index];
-          psiyy[index] /= psi[index];
-          psiyz[index] /= psi[index];
-          psizz[index] /= psi[index];
-        }
-      }
+      tmp1 = 1 / psi[i];
+
+      psix[i]  *= tmp1;
+      psiy[i]  *= tmp1;
+      psiz[i]  *= tmp1;
+      psixx[i] *= tmp1;
+      psixy[i] *= tmp1;
+      psixz[i] *= tmp1;
+      psiyy[i] *= tmp1;
+      psiyz[i] *= tmp1;
+      psizz[i] *= tmp1;
     }
   }
 
@@ -195,67 +168,37 @@ void BrillLindquist(CCTK_ARGUMENTS)
 
   if (*conformal_state == CONFORMAL_METRIC)
   {
-    for(k=0; k < nz; k++)
+    for (i = 0; i < npoints; i++)
     {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-          
-          gxx[index] = 1.0;
-          gyy[index] = 1.0;
-          gzz[index] = 1.0;
-          gxy[index] = 0.0;
-          gxz[index] = 0.0;
-          gyz[index] = 0.0;
-        }
-      }
+      gxx[i] = 1.0;
+      gyy[i] = 1.0;
+      gzz[i] = 1.0;
+      gxy[i] = 0.0;
+      gxz[i] = 0.0;
+      gyz[i] = 0.0;
     }
   }
-  else 
+  else
   {
-    for(k=0; k < nz; k++)
+    for (i = 0; i < npoints; i++)
     {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-          
-          gxx[index] = pow(psi[index],4);
-          gyy[index] = gxx[index];
-          gzz[index] = gxx[index];
-          gxy[index] = 0.0;
-          gxz[index] = 0.0;
-          gyz[index] = 0.0;
-        }
-      }
+      tmp1 = psi[i] * psi[i];
+      gxx[i] = tmp1 * tmp1;
+      gyy[i] = gxx[i];
+      gzz[i] = gxx[i];
+      gxy[i] = 0.0;
+      gxz[i] = 0.0;
+      gyz[i] = 0.0;
     }
   }
 
   /*     Time-symmetric data
    *     -------------------
    */
-
-  for(k=0; k < nz; k++)
-  {
-    for(j=0; j < ny; j++)
-    {
-      for(i=0; i < nx; i++)
-      {
-        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-          
-        kxx[index] = 0.0;
-        kyy[index] = 0.0;
-        kzz[index] = 0.0;
-        kxy[index] = 0.0;
-        kxz[index] = 0.0;
-        kyz[index] = 0.0;
-      }
-    }
-  }
-
-  return;
-
+  memset (kxx, 0, npoints * sizeof (CCTK_REAL));
+  memset (kyy, 0, npoints * sizeof (CCTK_REAL));
+  memset (kzz, 0, npoints * sizeof (CCTK_REAL));
+  memset (kxy, 0, npoints * sizeof (CCTK_REAL));
+  memset (kxz, 0, npoints * sizeof (CCTK_REAL));
+  memset (kyz, 0, npoints * sizeof (CCTK_REAL));
 }
diff --git a/src/Misner_multiple.c b/src/Misner_multiple.c
index aa50249..d53a3b9 100644
--- a/src/Misner_multiple.c
+++ b/src/Misner_multiple.c
@@ -3,11 +3,12 @@
    @date      
    @author    Carsten Gundlach
    @desc 
-      Set up initial data for multiple Misner black holes
+              Set up initial data for multiple Misner black holes
    @enddesc 
+   @version   $Id$
  @@*/
 
-#include <math.h>
+#include <string.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -15,6 +16,10 @@
 
 #include "CactusEinstein/Einstein/src/Einstein.h"
 
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusEinstein_IDAnalyticBH_Misner_multiple_c)
+
+
 #define SQR(a) ((a)*(a))
 
 void Misner_init(int n, CCTK_REAL mu, int terms);
@@ -26,13 +31,9 @@ void Misner_multiple(CCTK_ARGUMENTS);
    @date       
    @author     Carsten Gundlach
    @desc 
-      Set up initial data for multiple Misner black holes
+               Set up initial data for multiple Misner black holes
    @enddesc 
-   @calls     
-   @history 
- 
-   @endhistory 
-
+   @calls      MisnerEvalPsi
 @@*/
 void Misner_multiple(CCTK_ARGUMENTS)
 {
@@ -148,7 +149,7 @@ void Misner_multiple(CCTK_ARGUMENTS)
   {
     for(i = 0; i < npoints; i++)
     {
-      gxx[i] = pow(psi[i], 4);
+      gxx[i] = psi[i] * psi[i] * psi[i] * psi[i];
       gyy[i] = gxx[i];
       gzz[i] = gxx[i];
     }
diff --git a/src/Misner_standard.c b/src/Misner_standard.c
index ca2a11c..d037599 100644
--- a/src/Misner_standard.c
+++ b/src/Misner_standard.c
@@ -2,18 +2,19 @@
    @file      Misner_standard.c
    @date      March 1997
    @author    Joan Masso
-   @desc 
-      Set up initial data for two Misner black holes
-   @enddesc 
+   @desc
+              Set up initial data for two Misner black holes
+   @enddesc
    @history
    @hdate Sun Oct 17 11:05:48 1999 @hauthor Tom Goodale
    @hdesc Converted to C
    @endhistory
+   @version   $Id$
  @@*/
 
-#include <stdio.h>
-#include <stdlib.h>
 #include <math.h>
+#include <string.h>
+#include <stdlib.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
@@ -21,38 +22,38 @@
 
 #include "CactusEinstein/Einstein/src/Einstein.h"
 
-void Misner_standard(CCTK_ARGUMENTS);
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusEinstein_IDAnalyticBH_Misner_standard_c)
+
 
 #define SQR(a) ((a)*(a))
 
+void Misner_standard(CCTK_ARGUMENTS);
+
  /*@@
    @routine    Misner_standard
-   @date       
+   @date
    @author     Joan Masso, Ed Seidel
-   @desc 
-            Initialize the metric with a time symmetrical 
+   @desc
+            Initialize the metric with a time symmetrical
             black hole spacetime containing
-            two axially symmetric misner black holes with a 
+            two axially symmetric misner black holes with a
             mass/length parameter mu. The mass is computed.
             The spacetime line element has the form:
                  $$ ds^2 = -dt^2 + \Psi^4 (dx^2+dy^2+dz^2) $$
             and only $\Psi$ differs.
             (Conformal factor from Karen Camarda)
-   @enddesc 
-   @calls     
-   @history 
- 
-   @endhistory 
+   @enddesc
 
    @par      mu
-   @pdesc    Misner parameter. 
+   @pdesc    Misner parameter.
    @ptype    real
-   @pcomment Values less than 1.8 do not really correspond to two 
-            black holes, as there is an initial single event horizon 
+   @pcomment Values less than 1.8 do not really correspond to two
+            black holes, as there is an initial single event horizon
             surrounding the throats. So, with low values of mu we also
             have distorted single black holes.
    @endpar
- 
+
    @par      nmax
    @pdesc    Summation limit for the misner series in the 'twobh' case.
    @ptype    integer
@@ -80,7 +81,7 @@ void Misner_standard(CCTK_ARGUMENTS)
 
 
   /*     Initialize so we can accumulate
-   *     ------------------------------- 
+   *     -------------------------------
    */
   if (use_conformal_derivs)
   {
@@ -99,7 +100,7 @@ void Misner_standard(CCTK_ARGUMENTS)
   coth = csch + nmax + 1;
 
   /*     compute the ADM mass
-   *     -------------------- 
+   *     --------------------
    */
   mass = zero;
   for(n = 1; n <= nmax; n++)
@@ -123,22 +124,22 @@ void Misner_standard(CCTK_ARGUMENTS)
     {
       inv_r1 = one / sqrt(xy_squared+SQR(z[i]+coth[n]));
       inv_r2 = one / sqrt(xy_squared+SQR(z[i]-coth[n]));
-          
+
       psi[i] += csch[n]*(inv_r1 + inv_r2);
-          
+
       if (use_conformal_derivs)
       {
         inv_r1_cubed = inv_r1 * inv_r1 * inv_r1;
         inv_r2_cubed = inv_r2 * inv_r2 * inv_r2;
-        inv_r1_5     = pow(inv_r1, 5);
-        inv_r2_5     = pow(inv_r2, 5);
+        inv_r1_5     = inv_r1 * inv_r1 * inv_r1_cubed;
+        inv_r2_5     = inv_r2 * inv_r2 * inv_r2_cubed;
         psix[i]  +=  -x[i] * (inv_r2_cubed + inv_r1_cubed) * csch[n];
         psiy[i]  +=  -y[i] * (inv_r2_cubed + inv_r1_cubed) * csch[n];
         psiz[i]  +=  (-(z[i]-coth[n])*inv_r2_cubed - (z[i]+coth[n])*inv_r1_cubed) * csch[n];
         psixx[i] +=  (three*x_squared*(inv_r1_5 + inv_r2_5)
                       - inv_r1_cubed - inv_r2_cubed) * csch[n];
         psixy[i] +=  three*x[i]*y[i]*(inv_r1_5 + inv_r2_5) * csch[n];
-        psixz[i] +=  (three*x[i]*(z[i] - coth[n])*inv_r2_5 
+        psixz[i] +=  (three*x[i]*(z[i] - coth[n])*inv_r2_5
                       + three*x[i]*(z[i] + coth[n])*inv_r1_5) * csch[n];
         psiyy[i] +=  (three*y_squared*(inv_r1_5 + inv_r2_5)
                       - inv_r1_cubed - inv_r2_cubed) * csch[n];
@@ -167,15 +168,15 @@ void Misner_standard(CCTK_ARGUMENTS)
       psizz[i] *= inv_psi;
     }
   }
-  
+
   /*     Should initialize lapse to Cadez value if possible
    *     --------------------------------------------------
    */
 
   if (CCTK_Equals(initial_lapse,"cadez"))
-  { 
+  {
     CCTK_INFO("Initialise with cadez lapse");
-    
+
     for(i = 0; i < npoints; i++)
     {
       xy_squared = SQR(x[i]) + SQR(y[i]);
@@ -196,7 +197,7 @@ void Misner_standard(CCTK_ARGUMENTS)
       alp[i] /= psi[i];
     }
   }
-  
+
   /*     Metric depends on conformal state
    *     ---------------------------------
    */
@@ -210,11 +211,11 @@ void Misner_standard(CCTK_ARGUMENTS)
       gzz[i] = one;
     }
   }
-  else 
+  else
   {
     for(i = 0; i < npoints; i++)
     {
-      gxx[i] = pow(psi[i],4);
+      gxx[i] = psi[i] * psi[i] * psi[i] * psi[i];
       gyy[i] = gxx[i];
       gzz[i] = gxx[i];
     }
diff --git a/src/Schwarzschild.c b/src/Schwarzschild.c
index bd29ef3..85c0b88 100644
--- a/src/Schwarzschild.c
+++ b/src/Schwarzschild.c
@@ -2,21 +2,24 @@
    @file      Schwarzschild.c
    @date      Sun Oct 17 10:35:41 1999
    @author    Tom Goodale
-   @desc 
-   C version of Scwhwarzschild lapse routine
-   @enddesc 
+   @desc
+              C version of Scwhwarzschild lapse routine
+   @enddesc
+   @version   $Id$
  @@*/
 
-#include <math.h>
+#include <string.h>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-
 /* Need include file from Einstein */
 #include "CactusEinstein/Einstein/src/Einstein.h"
 
+static char *rcsid = "$Header$";
+CCTK_FILEVERSION(CactusEinstein_IDAnalyticBH_Schwarzschild_c)
+
 void Schwarzschild(CCTK_ARGUMENTS);
 
 
@@ -25,134 +28,81 @@ void Schwarzschild(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-  CCTK_REAL zero,one,two,three;
-  CCTK_REAL tmp;
-  CCTK_REAL r_squared;
+  const CCTK_REAL zero = 0.0, one = 1.0, two = 2.0, three = 3.0;
+  CCTK_REAL tmp, r_squared, r_cubed;
+  int i, npoints;
 
-  int i,j,k;
-  int nx,ny,nz;
-  int index;
 
-  zero = 0.0;
-  one = 1.0;
-  two = 2.0;
-  three = 3.0;
+  npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
 
-  nx = cctk_lsh[0];
-  ny = cctk_lsh[1];
-  nz = cctk_lsh[2];
-  
   /*     conformal metric flag */
   if(use_conformal == 1)
   {
-  
-    *conformal_state = CONFORMAL_METRIC;
 
+    *conformal_state = CONFORMAL_METRIC;
 
-    for(k=0; k < nz; k++)
+    for (i = 0; i < npoints; i++)
     {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-
-          /*        Compute conformal factor */
-          psi[index] = ( one + mass/two/r[index]);
-         
-
-          /*        derivatives of psi / psi */
-
-          tmp = mass/two/pow(r[index],3) / psi[index];
-          psix[index] = -x[index]*tmp;
-          psiy[index] = -y[index]*tmp;
-          psiz[index] = -z[index]*tmp;
-         
-          tmp = mass/two/pow(r[index],5)/psi[index];
-          psixy[index] = three*x[index]*y[index]*tmp;
-          psixz[index] = three*x[index]*z[index]*tmp;
-          psiyz[index] = three*y[index]*z[index]*tmp;
-
-          r_squared = r[index]*r[index];
-          psixx[index]  = (three*x[index]*x[index] - r_squared)*tmp;
-          psiyy[index]  = (three*y[index]*y[index] - r_squared)*tmp;
-          psizz[index]  = (three*z[index]*z[index] - r_squared)*tmp;
-         
-          gxx[index] = one;
-          gyy[index] = one;
-          gzz[index] = one;
-          gxy[index] = zero; 
-          gxz[index] = zero;  
-          gyz[index] = zero;  
-        }
-      }
+      /*        Compute conformal factor */
+      psi[i] = ( one + mass/two/r[i]);
+
+      /*        derivatives of psi / psi */
+      r_squared = r[i]*r[i];
+      r_cubed   = r[i]*r_squared;
+      tmp = mass/two/r_cubed/psi[i];
+      psix[i] = -x[i]*tmp;
+      psiy[i] = -y[i]*tmp;
+      psiz[i] = -z[i]*tmp;
+
+      tmp = mass/two/(r_squared*r_cubed)/psi[i];
+      psixy[i] = three*x[i]*y[i]*tmp;
+      psixz[i] = three*x[i]*z[i]*tmp;
+      psiyz[i] = three*y[i]*z[i]*tmp;
+
+      psixx[i]  = (three*x[i]*x[i] - r_squared)*tmp;
+      psiyy[i]  = (three*y[i]*y[i] - r_squared)*tmp;
+      psizz[i]  = (three*z[i]*z[i] - r_squared)*tmp;
+
+      gxx[i] = one;
+      gyy[i] = one;
+      gzz[i] = one;
+      gxy[i] = zero;
+      gxz[i] = zero;
+      gyz[i] = zero;
     }
-  }     
+  }
   else
   {
     *conformal_state = NOCONFORMAL_METRIC;
 
-    for(k=0; k < nz; k++)
+    for (i = 0; i < npoints; i++)
     {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-
-          gxx[index] = pow(( one + mass/two/r[index]), 4);
-          gyy[index] = gxx[index];
-          gzz[index] = gxx[index];
-          gxy[index] = zero;
-          gxz[index] = zero;
-          gyz[index] = zero;
-        }
-      }
+      r_squared = r[i] * r[i];
+      gxx[i] = one + mass/two/(r_squared * r_squared);
+      gyy[i] = gxx[i];
+      gzz[i] = gxx[i];
+      gxy[i] = zero;
+      gxz[i] = zero;
+      gyz[i] = zero;
     }
   }
 
   /*     If the initial lapse is not one ... */
   if (CCTK_Equals(initial_lapse,"schwarz"))
-  { 
+  {
     CCTK_INFO("Initialise with schwarzschild lapse");
 
-    for(k=0; k < nz; k++)
-    {
-      for(j=0; j < ny; j++)
-      {
-        for(i=0; i < nx; i++)
-        {
-          index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-
-          alp[index] = (2.*r[index] - mass)/(2.*r[index]+mass);
-        }
-      }
-    }
-  }
-      
-  /*     time symmetric initial slice */
-
-  for(k=0; k < nz; k++)
-  {
-    for(j=0; j < ny; j++)
+    for (i = 0; i < npoints; i++)
     {
-      for(i=0; i < nx; i++)
-      {
-        index = CCTK_GFINDEX3D(cctkGH, i,j,k);
-        
-        kxx[index] = zero;
-        kxy[index] = zero;
-        kxz[index] = zero;
-        kyy[index] = zero;
-        kyz[index] = zero;
-        kzz[index] = zero;
-      }
+      alp[i] = (2.*r[i] - mass)/(2.*r[i]+mass);
     }
   }
-    
-  return;
-
-  
-  
 
+  /*     time symmetric initial slice */
+  memset (kxx, 0, npoints * sizeof (CCTK_REAL));
+  memset (kxy, 0, npoints * sizeof (CCTK_REAL));
+  memset (kxz, 0, npoints * sizeof (CCTK_REAL));
+  memset (kyy, 0, npoints * sizeof (CCTK_REAL));
+  memset (kyz, 0, npoints * sizeof (CCTK_REAL));
+  memset (kzz, 0, npoints * sizeof (CCTK_REAL));
 }