MODULE Gauss_Weights

!      ALGORITHM 125, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM,
!      VOL. 5, NO. 10, October, 1962, P.  510--511

! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-03  Time: 14:03:32

! Original Fortran translation donated by

! M.Dow@anu.edu.au
! ANUSF,  Australian National University
! Canberra, Australia

! Tidied up to use workspace arrays, dimension arrays to
! required lengths, use 0: for dimensions of A, add comments,
! and NAG tools to layout source

! trh (20/07/97)

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)

CONTAINS


SUBROUTINE weightc(n, q, e, eps, w, x)

! Computes the abscissae x(i) and the weight coefficients w(i) for a
! Gaussian quadrature method
!      \int_0^b w(x)f(x) dx \approx \sum_{i=1}^{n}w_if(x_i) where
! \int_0^b w(x) dx = 1 and w(x)>= 0. The method requires the order n, a
! tolerance eps and the 2n-1 first coefficients of the continued fraction
!      \int_0^b {w(x) \over z-x} = { 1| \over |z} - {q_1 | \over |1} -
!      {e_1 | \over |z} - {q_2 | \over |1} -
!      {e_2 | \over |z} - \cdots
! to be given, the latter in the two arrays q(n) and e(n-1) all components
! of which are automatically positive by virtue of the condition w(x)>= 0.
! The method works as well if the upper bound b is actually infinity
! (note that b does not appear directly as an argument!) or if the density
! function w(x) dx is replaced by da(x) with a monotonically increasing a(x)
! with at least n points of of variation. The tolerance eps should be given
! in accordance to the machines accuracy (preferably by using the value of
! d1mach(4)).  The result is delivered as two arrays w(n) (the weight
! coefficients) and x(n) (the abscissae).  For a description of the method see
! H Rutishauser, ``On a modification of the QD-algorithm with Graeffe-type
! convergence'' [Proceedings of the IFIPS Congress, Munich, 1962].

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: q(n)
REAL (dp), INTENT(IN)      :: e(n-1)
REAL (dp), INTENT(IN OUT)  :: eps
REAL (dp), INTENT(OUT)     :: w(n)
REAL (dp), INTENT(OUT)     :: x(n)

! Allocate workspaces G & A
REAL (dp)  :: g(n), a(0:n,0:7)
!     ..
!     .. Local Scalars ..
REAL (dp)  :: m, p
INTEGER    :: k
!     ..
!     .. Intrinsic Functions ..
! INTRINSIC ABS,EXP,LOG
!     ..

x(1) = q(1) + e(1)
DO k = 2,n
  g(k-1) = e(k-1)*q(k)/x(k-1)
  IF (k == n) THEN
    x(k) = q(k) - g(k-1)
    
  ELSE
    x(k) = q(k) + e(k) - g(k-1)
  END IF
  
  g(k-1) = g(k-1)/x(k)
  w(k-1) = x(k)/x(k-1)
  x(k-1) = LOG(x(k-1))
END DO
x(n) = LOG(x(n))
p = 1
30  DO k = 1,n - 1
  IF (ABS(g(k)*w(k)) > eps) GO TO 40
END DO
GO TO 50

40 CALL qdgraeffe(n, x, g, w, a)
p = 2*p
GO TO 30

! What follows is a peculiar method to compute the w(k) from
! the given ratios g_k = w_{k+1}/w_k such that
!  \sum_{k=1}^n w_k = 1, but the straightforward formulae to do
! this might well produce overflow of exponent

50 w(1) = 1
m = 0
DO k = 1,n - 1
  w(k+1) = w(k)*g(k)
  IF (w(k) > m) m = w(k)
END DO
! WRITE (6,FMT=*) 'm=',m
m = SUM( w(1:n) )
DO k = 1,n
  w(k) = w(k)/m
  x(k) = EXP(x(k)/p)
END DO
RETURN

END SUBROUTINE weightc



SUBROUTINE red(a, f, n)

!  Subroutine RED reduces a heptadiagonal matrix a to tridiagonal form as
!  described in the paper referenced above.

INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN OUT)  :: a(0:n,0:7)
REAL (dp), INTENT(IN)      :: f(n)

!     ..
REAL (dp)  :: c
INTEGER    :: j, k
!     ..
!     .. Intrinsic Functions ..
! INTRINSIC ABS
!     ..
DO k = 1,n - 1
  DO j = k,n - 1
    
    c = -f(j)*a(j,7)/a(j,2)
    a(j,7) = 0
    a(j+1,2) = a(j+1,2) + c*a(j,5)
    a(j,1) = a(j,1) - c*f(j)*a(j,4)
    a(j,6) = a(j,6) - c*a(j+1,1)
    a(j+1,3) = a(j+1,3) - c*a(j+1,6)
  END DO
  DO j = k,n - 1
    c = -f(j)*a(j,4)/a(j,1)
    a(j,4) = 0
    a(j+1,1) = a(j+1,1) + c*a(j,6)
    a(j+1,6) = a(j+1,6) + c*a(j+1,3)
    a(j,5) = a(j,5) - c*a(j+1,2)
    a(j+1,0) = a(j+1,0) - c*a(j+1,5)
  END DO
  DO j = k + 1,n - 1
    c = -a(j,3)/a(j-1,6)
    a(j,3) = 0
    a(j,6) = a(j,6) + c*a(j,1)
    a(j-1,5) = a(j-1,5) - c*f(j)*f(j)*a(j,0)
    a(j,2) = a(j,2) - c*f(j)*f(j)*a(j,5)
    a(j,7) = a(j,7) - c*f(j)*a(j+1,2)
  END DO
  DO j = k + 1,n - 1
    c = -a(j,0)/a(j-1,5)
    a(j,0) = 0
    a(j+1,2) = a(j+1,2) + c*f(j)*a(j,7)
    a(j,5) = a(j,5) + c*a(j,2)
    a(j,1) = a(j,1) - c*f(j)*f(j)*a(j,6)
    a(j,4) = a(j,4) - c*f(j)*a(j+1,1)
  END DO
END DO
RETURN

END SUBROUTINE red



SUBROUTINE qdgraeffe(n, h, g, f, a)

! Subroutine QDGRAEFFE computes for a given continued fraction
!     f(z) = { 1| \over |z} - {q_1 | \over |1} -
!      {e_1 | \over |z} - {q_2 | \over |1} -
!      {e_2 | \over |z} - \cdots - {q_n | \over |1}
! another one, the poles of which are the squares of the poles of f(z)
! However QDGRAEFFE uses not the coefficients q_1 ... q_n and
! e_1 ... e_{n-1} of f(z) but the quotients f_k = q_{k+1}/q_k and
! g_k = e_k/q_{k+1} for k=1,2,...,n-1 and the h_k = ln(abs(q_k)) for
! k=1,2,...,n, and the results are delivered in the same form.  Routine
! QDGRAEFFE can be used independently, but requires subroutine RED

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(OUT)  :: h(n)
REAL (dp), INTENT(OUT)  :: g(n)
REAL (dp), INTENT(OUT)  :: f(n)
REAL (dp), INTENT(OUT)  :: a(0:n,0:7)

!     ..
!     .. Local Scalars ..
INTEGER  :: k
!     ..
!     .. External Subroutines ..
! EXTERNAL red
!     ..

g(n) = 0
f(n) = 0
DO k = 1,n
  a(k-1,4) = 1.0_dp
  a(k-1,5) = 1.0_dp
  a(k,2) = 1 + g(k)*f(k)
  a(k,1) = 1 + g(k)*f(k)
  a(k,6) = g(k)
  a(k,7) = g(k)
  a(k,0) = 0.0_dp
  a(k,3) = 0.0_dp
END DO
a(n,5) = 0.0_dp

! The array a represents the heptadiagonal matrix Q of the paper
! cited above, but with the modifications needed to avoid large numbers
! and with a peculiar arrangement.

CALL red(a, f, n)
DO k = 1,n
  h(k) = 2*h(k) + LOG(ABS(a(k,1)*a(k,2)))
END DO
DO k = 1,n - 1
  f(k) = f(k)*f(k)*a(k+1,2)*a(k+1,1)/ (a(k,1)*a(k,2))
  g(k) = a(k,5)*a(k,6)/ (a(k+1,1)*a(k+1,2))
END DO
RETURN

END SUBROUTINE qdgraeffe

END MODULE Gauss_Weights



PROGRAM main
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-03  Time: 14:01:08
 
! Driver for  CACM Alg 125 Rutishauser
! Author M.Dow@anu.edu.au
!        ANUSF,  Australian National University
!        Canberra, Australia

! Tidied up to use workspace arrays, real parameter values
! and put through nag tools

! trh (20/07/97)

USE Gauss_Weights
IMPLICIT NONE

!     ..
!     .. Parameters ..
INTEGER, PARAMETER    :: nm = 100
REAL (dp), PARAMETER  :: zero=0.0D0, one=1.0D0, two=2.0D0, three=3.0D0
!     .. Local Scalars ..
REAL (dp)  :: a, b, eps, exact, s
INTEGER    :: i, n, p
!     ..
!     .. Local Arrays ..
REAL (dp) :: e(nm), q(nm), w(nm), x(nm)
!     ..
!     .. External Subroutines ..
! EXTERNAL weightcoeff
!     ..

10 WRITE (6, FMT=*) ' Enter n, eps: '
READ (5, FMT=*, END=20) n, eps
IF (eps <= zero) eps = 1.D-15
IF (n >= nm-1) STOP
WRITE (6, FMT=*) 'n=', n, ' eps=', eps

!  q,e for w=1/2 interval (0,2)
! Ref: W.Jones & W.Thron Continued Fractions ...
!      Encyclopedia of maths...vol 11 p24
!      Continued Fraction expansion of log(1+x),
! transformed to Rutishauser form p33

DO i = 2,nm
  q(i) = two*i*i/ (two*i* (two*i-one))
  e(i) = two*i*i/ (two*i* (two*i+one))
END DO
q(1) = one
e(1) = one/three
CALL weightc(n, q, e, eps, w, x)

! Adjust weights, zeros for w=1, interval (-1,1)
DO i = 1, n
  w(i) = two*w(i)
  x(i) = x(i) - one
END DO

DO i = 1, MIN(n,10)
  WRITE (6, FMT=9000) w(i), x(i)
END DO

! Check for x^4
p = 4
a = -one
b = one
exact = (b** (p+1)-a** (p+1))/ (p+1)
s = zero
DO i = 1, n
  s = s + w(i)*x(i)**p
END DO
WRITE (6, FMT=9010) p, exact, s, exact - s
GO TO 10

20 STOP

9000 FORMAT (f20.16, '  ', f20.16)
9010 FORMAT (i3, ' Exact=', f20.16, ' Quadrature=', f20.16, ' Error=', e14.7)
END PROGRAM main