src/Common/gauss-quad.c

/* Derived from go/gaussq.f at www.netlib.org. The stuff seems to be
 * in the public domain.
 */

/*           this set of routines computes the nodes t(j) and weights */
/*        w(j) for gaussian-type quadrature rules with pre-assigned */
/*        nodes.  these are used when one wishes to approximate */

/*                 integral (from a to b)  f(x) w(x) dx */

/*                              n */
/*        by                   sum w  f(t ) */
/*                             j=1  j    j */

/*        (note w(x) and w(j) have no connection with each other.) */
/*        here w(x) is one of six possible non-negative weight */
/*        functions (listed below), and f(x) is the */
/*        function to be integrated.  gaussian quadrature is particularly */
/*        useful on infinite intervals (with appropriate weight */
/*        functions), since then other techniques often fail. */

/*           associated with each weight function w(x) is a set of */
/*        orthogonal polynomials.  the nodes t(j) are just the zeroes */
/*        of the proper n-th degree polynomial. */

/*     input parameters (all real numbers are in double precision) */

/*        kind     an int between 1 and 6 giving the type of */
/*                 quadrature rule: */

/*        kind = 1:  legendre quadrature, w(x) = 1 on (-1, 1) */
/*        kind = 2:  chebyshev quadrature of the first kind */
/*                   w(x) = 1/sqrt(1 - x*x) on (-1, +1) */
/*        kind = 3:  chebyshev quadrature of the second kind */
/*                   w(x) = sqrt(1 - x*x) on (-1, 1) */
/*        kind = 4:  hermite quadrature, w(x) = exp(-x*x) on */
/*                   (-infinity, +infinity) */
/*        kind = 5:  jacobi quadrature, w(x) = (1-x)**alpha * (1+x)** */
/*                   beta on (-1, 1), alpha, beta .gt. -1. */
/*                   note: kind=2 and 3 are a special case of this. */
/*        kind = 6:  generalized laguerre quadrature, w(x) = exp(-x)* */
/*                   x**alpha on (0, +infinity), alpha .gt. -1 */

/*        n        the number of points used for the quadrature rule */
/*        alpha    real parameter used only for gauss-jacobi and gauss- */
/*                 laguerre quadrature (otherwise use 0.d0). */
/*        beta     real parameter used only for gauss-jacobi quadrature-- */
/*                 (otherwise use 0.d0) */
/*        kpts     (int) normally 0, unless the left or right end- */
/*                 point (or both) of the interval is required to be a */
/*                 node (this is called gauss-radau or gauss-lobatto */
/*                 quadrature).  then kpts is the number of fixed */
/*                 endpoints (1 or 2). */
/*        endpts   real array of length 2.  contains the values of */
/*                 any fixed endpoints, if kpts = 1 or 2. */
/*        b        real scratch array of length n */

/*     output parameters (both double precision arrays of length n) */

/*        t        will contain the desired nodes. */
/*        w        will contain the desired weights w(j). */

/*     underflow may sometimes occur, but is harmless. */

/*     references */
/*        1.  golub, g. h., and welsch, j. h., "calculation of gaussian */
/*            quadrature rules," mathematics of computation 23 (april, */
/*            1969), pp. 221-230. */
/*        2.  golub, g. h., "some modified matrix eigenvalue problems," */
/*            siam review 15 (april, 1973), pp. 318-334 (section 7). */
/*        3.  stroud and secrest, gaussian quadrature formulas, prentice- */
/*            hall, englewood cliffs, n.j., 1966. */

/*        original version 20 jan 1975 from stanford */
/*        modified 21 dec 1983 by eric grosse */
/*          imtql2 => gausq2 */
/*          hex constant => d1mach (from core library) */
/*          compute pi using datan */
/*          removed accuracy claims, description of method */
/*          added single precision version */

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include "alberta_intern.h" /* to keep the proto-type in sync */
#include "alberta.h"

static REAL classical(int kind, int n, REAL alpha, REAL beta, REAL b[], REAL t[]);
static REAL solve(REAL shift, int n, REAL a[/*n*/], REAL b[/*n*/]);
static int gaussq2(int n, REAL d[/*n*/], REAL e[/*n*/], REAL z[/*n*/]);

#if 0
static REAL mypow(REAL a, REAL b)
{
  int e = (int)b;
  REAL res;

  if ((REAL)e == b) {
    res = 1.0;
    while (e--) {
      res *= a;
    }
    return res;
  } else {
    return pow(a, b);
  }
}

REAL mygamma(REAL x)
{
  int xi = (int)x, fac;
  REAL res;

  if ((REAL)xi == x && xi >= 1) {
    res = 1.0;
    fac = 1;
    while (--xi) {
      res *= (REAL)(fac++);
    }
    return res;
  } else {
    return tgamma(x);
  }
}
#endif

void _AI_gauss_quad(int kind, int n, REAL alpha,
		    REAL beta, int kpts, REAL endpts[2],
		    REAL t[/*n*/], REAL w[/*n*/])
{
  /* Local variables */
  REAL b[n]; /* scratch */
  int i;
  REAL t1, gam;
  REAL muzero;

  /* Function Body */
/*<       call class (kind, n, alpha, beta, b, t, muzero) >*/
  muzero = classical(kind, n, alpha, beta, &b[0], &t[0]);

  /* the matrix of coefficients is assumed to be symmetric.  the array
   * t contains the diagonal elements, the array b the off-diagonal
   * elements.  make appropriate changes in the lower right 2 by 2
   * submatrix.
   */
/*<       if (kpts.eq.0)  go to 100 >*/
  switch (kpts) {
  case 0:
    break;
  case 1:
    /* if kpts=1, only t(n) must be changed */
/*<       t(n) = solve(endpts(1), n, t, b)*b(n-1)**2 + endpts(1) >*/
    t[n-1] = solve(endpts[0], n, &t[0], &b[0]) * SQR(b[n-2]) + endpts[0];
/*<       go to 100 >*/
    break;
  case 2:
/*<       if (kpts.eq.2)  go to  50 >*/
    /* if kpts=2, t(n) and b(n-1) must be recomputed */
    /*<    50 gam = solve(endpts(1), n, t, b) >*/
    gam = solve(endpts[0], n, &t[0], &b[0]);
/*<       t1 = ((endpts(1) - endpts(2))/(solve(endpts(2), n, t, b) - gam)) >*/
    t1 = (endpts[0] - endpts[1]) / (solve(endpts[1], n, &t[0], &b[0]) - gam);
/*<       b(n-1) = dsqrt(t1) >*/
    b[n - 1] = sqrt(t1);
/*<       t(n) = endpts(1) + gam*t1 >*/
    t[n-1] = endpts[0] + gam * t1;
    break;
  }
  /* note that the indices of the elements of b run from 0 to n-2 and
   * thus the value of b[n-1] is arbitrary.  now compute the
   * eigenvalues of the symmetric tridiagonal matrix, which has been
   * modified as necessary.  the method used is a ql-type method with
   * origin shifting
   */
/*<   100 w(1) = 1.0d0 >*/
  w[0] = 1.0;
/*<       do 105 i = 2, n >*/
  for (i = 1; i < n; ++i) {
/*<   105    w(i) = 0.0d0 >*/
    w[i] = 0.0;
  }

/*<       call gausq2 (n, t, b, w, ierr) >*/
  gaussq2(n, &t[0], &b[0], &w[0]);
/*<       do 110 i = 1, n >*/
  for (i = 0; i < n; ++i) {
/*<   110    w(i) = muzero * w(i) * w(i) >*/
    w[i] = muzero * SQR(w[i]);
  }
} /* gaussq */

/*<       double precision function solve(shift, n, a, b) >*/
/*       this procedure performs elimination to solve for the */
/*       n-th component of the solution delta to the equation */

/*             (jn - shift*identity) * delta  = en, */

/*       where en is the vector of all zeroes except for 1 in */
/*       the n-th position. */

/*       the matrix jn is symmetric tridiagonal, with diagonal */
/*       elements a(i), off-diagonal elements b(i).  this equation */
/*       must be solved to obtain the appropriate changes in the lower */
/*       2 by 2 submatrix of coefficients for orthogonal polynomials. */
REAL solve(REAL shift, int n, REAL a[/*n*/], REAL b[/*n*/])
{
  REAL ret_val;
  /* Local variables */
  int i;
  REAL alpha;

/*<       double precision shift, a(n), b(n), alpha >*/

/*<       alpha = a(1) - shift >*/
  /* Function Body */
  alpha = a[0] - shift;
/*<       n-1 = n - 1 >*/
/*<       do 10 i = 2, n-1 >*/
  for (i = 1; i < n-1; ++i) {
/*<    10    alpha = a(i) - shift - b(i-1)**2/alpha >*/
/* L10: */
    alpha = a[i] - shift - SQR(b[i-1]) / alpha;
  }
/*<       solve = 1.0d0/alpha >*/
  ret_val = 1.0 / alpha;
  return ret_val;
/*<       end >*/
} /* solve */

/*<       subroutine classical(kind, n, alpha, beta, b, a, muzero) >*/
/*           this procedure supplies the coefficients a(j), b(j) of the */
/*        recurrence relation */

/*             b p (x) = (x - a ) p   (x) - b   p   (x) */
/*              j j            j   j-1       j-1 j-2 */

/*        for the various classical (normalized) orthogonal polynomials, */
/*        and the zero-th moment */

/*             muzero = integral w(x) dx */

/*        of the given polynomial's weight function w(x).  since the */
/*        polynomials are orthonormalized, the tridiagonal matrix is */
/*        guaranteed to be symmetric. */

/*           the input parameter alpha is used only for laguerre and */
/*        jacobi polynomials, and the parameter beta is used only for */
/*        jacobi polynomials.  the laguerre and jacobi polynomials */
/*        require the gamma function. */
REAL classical(int kind, int n, REAL alpha, REAL beta, REAL b[/*n*/], REAL a[/*n*/])
{
  /* Local variables */
  int i;
  REAL ab;
  REAL a2b2, abi, muzero;

  /* Function Body */
/*<       n-1 = n - 1 >*/
/*<       go to (10, 20, 30, 40, 50, 60), kind >*/
  switch (kind) {
  case 1:
    /* kind = 1:  legendre polynomials p(x) on (-1, +1), w(x) = 1. */
/*<    10 muzero = 2.0d0 >*/
    muzero = 2.0;
/*<       do 11 i = 1, n-1 >*/
    for (i = 0; i < n-1; ++i) {
/*<          a(i) = 0.0d0 >*/
      a[i] = 0.0;
/*<          abi = i >*/
      abi = (REAL)(i+1);
/*<    11    b(i) = abi/dsqrt(4*abi*abi - 1.0d0) >*/
/* L11: */
      b[i] = abi / sqrt(4.0 * SQR(abi) - 1.0);
    }
/*<       a(n) = 0.0d0 >*/
    a[n-1] = 0.0;
/*<       return >*/
    break;
  case 2:
    /* kind = 2: chebyshev polynomials of the first kind t(x) on (-1,
     * +1), w(x) = 1 / sqrt(1 - x*x)
     */
/*<    20 muzero = pi >*/
    muzero = M_PI;
/*<       do 21 i = 1, n-1 >*/
    for (i = 0; i < n-1; ++i) {
/*<          a(i) = 0.0d0 >*/
      a[i] = 0.0;
/*<    21    b(i) = 0.5d0 >*/
      b[i] = 0.5;
    }
/*<       b(1) = dsqrt(0.5d0) >*/
    b[0] = M_SQRT1_2;
/*<       a(n) = 0.0d0 >*/
    a[n-1] = 0.0;
    break;
  case 3:
    /* kind = 3: chebyshev polynomials of the second kind u(x) on (-1,
     * +1), w(x) = sqrt(1 - x*x)
     */
/*<    30 muzero = pi/2.0d0 >*/
    muzero = M_PI_2;
/*<       do 31 i = 1, n-1 >*/
    for (i = 0; i < n-1; ++i) {
/*<          a(i) = 0.0d0 >*/
      a[i] = 0.0;
/*<    31    b(i) = 0.5d0 >*/
      b[i] = 0.5;
    }
/*<       a(n) = 0.0d0 >*/
    a[n-1] = 0.0;
    break;
  case 4:
    /* kind = 4: hermite polynomials h(x) on (-infinity, +infinity),
     * w(x) = exp(-x**2)
     */
/*<    40 muzero = dsqrt(pi) >*/
    muzero = sqrt(M_PI);
/*<       do 41 i = 1, n-1 >*/
    for (i = 0; i < n-1; ++i) {
/*<          a(i) = 0.0d0 >*/
      a[i] = 0.0;
/*<    41    b(i) = dsqrt(i/2.0d0) >*/
      b[i] = sqrt((REAL)(i+1) / 2.0);
    }
/*<       a(n) = 0.0d0 >*/
    a[n-1] = 0.0;
    break;
  case 5:
    /* kind = 5: jacobi polynomials p(alpha, beta)(x) on (-1, +1),
     * w(x) = (1-x)**alpha * (1+x)**beta, alpha and beta greater than
     * -1 */
/*<    50 ab = alpha + beta >*/
    ab = alpha + beta;
/*<       abi = 2.0d0 + ab >*/
    abi = ab + 2.0;
/*<     muzero = 2.0d0 ** (ab + 1.0d0) * dgamma(alpha + 1.0d0) * dgamma(
     x beta + 1.0d0) / dgamma(abi)  >*/
    muzero =
      pow(2.0, ab + 1.0) * tgamma(alpha + 1.0) * tgamma(beta + 1.0)/tgamma(abi);
/*<       a(1) = (beta - alpha)/abi >*/
    a[0] = (beta - alpha) / abi;
/*<  b(1) = dsqrt(4.0d0*(1.0d0 + alpha)*(1.0d0 + beta)/((abi + 1.0d0)*
     1  abi*abi)) >*/
    b[0] = sqrt(4.0*(1.0+alpha)*(1.0+beta) / ((abi + 1.0) * SQR(abi)));
/*<       a2b2 = beta*beta - alpha*alpha >*/
    a2b2 = SQR(beta) - SQR(alpha);
/*<       do 51 i = 2, n-1 >*/
    for (i = 1; i < n-1; ++i) {
/*<          abi = 2.0d0*i + ab >*/
      abi = 2.0*(REAL)(i+1) + ab;
/*<          a(i) = a2b2/((abi - 2.0d0)*abi) >*/
      a[i] = a2b2 / ((abi - 2.0) * abi);
/*< 51    b(i) = dsqrt (4.0d0*i*(i + alpha)*(i + beta)*(i + ab)/
     1   ((abi*abi - 1)*abi*abi))   >*/
      b[i] = sqrt(4.0*(REAL)(i+1)
		  *
		  ((REAL)(i+1)+alpha)*((REAL)(i+1) + beta) * ((REAL)(i+1) + ab)
		  /
		  ((SQR(abi) - 1.0) * SQR(abi)));
    }
/*<       abi = 2.0d0*n + ab >*/
    abi = 2.0*(REAL)n + ab;
/*<       a(n) = a2b2/((abi - 2.0d0)*abi) >*/
    a[n-1] = a2b2 / ((abi - 2.0) * abi);
    break;
  case 6:
    /* kind = 6: laguerre polynomials l(alpha)(x) on (0, +infinity),
     * w(x) = exp(-x) * x**alpha, alpha greater than -1.
     */
/*<    60 muzero = dgamma(alpha + 1.0d0) >*/
    muzero = tgamma(alpha + 1.0);
/*<       do 61 i = 1, n-1 >*/
    for (i = 0; i < n-1; ++i) {
/*<          a(i) = 2.0d0*i - 1.0d0 + alpha >*/
      a[i] = 2.0*(REAL)(i+1) - 1.0 + alpha;
/*<    61    b(i) = dsqrt(i*(i + alpha)) >*/
      b[i] = sqrt((REAL)(i+1) * ((REAL)(i+1) + alpha));
    }
/*<       a(n) = 2.0d0*n - 1 + alpha >*/
    a[n-1] = 2.0*(REAL)n - 1.0 + alpha;
    break;
  default:
    muzero = 0.0;
  }
  return muzero;
/*<       end >*/
} /* class_ */

/*     this subroutine is a translation of an algol procedure, */
/*     num. math. 12, 377-383(1968) by martin and wilkinson, */
/*     as modified in num. math. 15, 450(1970) by dubrulle. */
/*     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). */
/*     this is a modified version of the 'eispack' routine imtql2. */

/*     this subroutine finds the eigenvalues and first components of the */
/*     eigenvectors of a symmetric tridiagonal matrix by the implicit ql */
/*     method. */

/*     on input: */

/*        n is the order of the matrix; */

/*        d contains the diagonal elements of the input matrix; */

/*        e contains the subdiagonal elements of the input matrix */
/*          in its first n-1 positions.  e(n) is arbitrary; */

/*        z contains the first row of the identity matrix. */

/*      on output: */

/*        d contains the eigenvalues in ascending order.  if an */
/*          error exit is made, the eigenvalues are correct but */
/*          unordered for indices 1, 2, ..., ierr-1; */

/*        e has been destroyed; */

/*        z contains the first components of the orthonormal eigenvectors */
/*          of the symmetric tridiagonal matrix.  if an error exit is */
/*          made, z contains the eigenvectors associated with the stored */
/*          eigenvalues; */

/*        ierr is set to */
/*          zero       for normal return, */
/*          j          if the j-th eigenvalue has not been */
/*                     determined after 30 iterations. */

/*< ierr = gaussq2(n, d, e, z) >*/
int gaussq2(int n, REAL d[/*n*/], REAL e[/*n*/], REAL z[/*n*/])
{
  FUNCNAME("gaussq2");
  /* Local variables */
  REAL b, c, f, g;
  int i, j, k, l, m;
  REAL p, r, s;
  int ii, ierr;

/*<       int i, j, k, l, m, n, ii, mml, ierr >*/
/*<       real*8 d(n), e(n), z(n), b, c, f, g, p, r, s >*/
/*<       real*8 dsqrt, dabs, dsign, d1mach >*/

/*<       ierr = 0 >*/
  ierr = 0;
/*<       if (n .eq. 1) go to 1001 >*/
  if (n == 1) {
    return ierr;
  }

/*<       e(n) = 0.0d0 >*/
  e[n-1] = 0.0;
/*<       do 240 l = 1, n >*/
  for (l = 0; l < n; ++l) {
/*<          j = 0 >*/
    j = 0;
/*     :::::::::: look for small sub-diagonal element :::::::::: */
    for (;;) {
/*<   105    do 110 m = l, n >*/
      for (m = l; m < n-1; ++m) {
/*<             if (m .eq. n) go to 120 >*/
      /* put into loop limit */
/*< if (dabs(e(m)) .le. machep * (dabs(d(m)) + dabs(d(m+1))))
     x         go to 120   >*/
	if (fabs(e[m]) <= SQR(REAL_EPSILON) * (fabs(d[m]) + fabs(d[m + 1]))) {
	  break;
	}
/*<   110    continue >*/
/* L110: */
      }
/*<   120    p = d(l) >*/
      p = d[l];
/*<          if (m .eq. l) go to 240 >*/
      if (m == l) {
	break; /* out of the "infinite" loop */
      }
/*<          if (j .eq. 30) go to 1000 >*/
      if (j == 30) {
	ERROR_EXIT("Iteration limit %d reached\n", 30);
	return 1; /* iteration limit reached */
      }
/*<          j = j + 1 >*/
      ++j;
/*     :::::::::: form shift :::::::::: */
/*<          g = (d(l+1) - p) / (2.0d0 * e(l)) >*/
      g = (d[l + 1] - p) / (2.0 * e[l]);
/*<          r = dsqrt(g*g+1.0d0) >*/
      r = sqrt(SQR(g) + 1.0);
/*<          g = d(m) - p + e(l) / (g + dsign(r, g)) >*/
      g = d[m] - p + e[l] / (g + (g >= 0 ? fabs(r) : -fabs(r)));
/*<          s = 1.0d0 >*/
      s = 1.0;
/*<          c = 1.0d0 >*/
      c = 1.0;
/*<          p = 0.0d0 >*/
      p = 0.0;
/*<          mml = m - l >*/
/*      mml = m - l;*/
/*     :::::::::: for i=m-1 step -1 until l do -- :::::::::: */
/*<          do 200 ii = 1, mml >*/
      for (i = m-1; i >= l; --i) {
/*<             f = s * e(i) >*/
	f = s * e[i];
/*<             b = c * e(i) >*/
	b = c * e[i];
/*<             if (dabs(f) .lt. dabs(g)) go to 150 >*/
	if (fabs(f) < fabs(g)) {
/*<   150       s = f / g >*/
	  s = f / g;
/*<             r = dsqrt(s*s+1.0d0) >*/
	  r = sqrt(SQR(s) + 1.0);
/*<             e(i+1) = g * r >*/
	  e[i + 1] = g * r;
/*<             c = 1.0d0 / r >*/
	  c = 1.0 / r;
/*<             s = s * c >*/
	  s *= c;
	} else {
/*<             c = g / f >*/
	  c = g / f;
/*<             r = dsqrt(c*c+1.0d0) >*/
	  r = sqrt(SQR(c) + 1.0);
/*<             e(i+1) = f * r >*/
	  e[i + 1] = f * r;
/*<             s = 1.0d0 / r >*/
	  s = 1.0 / r;
/*<             c = c * s >*/
	  c *= s;
/*<             go to 160 >*/
	}
/*<   160       g = d(i+1) - p >*/
	g = d[i + 1] - p;
/*<             r = (d(i) - g) * s + 2.0d0 * c * b >*/
	r = (d[i] - g) * s + 2.0 * c * b;
/*<             p = s * r >*/
	p = s * r;
/*<             d(i+1) = g + p >*/
	d[i + 1] = g + p;
/*<             g = c * r - b >*/
	g = c * r - b;
/*     :::::::::: form first component of vector :::::::::: */
/*<             f = z(i+1) >*/
	f = z[i + 1];
/*<             z(i+1) = s * z(i) + c * f >*/
	z[i + 1] = s * z[i] + c * f;
/*<   200       z(i) = c * z(i) - s * f >*/
/* L200: */
	z[i] = c * z[i] - s * f;
      }

/*<          d(l) = d(l) - p >*/
      d[l] -= p;
/*<          e(l) = g >*/
      e[l] = g;
/*<          e(m) = 0.0d0 >*/
      e[m] = 0.0;
/*<          go to 105 >*/
    } /* for (;;) */
/*<   240 continue >*/
  }

/*     :::::::::: order eigenvalues and eigenvectors :::::::::: */
/*<       do 300 ii = 2, n >*/
  for (ii = 1; ii < n; ++ii) {
/*<          i = ii - 1 >*/
    i = ii - 1;
/*<          k = i >*/
    k = i;
/*<          p = d(i) >*/
    p = d[i];

/*<          do 260 j = ii, n >*/
    for (j = ii; j < n; ++j) {
/*<             if (d(j) .ge. p) go to 260 >*/
      if (d[j] >= p) {
	continue;
      }
/*<             k = j >*/
      k = j;
/*<             p = d(j) >*/
      p = d[j];
/*<   260    continue >*/
    }

/*<          if (k .eq. i) go to 300 >*/
    if (k == i) {
      continue;
    }
/*<          d(k) = d(i) >*/
    d[k] = d[i];
/*<          d(i) = p >*/
    d[i] = p;
/*<          p = z(i) >*/
    p = z[i];
/*<          z(i) = z(k) >*/
    z[i] = z[k];
/*<          z(k) = p >*/
    z[k] = p;
/*<   300 continue >*/
  }

  return 0;
/*     :::::::::: last card of gaussq2 :::::::::: */
/*<       end >*/
} /* gausq2_ */