C  Adaptive Simpson integration program with simple tests.
C
C  File:  adsimp.f
C
C  See CPS615 class notes for 14 Sept. 1992 for background
C  This could easily be extended to a form of Romberg integration.
C 
C  For specificity, all arrays are parallel with dimension (1:P).


C ---------------------------------------------------------------------------
C Adaptive Simpson integration

C  Calculate  Integral from a to b F(x) dx  using Simpson integration
C  to error < eps.  The program uses P virtual processors, with P
C  chosen to fit the context.  Each processor handles an interval
C  of width h0 = (b - a)/P .

      real function Simpson(a, b, F, eps)

	integer, parameter:: P = 8192	! 8K
	integer, parameter:: Lim = 1024	! Maximum refinement

	real a, b, eps
	INTERFACE
	  function F(x)
	  integer, parameter:: P = 8192	! 8K
	  real, array(1:P) :: F, x
CMPF      ONDPU x
CMF$      LAYOUT x(,)
	END INTERFACE

	real, array(1:P) :: x0, x1, VA, VC	! Working arrays
	real SA, SB, SC, h0, h, Inew, Iold, err	! Sums, steps, integrals, error
	integer j, k, i(1:P)			! count, depth, position
CMPF    ONDPU x0, x1, VA, VC

	i = [1:P]		! Processor coordinate, so to speak
	h0 = (b - a)/P		! Size of initial panel
	h = h0
	x0 = a + (i-1)*h	! Left side of interval
	x1 = cshift(x0, 1, 1)	! Right side of interval
        x1(P) = b		!    except for entry P
        VA = F(x0) + F(x1)  	! For trapezoidal approximation
	SA = sum(VA)		!    (intentionally inefficient)
	Inew = 0.5*h*SA		! The trapezoidal integral
	err = 1E6 + eps		! Forces one iteration of loop
	SB = 0.0		! Terms with weight 2
        k = 1

C	Adaptive iteration with Simpson's rule
C	SA = sum F(x) for initial points with weights 1 2 2 . . . 2 1
C	SB = sum F(x) which in future get weight 2
C	h0/k is current panel width
C	Inew is most recent approximation to the integral
C	The test on k forces termination if convergence fails.

	do while ( err > eps .and. k < Lim )
	  print 21, k, Inew, err	! For convergence study
 21	  FORMAT (" Step ", I4, "  Inew = ", E14.6, "  err = ", E14.6)
	  Iold = Inew			! For error calculation
	  k = 2*k
	  h = h0/k			! New panel width
	  VC = 0.0			! To accumulate sum
	  do, j = 1, k, 2		! Sum new function values
	    VC = VC + F(x0 + h*j)
	  end do
	  SC = sum(VC)			! Total of new function values
	  Inew = h*(SA + 2.0*SB + 4.0*SC)/3.0
	  SB = SB + SC			! Update SB for next cycle
	  err = abs(Inew - Iold)
	end do

	if ( err > eps ) then		! Report non-convergence
	  print *, " Adaptive integration failed to converge."
	end if

	Simpson = Inew

      end function Simpson


C ---------------------------------------------------------------------------
C  A trignometric test case
C
C   Ftrig(x) = x*(sin(x)**2)
C
C  See Itrig below for indefinite integral

      function Ftrig(x)
	integer, parameter:: P=8192	! 8K
	real, array(1:P) :: Ftrig, x
CMPF    ONDPU x
	Ftrig = x*(sin(x))**2
      end function Ftrig


C  Indefinite integral of Ftrig
C
C                               x**2   x*sin(2x)   cos(2x)
C    Itrig = Integeral  Ftrig = ---- - --------- - -------
C                                4         4          8
C
C  Note that this is a scalar function

      real function Itrig(x)
	real x
	Itrig = (x*(x - sin(2*x)) - 0.5*cos(2*x))*0.25
      end function Itrig


C  Announce the trigonometric function

      subroutine Atrig()
	print '(/ "  = x*(sin(x)**2)")'
      end subroutine Atrig


C ---------------------------------------------------------------------------
C  An algebraic example with a singularity at x = 1
C
C                  x**3
C   Falg = -------------------
C          sqrt((x**2 - 1)**3)
C
C  See Ialg below for indefinite integral

      function Falg(x)
	integer, parameter:: P = 8192	! 8K
	real, array(1:P) :: Falg, x
CMPF    ONDPU x
	Falg = (x**3)/sqrt((x*x - 1)**3)
      end function Falg

C  Indefinite integral of Falg
C
C                                                   1
C    Ialg = Integral Falg = sqrt(x**2 - 1) - --------------
C                                            sqrt(x**2 - 1)
C
C  Note that this is a scalar function

      real function Ialg(x)
	real x, s
	s = sqrt(x*x - 1)
	Ialg = s - 1.0/s
      end function Ialg


C Announce the algebraic function

      subroutine Aalg()
	print '(/ "            x**3")'
	print '("  = -------------------")'
	print '("    sqrt((x**2 - 1)**3)")'
	print '(/ " Use with x > 1 .")'
      end subroutine Aalg


C ---------------------------------------------------------------------------
C  Procedure to conduct interactive experiment with integrand Fx (parallel)
C  whose indefinite integral is Ix (scalar)

      subroutine Test(Fx, Ix)

	EXTERNAL Fx

	INTERFACE
	  real function Ix(x)
	  real x
	END INTERFACE

	INTERFACE
	  real function Simpson(a, b, F, eps)
	  real a, b, eps
	  INTERFACE
	    function F(x)
	    integer, parameter:: P = 8192	! 8K
	    real, array(1:P) :: F, x
CMPF        ONDPU x
	  END INTERFACE
	END INTERFACE


	real a, b, Ir, Ia, err, rel
	character*10 yn

C	Loop to conduct trials

  10	print '(" More trials (y or n)? "$)'
	read *, yn
        if ( (yn(1:1) .ne. 'y') .and. (yn(1:1) .ne. 'Y') ) RETURN
	print '(" a, b, eps > "$)'
	read *, a, b, eps
	Ir = Simpson(a, b, Fx, eps)
	Ia = Ix(b) - Ix(a)
	err = Ir - Ia
	rel = abs(err/Ia)
	print 22, Ia, Ir, err, rel
  22	FORMAT (" Analytic = ", E14.6, "  Calculated = ", E14.6 /
     *          " Abs. err = ", E14.6, "  Relativerr = ", E14.6)
	 
	go to 10 

      end subroutine Test


C ---------------------------------------------------------------------------
C Program to try trigonometric and algebraic experiments

      program TrySimpson

	EXTERNAL Ftrig

	EXTERNAL Itrig

	EXTERNAL Falg

	EXTERNAL Ialg

	INTERFACE
	  subroutine Test(Fx, Ix)
	  INTERFACE
	    function Fx(x)
	    integer, parameter:: P = 8192	! 8K
	    real, array(1:P) :: Fx, x
CMPF        ONDPU x
	  END INTERFACE
	  INTERFACE
	    real function Ix(x)
	    real x
	  END INTERFACE
	END INTERFACE

	print '(" Adaptive Simpson integration experiments" // " integrand")'
	call Atrig
	call Test(Ftrig, Itrig)
	print '(// " integrand")'
	call Aalg
	call Test(Falg,  Ialg)

      end program TrySimpson
