remove some code from numerical recipes and replace by code using the gnu

scientific library


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

4 files changed, 81 insertions, 243 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..21d4fd8
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn TwoPunctures
+
+REQUIRES GSL
diff --git a/src/Newton.c b/src/Newton.c
index 9140c08..54fad9e 100644
--- a/src/Newton.c
+++ b/src/Newton.c
@@ -5,6 +5,8 @@
 #include <string.h>
 #include <math.h>
 #include <ctype.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
 #include "cctk_Parameters.h"
 #include "TP_utilities.h"
 #include "TwoPunctures.h"
@@ -99,52 +101,53 @@ LineRelax_al (CCTK_REAL * restrict const dv,
               CCTK_REAL const * restrict const * restrict const JFD)
 {
   int i, m, Ic, Ip, Im, col, ivar;
-  CCTK_REAL *a, *b, *c, *r, *u;
-  a = dvector (0, n1 - 1);
-  b = dvector (0, n1 - 1);
-  c = dvector (0, n1 - 1);
-  r = dvector (0, n1 - 1);
-  u = dvector (0, n1 - 1);
+
+  gsl_vector *diag = gsl_vector_alloc(n1);
+  gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
+  gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
+  gsl_vector *b = gsl_vector_alloc(n1);   /* rhs */
+  gsl_vector *x = gsl_vector_alloc(n1);   /* solution vector */
 
   for (ivar = 0; ivar < nvar; ivar++)
   {
+    gsl_vector_set_zero(diag);
+    gsl_vector_set_zero(e);
+    gsl_vector_set_zero(f);
     for (i = 0; i < n1; i++)
     {
       Ip = Index (ivar, i + 1, j, k, nvar, n1, n2, n3);
       Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
       Im = Index (ivar, i - 1, j, k, nvar, n1, n2, n3);
-      a[i] = 0;
-      b[i] = 0;
-      c[i] = 0;
-      r[i] = rhs[Ic];
+      gsl_vector_set(b,i,rhs[Ic]);
       for (m = 0; m < ncols[Ic]; m++)
       {
 	col = cols[Ic][m];
 	if (col != Ip && col != Ic && col != Im)
-	  r[i] -= JFD[Ic][m] * dv[col];
+          *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
 	else
 	{
-	  if (col == Im)
-	    a[i] = JFD[Ic][m];
+	  if (col == Im && i > 0)
+            gsl_vector_set(f,i-1,JFD[Ic][m]);
 	  if (col == Ic)
-	    b[i] = JFD[Ic][m];
-	  if (col == Ip)
-	    c[i] = JFD[Ic][m];
+            gsl_vector_set(diag,i,JFD[Ic][m]);
+	  if (col == Ip && i < n1-1)
+            gsl_vector_set(e,i,JFD[Ic][m]);
 	}
       }
     }
-    tridag (a, b, c, r, u, n1);
+    gsl_linalg_solve_tridiag(diag, e, f, b, x);
     for (i = 0; i < n1; i++)
     {
       Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
-      dv[Ic] = u[i];
+      dv[Ic] = gsl_vector_get(x, i);
     }
   }
-  free_dvector (a, 0, n1 - 1);
-  free_dvector (b, 0, n1 - 1);
-  free_dvector (c, 0, n1 - 1);
-  free_dvector (r, 0, n1 - 1);
-  free_dvector (u, 0, n1 - 1);
+
+  gsl_vector_free(diag);
+  gsl_vector_free(e);
+  gsl_vector_free(f);
+  gsl_vector_free(b);
+  gsl_vector_free(x);
 }
 
 /* --------------------------------------------------------------------------*/
@@ -158,52 +161,52 @@ LineRelax_be (CCTK_REAL * restrict const dv,
               CCTK_REAL const * restrict const * restrict const JFD)
 {
   int j, m, Ic, Ip, Im, col, ivar;
-  CCTK_REAL *a, *b, *c, *r, *u;
-  a = dvector (0, n2 - 1);
-  b = dvector (0, n2 - 1);
-  c = dvector (0, n2 - 1);
-  r = dvector (0, n2 - 1);
-  u = dvector (0, n2 - 1);
+
+  gsl_vector *diag = gsl_vector_alloc(n2);
+  gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
+  gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
+  gsl_vector *b = gsl_vector_alloc(n2);   /* rhs */
+  gsl_vector *x = gsl_vector_alloc(n2);   /* solution vector */
 
   for (ivar = 0; ivar < nvar; ivar++)
   {
+    gsl_vector_set_zero(diag);
+    gsl_vector_set_zero(e);
+    gsl_vector_set_zero(f);
     for (j = 0; j < n2; j++)
     {
       Ip = Index (ivar, i, j + 1, k, nvar, n1, n2, n3);
       Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
       Im = Index (ivar, i, j - 1, k, nvar, n1, n2, n3);
-      a[j] = 0;
-      b[j] = 0;
-      c[j] = 0;
-      r[j] = rhs[Ic];
+      gsl_vector_set(b,j,rhs[Ic]);
       for (m = 0; m < ncols[Ic]; m++)
       {
 	col = cols[Ic][m];
 	if (col != Ip && col != Ic && col != Im)
-	  r[j] -= JFD[Ic][m] * dv[col];
+          *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
 	else
 	{
-	  if (col == Im)
-	    a[j] = JFD[Ic][m];
+	  if (col == Im && j > 0)
+            gsl_vector_set(f,j-1,JFD[Ic][m]);
 	  if (col == Ic)
-	    b[j] = JFD[Ic][m];
-	  if (col == Ip)
-	    c[j] = JFD[Ic][m];
+            gsl_vector_set(diag,j,JFD[Ic][m]);
+	  if (col == Ip && j < n2-1)
+            gsl_vector_set(e,j,JFD[Ic][m]);
 	}
       }
     }
-    tridag (a, b, c, r, u, n2);
+    gsl_linalg_solve_tridiag(diag, e, f, b, x);
     for (j = 0; j < n2; j++)
     {
       Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
-      dv[Ic] = u[j];
+      dv[Ic] = gsl_vector_get(x, j);
     }
   }
-  free_dvector (a, 0, n2 - 1);
-  free_dvector (b, 0, n2 - 1);
-  free_dvector (c, 0, n2 - 1);
-  free_dvector (r, 0, n2 - 1);
-  free_dvector (u, 0, n2 - 1);
+  gsl_vector_free(diag);
+  gsl_vector_free(e);
+  gsl_vector_free(f);
+  gsl_vector_free(b);
+  gsl_vector_free(x);
 }
 
 /* --------------------------------------------------------------------------*/
diff --git a/src/TP_utilities.c b/src/TP_utilities.c
index 2d44511..7d577c1 100644
--- a/src/TP_utilities.c
+++ b/src/TP_utilities.c
@@ -6,16 +6,7 @@
 #include <stdlib.h>
 #include "TP_utilities.h"
 
-/*---------------------------------------------------------------------------*/
-void
-nrerror (char error_text[])
-/* Numerical Recipes standard error handler */
-{
-  fprintf (stderr, "Numerical Recipes run-time error...\n");
-  fprintf (stderr, "%s\n", error_text);
-  fprintf (stderr, "...now exiting to system...\n");
-  exit (1);
-}
+#include "cctk_Functions.h"
 
 /*---------------------------------------------------------------------------*/
 int *
@@ -26,7 +17,7 @@ ivector (long nl, long nh)
 
   v = (int *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (int)));
   if (!v)
-    nrerror ("allocation failure in ivector()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure in ivector()");
   return v - nl + NR_END;
 }
 
@@ -39,7 +30,7 @@ dvector (long nl, long nh)
 
   v = (CCTK_REAL *) malloc ((size_t) ((nh - nl + 1 + NR_END) * sizeof (CCTK_REAL)));
   if (!v)
-    nrerror ("allocation failure in dvector()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure in dvector()");
   return v - nl + NR_END;
 }
 
@@ -54,7 +45,7 @@ imatrix (long nrl, long nrh, long ncl, long nch)
   /* allocate pointers to rows */
   m = (int **) malloc ((size_t) ((nrow + NR_END) * sizeof (int *)));
   if (!m)
-    nrerror ("allocation failure 1 in matrix()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 1 in matrix()");
   m += NR_END;
   m -= nrl;
 
@@ -62,7 +53,7 @@ imatrix (long nrl, long nrh, long ncl, long nch)
   /* allocate rows and set pointers to them */
   m[nrl] = (int *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (int)));
   if (!m[nrl])
-    nrerror ("allocation failure 2 in matrix()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 2 in matrix()");
   m[nrl] += NR_END;
   m[nrl] -= ncl;
 
@@ -84,7 +75,7 @@ dmatrix (long nrl, long nrh, long ncl, long nch)
   /* allocate pointers to rows */
   m = (CCTK_REAL **) malloc ((size_t) ((nrow + NR_END) * sizeof (CCTK_REAL *)));
   if (!m)
-    nrerror ("allocation failure 1 in matrix()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 1 in matrix()");
   m += NR_END;
   m -= nrl;
 
@@ -92,7 +83,7 @@ dmatrix (long nrl, long nrh, long ncl, long nch)
   m[nrl] =
     (CCTK_REAL *) malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (CCTK_REAL)));
   if (!m[nrl])
-    nrerror ("allocation failure 2 in matrix()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 2 in matrix()");
   m[nrl] += NR_END;
   m[nrl] -= ncl;
 
@@ -114,7 +105,7 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
   /* allocate pointers to pointers to rows */
   t = (CCTK_REAL ***) malloc ((size_t) ((nrow + NR_END) * sizeof (CCTK_REAL **)));
   if (!t)
-    nrerror ("allocation failure 1 in f3tensor()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 1 in f3tensor()");
   t += NR_END;
   t -= nrl;
 
@@ -123,7 +114,7 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
     (CCTK_REAL **)
     malloc ((size_t) ((nrow * ncol + NR_END) * sizeof (CCTK_REAL *)));
   if (!t[nrl])
-    nrerror ("allocation failure 2 in f3tensor()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 2 in f3tensor()");
   t[nrl] += NR_END;
   t[nrl] -= ncl;
 
@@ -132,7 +123,7 @@ d3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
     (CCTK_REAL *)
     malloc ((size_t) ((nrow * ncol * ndep + NR_END) * sizeof (CCTK_REAL)));
   if (!t[nrl][ncl])
-    nrerror ("allocation failure 3 in f3tensor()");
+    CCTK_WARN (CCTK_WARN_ABORT, "allocation failure 3 in f3tensor()");
   t[nrl][ncl] += NR_END;
   t[nrl][ncl] -= ndl;
 
@@ -525,6 +516,7 @@ Ccoth (dcomplex z)
 /*--------------------------------------------------------------------------*/
 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;
@@ -572,6 +564,7 @@ chebft_Zeros (CCTK_REAL u[], int n, int inv)
 
 void
 chebft_Extremes (CCTK_REAL u[], int n, int inv)
+    /* eq. 5.8.7 and 5.8.8 at x = (5.8.5) of 2nd edition C++ NR */
 {
   int k, j, isignum, N = n - 1;
   CCTK_REAL fac, sum, PioN, *c;
@@ -627,23 +620,31 @@ chder (CCTK_REAL *c, CCTK_REAL *cder, int n)
 /* --------------------------------------------------------------------------*/
 CCTK_REAL
 chebev (CCTK_REAL a, CCTK_REAL b, CCTK_REAL c[], int m, CCTK_REAL x)
+    /* eq. 5.8.11 of C++ NR (2nd ed) */
 {
-  CCTK_REAL d = 0.0, dd = 0.0, sv, y, y2;
   int j;
-
-  y2 = 2.0 * (y = (2.0 * x - a - b) / (b - a));
-  for (j = m - 1; j >= 1; j--)
-  {
-    sv = d;
-    d = y2 * d - dd + c[j];
-    dd = sv;
+  CCTK_REAL djp2, djp1, dj; /* d_{j+2}, d_{j+1} and d_j */
+  CCTK_REAL y;
+
+  /* rescale input to lie within [-1,1] */
+  y = 2*(x - 0.5*(b+a))/(b-a);
+
+  dj = djp1 = 0;
+  for(j = m-1 ; j >= 1; j--)
+  { 
+    /* advance the coefficients */
+    djp2 = djp1; 
+    djp1 = dj;
+    dj   = 2*y*djp1 - djp2 + c[j];
   }
-  return y * d - dd + 0.5 * c[0];
+
+  return y*dj - djp1 + 0.5*c[0];
 }
 
 /* --------------------------------------------------------------------------*/
 void
 fourft (CCTK_REAL *u, int N, int inv)
+    /* a (slow) Fourier transform, seems to be just eq. 12.1.6 and 12.1.9 of C++ NR (2nd ed) */
 {
   int l, k, iy, M;
   CCTK_REAL x, x1, fac, Pi_fac, *a, *b;
@@ -754,128 +755,6 @@ fourev (CCTK_REAL *u, int N, CCTK_REAL x)
   return result;
 }
 
-/* --------------------------------------------------------------------------*/
-void
-ludcmp (CCTK_REAL **a, int n, int *indx, CCTK_REAL *d)
-{				/* Version of 'ludcmp' of the numerical recipes for*/
-  /* matrices a[0:n-1][0:n-1]*/
-  int i, imax, j, k;
-  CCTK_REAL big, dum, sum, temp;
-  CCTK_REAL *vv;
-
-  vv = dvector (0, n - 1);
-  *d = 1.0;
-  for (i = 0; i < n; i++)
-  {
-    big = 0.0;
-    for (j = 0; j < n; j++)
-      if ((temp = fabs (a[i][j])) > big)
-	big = temp;
-    if (big == 0.0)
-      nrerror ("Singular matrix in routine ludcmp");
-    vv[i] = 1.0 / big;
-  }
-  for (j = 0; j < n; j++)
-  {
-    for (i = 0; i < j; i++)
-    {
-      sum = a[i][j];
-      for (k = 0; k < i; k++)
-	sum -= a[i][k] * a[k][j];
-      a[i][j] = sum;
-    }
-    big = 0.0;
-    for (i = j; i < n; i++)
-    {
-      sum = a[i][j];
-      for (k = 0; k < j; k++)
-	sum -= a[i][k] * a[k][j];
-      a[i][j] = sum;
-      if ((dum = vv[i] * fabs (sum)) >= big)
-      {
-	big = dum;
-	imax = i;
-      }
-    }
-    if (j != imax)
-    {
-      for (k = 0; k < n; k++)
-      {
-	dum = a[imax][k];
-	a[imax][k] = a[j][k];
-	a[j][k] = dum;
-      }
-      *d = -(*d);
-      vv[imax] = vv[j];
-    }
-    indx[j] = imax;
-    if (a[j][j] == 0.0)
-      a[j][j] = TINY;
-    if (j != n)
-    {
-      dum = 1.0 / (a[j][j]);
-      for (i = j + 1; i < n; i++)
-	a[i][j] *= dum;
-    }
-  }
-  free_dvector (vv, 0, n - 1);
-}
-
-/* --------------------------------------------------------------------------*/
-void
-lubksb (CCTK_REAL **a, int n, int *indx, CCTK_REAL b[])
-{				/* Version of 'lubksb' of the numerical recipes for*/
-  /* matrices a[0:n-1][0:n-1] and vectors b[0:n-1]*/
-
-  int i, ii = 0, ip, j;
-  CCTK_REAL sum;
-
-  for (i = 0; i < n; i++)
-  {
-    ip = indx[i];
-    sum = b[ip];
-    b[ip] = b[i];
-    if (ii)
-      for (j = ii; j <= i - 1; j++)
-	sum -= a[i][j] * b[j];
-    else if (sum)
-      ii = i;
-    b[i] = sum;
-  }
-  for (i = n - 1; i >= 0; i--)
-  {
-    sum = b[i];
-    for (j = i + 1; j <= n; j++)
-      sum -= a[i][j] * b[j];
-    b[i] = sum / a[i][i];
-  }
-}
-
-/* -------------------------------------------------------------------------*/
-void
-tridag (CCTK_REAL a[], CCTK_REAL b[], CCTK_REAL c[], CCTK_REAL r[], CCTK_REAL u[], int n)
-{				/* Version of 'tridag' of the numerical recipes for*/
-  /* vectors a, b, c, r, u with indices in the range [0:n-1]*/
-  int j;
-  CCTK_REAL bet, *gam;
-
-  gam = dvector (0, n - 1);
-  if (b[0] == 0.0)
-    nrerror ("Error 1 in tridag");
-  u[0] = r[0] / (bet = b[0]);
-  for (j = 1; j < n; j++)
-  {
-    gam[j] = c[j - 1] / bet;
-    bet = b[j] - a[j] * gam[j];
-    if (bet == 0.0)
-      nrerror ("Error 2 in tridag");
-    u[j] = (r[j] - a[j] * u[j - 1]) / bet;
-  }
-  for (j = (n - 2); j >= 0; j--)
-    u[j] -= gam[j + 1] * u[j + 1];
-  free_dvector (gam, 0, n - 1);
-}
-
 /* ------------------------------------------------------------------------*/
 CCTK_REAL
 norm1 (CCTK_REAL *v, int n)
@@ -917,44 +796,3 @@ scalarproduct (CCTK_REAL *v, CCTK_REAL *w, int n)
 }
 
 /* -------------------------------------------------------------------------*/
-CCTK_REAL
-plgndr (int l, int m, CCTK_REAL x)
-{
-  void nrerror (char error_text[]);
-  CCTK_REAL fact, pll, pmm, pmmp1, somx2;
-  int i, ll;
-
-  if (m < 0 || m > l || fabs (x) > 1.0)
-    nrerror ("Bad arguments in routine plgndr");
-  pmm = 1.0;
-  if (m > 0)
-  {
-    somx2 = sqrt ((1.0 - x) * (1.0 + x));
-    fact = 1.0;
-    for (i = 1; i <= m; i++)
-    {
-      pmm *= -fact * somx2;
-      fact += 2.0;
-    }
-  }
-  if (l == m)
-    return pmm;
-  else
-  {
-    pmmp1 = x * (2 * m + 1) * pmm;
-    if (l == (m + 1))
-      return pmmp1;
-    else
-    {
-      for (ll = m + 2; ll <= l; ll++)
-      {
-	pll = (x * (2 * ll - 1) * pmmp1 - (ll + m - 1) * pmm) / (ll - m);
-	pmm = pmmp1;
-	pmmp1 = pll;
-      }
-      return pll;
-    }
-  }
-}
-
-/* -------------------------------------------------------------------------*/
diff --git a/src/TP_utilities.h b/src/TP_utilities.h
index 95c2045..7ec7201 100644
--- a/src/TP_utilities.h
+++ b/src/TP_utilities.h
@@ -86,12 +86,6 @@ void fourder2 (CCTK_REAL u[], CCTK_REAL d2u[], int N);
 CCTK_REAL fourev (CCTK_REAL *u, int N, CCTK_REAL x);
 
 
-void ludcmp (CCTK_REAL **a, int n, int *indx, CCTK_REAL *d);
-void lubksb (CCTK_REAL **a, int n, int *indx, CCTK_REAL b[]);
-void tridag (CCTK_REAL a[], CCTK_REAL b[], CCTK_REAL c[], CCTK_REAL r[], CCTK_REAL u[],
-	     int n);
 CCTK_REAL norm1 (CCTK_REAL *v, int n);
 CCTK_REAL norm2 (CCTK_REAL *v, int n);
 CCTK_REAL scalarproduct (CCTK_REAL *v, CCTK_REAL *w, int n);
-
-CCTK_REAL plgndr (int l, int m, CCTK_REAL x);