20% speedup by taking the sinh and tanh calculations out of the loop.

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

1 files changed, 40 insertions, 34 deletions

diff --git a/src/Misner_standard.c b/src/Misner_standard.c
index 9b9a58a..ce39382 100644
--- a/src/Misner_standard.c
+++ b/src/Misner_standard.c
@@ -62,9 +62,9 @@ void Misner_standard(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-  int i,n;
+  int i, n;
   int npoints;
-  CCTK_REAL csch, coth, inv_r1, inv_r2;
+  CCTK_REAL *csch, *coth, inv_r1, inv_r2;
   CCTK_REAL x_squared, y_squared, xy_squared;
   CCTK_REAL inv_r1_cubed, inv_r2_cubed;
   CCTK_REAL inv_r1_5, inv_r2_5;
@@ -93,6 +93,22 @@ void Misner_standard(CCTK_ARGUMENTS)
     memset (psizz, 0, npoints * sizeof (CCTK_REAL));
   }
 
+  csch = (CCTK_REAL *) malloc (2 * (nmax + 1) * sizeof (CCTK_REAL));
+  coth = csch + nmax + 1;
+
+  /*     compute the ADM mass
+   *     -------------------- 
+   */
+  mass = zero;
+  for(n = 1; n <= nmax; n++)
+  {
+    csch[n] = one / sinh(mu*n);
+    coth[n] = one / tanh(mu*n);
+    mass   += 4.0 * csch[n];
+  }
+  CCTK_VInfo(CCTK_THORNSTRING, "ADM mass is %f", (double) mass);
+
+
   for(i = 0; i < npoints; i++)
   {
     psi [i] = one;
@@ -103,12 +119,10 @@ void Misner_standard(CCTK_ARGUMENTS)
 
     for(n = nmax; n >= 1; n--)
     {
-      csch = one/sinh(mu*n);
-      coth = one/tanh(mu*n);
-      inv_r1 = one / sqrt(xy_squared+SQR(z[i]+coth));
-      inv_r2 = one / sqrt(xy_squared+SQR(z[i]-coth));
+      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*(inv_r1 + inv_r2);
+      psi[i] += csch[n]*(inv_r1 + inv_r2);
           
       if (use_conformal_derivs)
       {
@@ -116,20 +130,20 @@ void Misner_standard(CCTK_ARGUMENTS)
         inv_r2_cubed = inv_r2 * inv_r2 * inv_r2;
         inv_r1_5     = pow(inv_r1, 5);
         inv_r2_5     = pow(inv_r2, 5);
-        psix[i]  +=  -x[i] * (inv_r2_cubed + inv_r1_cubed) * csch;
-        psiy[i]  +=  -y[i] * (inv_r2_cubed + inv_r1_cubed) * csch;
-        psiz[i]  +=  (-(z[i]-coth)*inv_r2_cubed - (z[i]+coth)*inv_r1_cubed) * csch;
+        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;
-        psixy[i] +=  three*x[i]*y[i]*(inv_r1_5 + inv_r2_5) * csch;
-        psixz[i] +=  (three*x[i]*(z[i] - coth)*inv_r2_5 
-                      + three*x[i]*(z[i] + coth)*inv_r1_5) * csch;
+                      - 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 
+                      + 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;
-        psiyz[i] +=  (three*y[i]*(z[i] - coth)*inv_r2_5
-                      + three*y[i]*(z[i] + coth)*inv_r1_5) * csch;
-        psizz[i] += (-inv_r2_cubed+three*SQR(z[i] - coth)*inv_r2_5
-                      + three*SQR(z[i] + coth)*inv_r1_5 - inv_r1_cubed) * csch;
+                      - inv_r1_cubed - inv_r2_cubed) * csch[n];
+        psiyz[i] +=  (three*y[i]*(z[i] - coth[n])*inv_r2_5
+                      + three*y[i]*(z[i] + coth[n])*inv_r1_5) * csch[n];
+        psizz[i] += (-inv_r2_cubed+three*SQR(z[i] - coth[n])*inv_r2_5
+                      + three*SQR(z[i] + coth[n])*inv_r1_5 - inv_r1_cubed) * csch[n];
       }
     }
 
@@ -152,17 +166,6 @@ void Misner_standard(CCTK_ARGUMENTS)
     }
   }
   
-  /*     compute the ADM mass
-   *     -------------------- 
-   */
-  mass = zero;
-  
-  for(n = 1; n <= nmax; n++)
-  {
-    mass += 4.0 / sinh(n*mu);
-  }
-  CCTK_VInfo(CCTK_THORNSTRING, "ADM mass is %f", mass);
-
   /*     Should initialize lapse to Cadez value if possible
    *     --------------------------------------------------
    */
@@ -181,12 +184,11 @@ void Misner_standard(CCTK_ARGUMENTS)
 
       for(n = 1; n <= nmax; n++)
       {
-        coth = one/tanh(mu*n);
-        inv_r1 = one / sqrt(xy_squared+SQR(z[i]+coth));
-        inv_r2 = one / sqrt(xy_squared+SQR(z[i]-coth));
+        inv_r1 = one / sqrt(xy_squared+SQR(z[i]+coth[n]));
+        inv_r2 = one / sqrt(xy_squared+SQR(z[i]-coth[n]));
         powfac = -powfac;
 
-        alp[i] += powfac * one/sinh(mu*n)*(inv_r1 + inv_r2);
+        alp[i] += powfac * csch[n] * (inv_r1 + inv_r2);
       }
 
       alp[i] /= psi[i];
@@ -229,4 +231,8 @@ void Misner_standard(CCTK_ARGUMENTS)
   memset (kxz, 0, npoints * sizeof (CCTK_REAL));
   memset (kyz, 0, npoints * sizeof (CCTK_REAL));
 
+  if (csch)
+  {
+    free (csch);
+  }
 }