[Commits] [svn:einsteintoolkit] TwoPunctures/branches/momentum_constraint/src/ (Rev. 114)

[roland.haas at physics.gatech.edu]

User: rhaas
Date: 2010/07/19 03:57 PM

Modified:
 /branches/momentum_constraint/src/
  TP_utilities.c

Log:
 TwoPunctures: add caching of cos() to chebft_Zeros
 
 For grid_setup_method = "evolution" and npoints_A=npoints_B=20, npoints_phi=4
 TwoPunctures spends about 90% (of about 2h) of the time in chebft_Zeros. With
 this patch it is about 77% (of 10min).
 
 For taylor-expansion most of the time is spend in LinRelax, which is fine.
 
 Best might be to actually construct the full set of spectral coefficients and
 then only to chebev/fourev instead of always going from collocation points to
 spectral rep. then evaluate (trades memory for speed).

File Changes:

Directory: /branches/momentum_constraint/src/
=============================================

File [modified]: TP_utilities.c
Delta lines: +33 -2
===================================================================
--- branches/momentum_constraint/src/TP_utilities.c	2010-07-19 13:46:41 UTC (rev 113)
+++ branches/momentum_constraint/src/TP_utilities.c	2010-07-19 20:57:25 UTC (rev 114)
@@ -244,12 +244,43 @@
 }
 
 /*--------------------------------------------------------------------------*/
+#define COSCACHE_N_MAX 256
+static CCTK_REAL *
+chebft_Zeros_cache (int n)
+  /* caches values of cos (Pi/n * j * (k + 0.5)) */
+{
+  int j, k;
+  CCTK_REAL Pion;
+  static CCTK_REAL *cache[COSCACHE_N_MAX];
+  assert(n <= COSCACHE_N_MAX);
+
+  if(cache[n-1] == NULL)
+  {
+#pragma omp critical
+    if(cache[n-1] == NULL) /* check again in case some other thread allocated storage in omp critical sync */
+    {
+      /* set up new cache */
+      cache[n-1] = dvector(0, n*n-1);
+      /* TODO: add method to free allocated memory */
+      Pion = Pi/n;
+      for (j = 0; j < n; j++)
+        for (k = 0; k < n; k++)
+          cache[n-1][j+k*n] = cos (Pion * j * (k + 0.5));
+    }
+  }
+
+  assert(cache[n-1]);
+  return cache[n-1];
+}
+
+/*--------------------------------------------------------------------------*/
 void
 chebft_Zeros (CCTK_REAL u[], int n, int inv)
     /* eq. 5.8.7 and 5.8.8 at x = (5.8.4) of 2nd edition C++ NR */
 {
   int k, j, isignum;
   CCTK_REAL fac, sum, Pion, *c;
+  CCTK_REAL *Cos = chebft_Zeros_cache(n);
 
   c = dvector (0, n);
   Pion = Pi / n;
@@ -261,7 +292,7 @@
     {
       sum = 0.0;
       for (k = 0; k < n; k++)
-	sum += u[k] * cos (Pion * j * (k + 0.5));
+	sum += u[k] * Cos[j+k*n]; /*cos (Pion * j * (k + 0.5))*/
       c[j] = fac * sum * isignum;
       isignum = -isignum;
     }
@@ -274,7 +305,7 @@
       isignum = 1;
       for (k = 0; k < n; k++)
       {
-	sum += u[k] * cos (Pion * (j + 0.5) * k) * isignum;
+	sum += u[k] * Cos[k+j*n] * isignum; /*cos (Pion * (j + 0.5) * k)*/
 	isignum = -isignum;
       }
       c[j] = sum;