From 3c3ec9e436cf8cbe6d8c158ad31387068a3f4400 Mon Sep 17 00:00:00 2001 From: rhaas Date: Wed, 24 Mar 2010 00:48:53 +0000 Subject: 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 --- configuration.ccl | 3 + src/Newton.c | 95 ++++++++++++----------- src/TP_utilities.c | 220 +++++++---------------------------------------------- src/TP_utilities.h | 6 -- 4 files changed, 81 insertions(+), 243 deletions(-) create mode 100644 configuration.ccl 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 #include #include +#include +#include #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 #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); -- cgit v1.2.3