Date: Thu, 02 Feb 95 15:00:37 +0000

This is a transcription of Algorithm 001 from the Collected Algorithms of
the ACM.

The original was in Algol, and I had to substitute subindices by the
more common notation (also valid in Algol) of square brackets.

I also include an example full implementation of the algorithm in C (not
that it is such a lot more younger computer language either) ;-)

In any case, it is a nice algorithm for a parallel implementation... And #1!

				jr

			Jose R. Valverde
		European Bioinformatics Institute
			txomsy@ebi.ac.uk

---------------------------------------------------------------------------
1. QuadI
	R. J. Herbold
	National Bureau of Standards, Washington 25, D. C.
	
comment	QuadI is useful when integration of several func-
	tions of same limits at same time using same
	point rule is desired. The interval (a, b) is di-
	vided into m equal subintervals for an n-point
	quadrature integration. p is the number of func-
	tions to be integrated. W[k] and u[k] are normalized
	weights and abscissas respectively, where
	k = 1, 2, 3, ..., n. u[k] must be in ascending order.
	P(B,j) =: (c) is a procedure which must be sup-
	plied by the programmer. It evaluates (c) the
	function (as indicated by j) for B. I[j] is the result
	of integration for function j.;

procedure QuadI (a, b, m, n, p, w[k], u[k], P(B,j) =:(c)) =: (I[j])
	  begin

QuadI:		h := (b - a) / m
	  for	j := 1 (1) p; I[j] := 0
	  	A := a - h / 2
	  for	i := 1(1)m
L1:	  begin	A := A + h
	  for	k := 1 (1) n
L2:	  begin B := A + (h / 2) * u[k]
	  for	j :=1 (1) p
L3:	  begin P(B,j) =:(c)
	  	I[j] := I[j] + W[k] * c	end L3 ; end L2
	  	  end L
	  for	j := 1 (1) p
	  	I[j] := (h / 2) * I[j]
	  return
	  integer (j, k, i)
	  end	QuadI


------------------------------------------------------------------------

Which I understand to be translated into something like this in C:

--------

/* QuadI:
 *	QuadI is useful when integration of several func-
 *	tions of same limits at same time using same
 *	point rule is desired. The interval (a, b) is di-
 *	vided into m equal subintervals for an n-point
 *	quadrature integration. p is the number of func-
 *	tions to be integrated. W[k] and u[k] are normalized
 *	weights and abscissas respectively, where
 *	k = 1, 2, 3, ..., n. u[k] must be in ascending order.
 *	P(B,j) =: (c) is a procedure which must be sup-
 *	plied by the programmer. It evaluates (c) the
 *	function (as indicated by j) for B. I[j] is the result
 *	of integration for function j.;
 */

/* Note that in this implementation QuadI allocates its result
 * from dynamic memory. Hence the caller should free() the
 * pointer (array) returned by QuadI when finished using it.
 *
 * Implementation notes:
 *	To implement function P, I have indexed all the
 *	functions to be integrated into an array. P only
 *	looks up the function to apply using the index
 *	and returns the result given by the function.
 *
 * Personal note: this is just but a parallel implementation of
 *	an integration method. Given the properties of the
 *	method I don't think much is saved using it in current
 *	uniprocessors. But it is a nice example of parallel
 *	processing and could do well in a multiprocessor or
 *	other parallel computer.
 *
 * Implemented by:
 *	Jose R. Valverde
 *	European Bioinformatics Institute.
 *	Jose.Valverde@ebi.ac.uk			(currently)
 *
 *	JRValverde@enlil.iib.uk			(oldies, still valid)
 *	JRamon@mvax.fmed.uam.es
 *
 * Disclaimer:
 *	I made this out of my wits. The European Bioinformatics Institute
 *	does not guarantee or otherwise support or accept any liability
 *	in any manner whatsoever (nor explicit or implicit or ever) with
 *	respect to this software.
 *
 *	Jose R. Valverde. 1 - February - 1995.
 */

/* Let's have first the following p functions which we want to
 * integrate over the same interval simultaneously: */

double function_1(double);
double function_2(double);
/*...*/
double function_p(double);

/* This is an array of functions that take a double argument and
 * return a double value, in which we store the addresses of the
 * functions.  */

double (*(array_of_functions[])(double)) = {
					function_1,
					function_2,
					/* ... */
					function_p}

double P(B, j)
double B;
int j;
{
	return (*(array_of_functions[j]))( B );
}

double *QuadI(a, b, m, n, p, w, u, P)
double a, b;			/* interval limits */
int m;				/* number of subintervals */
int n;				/* for an n-point quadrature integration */
double *w;			/* normalized weights: array of n values */
double *u;			/* abscissas: array of n values */
P;				/* function that evaluates c, the function */
				/* (as indicated by j) for B */
				
{				/* QuadI */

    double h;		/* subinterval size */
    double *I;		/* Result vector. Should be freed by caller! */
    double A, B, c;
    int i, j, k;
	
    I = (double *) malloc(n * sizeof double);
	
    h = (b - a) / m;	/* compute subinterval size */
    for (j = 0; j < p; j++)
	I[j] = 0.0;	/* Initialize result vector */
    A := a - h / 2;
    for (i = 0; i < m; i++) {		/* L1 */
	A += h;
	for (k = 0; k < n; k++) {	/* L2 */
	    B := A + (h / 2) * u[k];
	    for (j = 0; j < p; j++) {	/* L3 */
		c = P(B, j);
		I[j] = I[j] + w[k] * c;
	    }				/* L3 */
	}				/* L2 */
    }					/* L1 */
    for (j = 0; j < p; j++)
	I[j] = (h / 2) * I[j];

    return I;
}