C      ALGORITHM 332, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 11, NO. 6, June, 1968, PP.436--437.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Wed Jan 18 20:30:18 2006
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver
RESULTS= Res

Objs1= driver.o src.o
driver: $(Objs1)
	$(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS)
Res: driver 
	./driver >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c*******************************************************************************
c
cc TOMS332_PRB tests the JACOBI routine.
c
c  Modified:
c
c    04 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS332_PRB:'
      write ( *, '(a)' ) '  Test the ACM TOMS 332 Algorithm Jacobi.'

      call test01

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS332_PRB:'
      write ( *, '(a)' ) '  Normal end of execution.'

      write ( *, '(a)' ) ' '
      stop
      end
      subroutine test01

c*******************************************************************************
c
cc TEST01 tests JACOBI against values stored in JACOBI_POLY_VALUES.
c
c  Modified:
c
c    04 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer a
      double precision alfa
      integer b
      double precision beta
      double precision e
      double precision ed
      double precision fd
      integer flagf
      integer flagfd
      double precision fx
      double precision fx2
      integer n
      integer n_data
      double precision x

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01:'
      write ( *, '(a)' ) '  JACOBI_POLY_VALUES returns exact values of '
      write ( *, '(a)' ) '  the Jacobi polynomial.'
      write ( *, '(a)' ) '  JACOBI computes values of '
      write ( *, '(a)' ) '  the Jacobi polynomial.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '   N       A       B      X          Exact           Computed'
      write ( *, '(a)' ) ' '

      n_data = 0

10    continue

        call jacobi_poly_values ( n_data, n, a, b, x, fx )

        if ( n_data == 0 ) then
          go to 20
        end if

        alfa = dble ( a )
        beta = dble ( b )

        call jacobi ( n, alfa, beta, x, fx2, fd, e, ed, flagf, flagfd )

        write ( *, '(2x,i2,2x,i6,2x,i6,2x,f10.6,2x,g16.8,2x,g16.8)' ) 
     &    n, a, b, x, fx, fx2

      go to 10

20    continue

      return
      end
      subroutine jacobi_poly_values ( n_data, n, a, b, x, fx )

c*******************************************************************************
c
cc JACOBI_POLY_VALUES returns some values of the Jacobi polynomial.
c
c  Discussion:
c
c    In Mathematica, the function can be evaluated by:
c
c      JacobiP[ n, a, b, x ]
c
c  Modified:
c
c    04 January 2006
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz and Irene Stegun,
c    Handbook of Mathematical Functions,
c    US Department of Commerce, 1964.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Wolfram Media / Cambridge University Press, 1999.
c
c  Parameters:
c
c    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
c    first call.  On each call, the routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    Output, integer N, the degree of the polynomial.
c
c    Output, integer A, B, parameters of the function.
c
c    Output, double precision X, the argument of the function.
c
c    Output, double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 26 )

      integer a
      integer a_vec(n_max)
      integer b
      integer b_vec(n_max)
      double precision fx
      double precision fx_vec(n_max)
      integer n
      integer n_data
      integer n_vec(n_max)
      double precision x
      double precision x_vec(n_max)

      save a_vec
      save b_vec
      save fx_vec
      save n_vec
      save x_vec

      data a_vec /
     &    0, 0, 0, 0, 
     &    0, 0, 1, 2, 
     &    3, 4, 5, 0, 
     &    0, 0, 0, 0, 
     &    0, 0, 0, 0, 
     &    0, 0, 0, 0, 
     &    0, 0 /
      data b_vec /
     &   1, 1, 1, 1, 
     &   1, 1, 1, 1, 
     &   1, 1, 1, 2, 
     &   3, 4, 5, 1, 
     &   1, 1, 1, 1, 
     &   1, 1, 1, 1, 
     &   1, 1 /
      data fx_vec /
     &    0.1000000000000000D+01, 0.2500000000000000D+00, 
     &   -0.3750000000000000D+00, -0.4843750000000000D+00, 
     &   -0.1328125000000000D+00, 0.2753906250000000D+00, 
     &   -0.1640625000000000D+00, -0.1174804687500000D+01, 
     &   -0.2361328125000000D+01, -0.2616210937500000D+01, 
     &    0.1171875000000000D+00, 0.4218750000000000D+00, 
     &    0.5048828125000000D+00, 0.5097656250000000D+00, 
     &    0.4306640625000000D+00, -0.6000000000000000D+01, 
     &    0.3862000000000000D-01, 0.8118400000000000D+00, 
     &    0.3666000000000000D-01, -0.4851200000000000D+00, 
     &   -0.3125000000000000D+00, 0.1891200000000000D+00, 
     &    0.4023400000000000D+00, 0.1216000000000000D-01, 
     &   -0.4396200000000000D+00, 0.1000000000000000D+01 /
      data n_vec /
     &    0, 1, 2, 3, 
     &    4, 5, 5, 5, 
     &    5, 5, 5, 5, 
     &    5, 5, 5, 5, 
     &    5, 5, 5, 5, 
     &    5, 5, 5, 5, 
     &    5, 5 /
      data x_vec /
     &    0.5D+00, 0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00, 
     &    0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00, 
     &    0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00,  0.5D+00, 
     &   -1.0D+00, -0.8D+00, -0.6D+00, -0.4D+00, -0.2D+00, 
     &    0.0D+00,  0.2D+00,  0.4D+00,  0.6D+00,  0.8D+00, 
     &    1.0D+00 /

      if ( n_data < 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max < n_data ) then
        n_data = 0
        n = 0
        a = 0
        b = 0
        x = 0.0D+00
        fx = 0.0D+00
      else
        n = n_vec(n_data)
        a = a_vec(n_data)
        b = b_vec(n_data)
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS332_PRB:
  Test the ACM TOMS 332 Algorithm Jacobi.
 
TEST01:
  JACOBI_POLY_VALUES returns exact values of 
  the Jacobi polynomial.
  JACOBI computes values of 
  the Jacobi polynomial.
 
   N       A       B      X          Exact           Computed
 
   0       0       1    0.500000     1.0000000         1.0000000    
   1       0       1    0.500000    0.25000000        0.25000000    
   2       0       1    0.500000   -0.37500000       -0.37500000    
   3       0       1    0.500000   -0.48437500       -0.48437500    
   4       0       1    0.500000   -0.13281250       -0.13281250    
   5       0       1    0.500000    0.27539062        0.27539062    
   5       1       1    0.500000   -0.16406250       -0.16406250    
   5       2       1    0.500000    -1.1748047        -1.1748047    
   5       3       1    0.500000    -2.3613281        -2.3613281    
   5       4       1    0.500000    -2.6162109        -2.6162109    
   5       5       1    0.500000    0.11718750        0.11718750    
   5       0       2    0.500000    0.42187500        0.42187500    
   5       0       3    0.500000    0.50488281        0.50488281    
   5       0       4    0.500000    0.50976562        0.50976562    
   5       0       5    0.500000    0.43066406        0.43066406    
   5       0       1   -1.000000    -6.0000000        -6.0000000    
   5       0       1   -0.800000    0.38620000E-01    0.38620000E-01
   5       0       1   -0.600000    0.81184000        0.81184000    
   5       0       1   -0.400000    0.36660000E-01    0.36660000E-01
   5       0       1   -0.200000   -0.48512000       -0.48512000    
   5       0       1    0.000000   -0.31250000       -0.31250000    
   5       0       1    0.200000    0.18912000        0.18912000    
   5       0       1    0.400000    0.40234000        0.40234000    
   5       0       1    0.600000    0.12160000E-01    0.12160000E-01
   5       0       1    0.800000   -0.43962000       -0.43962000    
   5       0       1    1.000000     1.0000000         1.0000000    
 
TOMS332_PRB:
  Normal end of execution.
 
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      SUBROUTINE JACOBI ( DEGREE, ALFA, BETA, X, F, FD, E, ED,
     &  FLAGF, FLAGFD )
C
      IMPLICIT NONE

      DOUBLE PRECISION A, B, C, D, E
      DOUBLE PRECISION ALF
      DOUBLE PRECISION ALFA
      DOUBLE PRECISION BET
      DOUBLE PRECISION BETA
      INTEGER DEGREE
      DOUBLE PRECISION ED
      DOUBLE PRECISION EG
      DOUBLE PRECISION E1
      DOUBLE PRECISION E2
      DOUBLE PRECISION F
      DOUBLE PRECISION FD
      INTEGER FLAGF
      INTEGER FLAGFD
      DOUBLE PRECISION G
      DOUBLE PRECISION H
      INTEGER I
      INTEGER J
      INTEGER K
      INTEGER M
      INTEGER N
      DOUBLE PRECISION P
      DOUBLE PRECISION PD
      DOUBLE PRECISION Q
      DOUBLE PRECISION QD
      DOUBLE PRECISION S
      DOUBLE PRECISION T1
      DOUBLE PRECISION T2
      DOUBLE PRECISION U
      DOUBLE PRECISION V
      DOUBLE PRECISION W
      DOUBLE PRECISION X
      DOUBLE PRECISION Y

      DIMENSION P(25)
      DIMENSION PD(25)
      DIMENSION Q(25)
      DIMENSION QD(25)
      DIMENSION U(25)
      DIMENSION V(25)
      DIMENSION W(25)

      SAVE A
      SAVE ALF
      SAVE B
      SAVE BET
      SAVE M
      SAVE U
      SAVE V
      SAVE W
      SAVE Y

      DATA ALF / -2.0D+00 /
      DATA BET / -2.0D+00 /
      DATA M / -2 /

      Y = epsilon(0.0d0)
      IF ( DEGREE .EQ. 0 ) GO TO 10

      IF ( ( ALFA .NE. ALF ) 
     & .OR. ( BETA .NE. BET ) ) GO TO 1

      IF ( DEGREE .LE. M ) GO TO 5

      I = M
      K = DEGREE - 1
      M = DEGREE

      IF ( I - 2 ) 2, 3, 3
C
C  CALCULATE THE U(J), V(J), W(J) IN 
C  THE RECURRENCE RELATION
C  P(J) = P(J-1) * ( U(J) + V(J) * X ) - P(J-2) * W(J)
C
1     M = DEGREE
      ALF = ALFA
      BET = BETA
      A = ALF + BET
      B = ALF - BET
      U(1) = B / 2.0D+00
      V(1) = 1.0D+00 + A / 2.0D+00
      W(1) = 0.0D+00

      IF ( DEGREE .EQ. 1 ) GO TO 5

2     U(2) = A * B * ( A + 3.0D+00 ) 
     &  / ( 4.0D+00 * ( A + 2.0D+00 )**2 )
      V(2) = ( A + 3.0D+00 ) * ( A + 4.0D+00 ) 
     &  / ( 4.0D+00 * ( A + 2.0D+00 ) )
      W(2) = ( 1.0D+00 + ALF ) * ( 1.0D+00 + BET ) * ( A + 4.0D+00 )
      W(2) = W(2) / ( 2.0D+00 * ( A + 2.0D+00 )**2 )
      I = 2
      K = DEGREE - 1

3     IF ( ( DEGREE .EQ. 2 )
     &  .OR. ( I .GT. K ) ) GO TO 5

      DO 4 J = I, K
        A = 2 * J + 2
        D = ALF + BET
        A = A + D
        B = D * ( A - 1.0D+00 ) * ( ALF - BET )
        C = J + 1
        C = 2.0D+00 * C * ( A - 2.0D+00 ) * ( C + D )
        U(J+1) = B / C
        D = A * ( A - 1.0D+00 ) * ( A - 2.0D+00 )
        V(J+1) = D / C
        D = J
        A = 2.0D+00 * ( D + ALF ) * ( D + BET ) * A
        W(J+1) = A / C
4     CONTINUE
C
C  FIND THE STARTING VALES FOR J = 1 
C  AND J = 2 FOR USE IN THE RECURSION.
C
5     T1 = V(1) * X
      P(1) = U(1) + T1
      S = Y * MAX ( ABS ( U(1) ), ABS ( T1 ) )
      Q(1) = P(1) + S
      PD(1) = V(1)
      QD(1) = V(1)

      IF ( DEGREE .EQ. 1 ) GO TO 7
      T1 = V(2) * X
      G = U(2) + T1
      EG = Y * MAX ( ABS ( U(2) ), ABS ( T1 ) )
      H = G + EG
      T1 = G * P(1)
      E1 = ABS ( EG * P(1) )
      P(2) = T1 - W(2)
      S = Y * ABS ( W(2) )
      S = MAX ( E1, S )
      Q(2) = H * Q(1) - W(2) + S
      PD(2) = G * PD(1) + V(2) * P(1)
      QD(2) = H * QD(1) + V(2) * Q(1)

      IF ( DEGREE .EQ. 2 ) GO TO 7
C
C  USE THE RECURSION.
C
      DO 6 J = 3, DEGREE
        T2 = V(J) * X
        G = U(J) + T2
        EG = Y * MAX ( ABS ( U(J) ), ABS ( T2 ) )
        H = G + EG
        T1 = G * P(J-1)
        T2 = W(J) * P(J-2)
        E1 = ABS ( EG * P(J-1) )
        E2 = ABS ( T2 ) * Y
        P(J) = T1 - T2
        S = MAX ( E1, E2 )
        Q(J)  = H * Q(J-1)  - W(J) * Q(J-2) + S
        PD(J) = G * PD(J-1) - W(J) * PD(J-2)
        QD(J) = H * QD(J-1) - W(J) * QD(J-2)
        PD(J) = PD(J) + V(J) * P(J-1)
        QD(J) = QD(J) + V(J) * Q(J-1)
6     CONTINUE
C
C  PREPARE THE OUTPUT.
C
7     N = DEGREE
      F = P(N)
      IF ( ABS ( F ) .LT. Y ) GO TO 8
      FLAGF = 0
      E = ABS ( 1.0D+00 - Q(N) / F )
      GO TO 9
8     E = ABS ( F - Q(N) )
      FLAGF = 1
9     FD = PD(N)
      IF ( ABS ( FD ) .LT. Y ) GO TO 11
      FLAGFD = 0
      ED = ABS ( 1.0D+00 - QD(N) / FD )
      GO TO 12
10    F = 1.0D+00
      E = 0.0D+00
      FD = 0.0D+00
      ED = 0.0D+00
      FLAGF = 2
      FLAGFD = 2
      GO TO 12
11    ED = ABS ( FD - QD(N) )
      FLAGFD = 1
12    RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0