C ALGORITHM 353, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C VOL. 12, NO. 8, August, 1969, PP.457--458.
#! /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:14 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 TOMS353_PRB tests FSPL2
c
c Modified:
c
c 17 January 2006
c
c Author:
c
c John Burkardt
c
implicit none
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TOMS353_PRB:'
write ( *, '(a)' ) ' Test ACM algorithm 353, which uses'
write ( *, '(a)' ) ' Filon quadrature to approximate integrals'
write ( *, '(a)' ) ' of functions of the form F(X)*COS(M*PI*X)'
write ( *, '(a)' ) ' or F(X)*SIN(M*PI*X).'
call test01
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TOMS353_PRB'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
stop
end
subroutine test01
c***********************************************************************
c
cc TEST01 tests FSER1 with integrands of the form F(X)*COS(M*PI*X).
c
implicit none
real a
real b
real c
real eps
real exact
real f1, f2, f3
external f1
external f2
external f3
integer lc
integer ls
integer m
integer max
real pi
real s
pi = 3.141592653589793E+00
a = 0.0E+00
b = 1.0E+00
eps = 0.00001E+00
max = 8
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'TEST01'
write ( *, '(a)' ) ' Use FSER1 to estimate the integrals of'
write ( *, '(a)' ) ' the form F(X) * COS ( M * PI * X )'
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Use integrand factors:'
write ( *, '(a)' ) ' F(X) = 1, X, X*X.'
write ( *, '(a)' ) ' '
write ( *, '(a,a)' )
& ' M Approximate Exact'
write ( *, '(a)' ) ' '
m = 1
lc = 1
ls = 0
call fser1 ( f1, eps, max, m, c, s, lc, ls )
exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' 1: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f2, eps, max, m, c, s, lc, ls )
exact = ( ( cos ( m * pi * b )
& + m * pi * b * sin ( m * pi * b ) )
& - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) )
& / ( m * pi )**2
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f3, eps, max, m, c, s, lc, ls )
exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b )
& + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) )
& - ( 2.0E+00 * m * pi * a * cos ( m * pi * a )
& + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) )
& / ( m * pi )**3
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X*X: ', m, c, exact
write ( *, '(a)' ) ' '
m = 2
lc = 1
ls = 0
call fser1 ( f1, eps, max, m, c, s, lc, ls )
exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' 1: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f2, eps, max, m, c, s, lc, ls )
exact = ( ( cos ( m * pi * b )
& + m * pi * b * sin ( m * pi * b ) )
& - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) )
& / ( m * pi )**2
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f3, eps, max, m, c, s, lc, ls )
exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b )
& + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) )
& - ( 2.0E+00 * m * pi * a * cos ( m * pi * a )
& + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) )
& / ( m * pi )**3
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X*X: ', m, c, exact
write ( *, '(a)' ) ' '
m = 3
lc = 1
ls = 0
call fser1 ( f1, eps, max, m, c, s, lc, ls )
exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' 1: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f2, eps, max, m, c, s, lc, ls )
exact = ( ( cos ( m * pi * b )
& + m * pi * b * sin ( m * pi * b ) )
& - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) )
& / ( m * pi )**2
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X: ', m, c, exact
lc = 1
ls = 0
call fser1 ( f3, eps, max, m, c, s, lc, ls )
exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b )
& + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) )
& - ( 2.0E+00 * m * pi * a * cos ( m * pi * a )
& + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) )
& / ( m * pi )**3
write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) ' X*X: ', m, c, exact
return
end
function f1 ( x )
c***********************************************************************
c
cc F1 evaluates the integrand factor F(X) = 1.
c
implicit none
real f1
real x
f1 = 1.0E+00
return
end
function f2 ( x )
c***********************************************************************
c
cc F2 evaluates the integrand factor F(X) = X.
c
implicit none
real f2
real x
f2 = x
return
end
function f3 ( x )
c***********************************************************************
c
cc F3 evaluates the integrand factor F(X) = X*X.
c
implicit none
real f3
real x
f3 = x*x
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'
TOMS353_PRB:
Test ACM algorithm 353, which uses
Filon quadrature to approximate integrals
of functions of the form F(X)*COS(M*PI*X)
or F(X)*SIN(M*PI*X).
TEST01
Use FSER1 to estimate the integrals of
the form F(X) * COS ( M * PI * X )
Use integrand factors:
F(X) = 1, X, X*X.
M Approximate Exact
1: 1 -0.428771E-07 -0.278275E-07
X: 1 -0.202642 -0.202642
X*X: 1 -0.202642 -0.202642
1: 2 0.255680E-08 0.278275E-07
X: 2 0.136169E-07 0.271765E-07
X*X: 2 0.506606E-01 0.506606E-01
1: 3 -0.215751E-08 -0.253054E-08
X: 3 -0.225158E-01 -0.225158E-01
X*X: 3 -0.225158E-01 -0.225158E-01
TOMS353_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 FSER1 ( F, EPS, MAX, M, C, S, LC, LS )
IMPLICIT NONE
REAL A
REAL B
REAL B1
REAL B2
REAL B3
REAL B4
REAL B5
REAL C
REAL C1
REAL COSIP
REAL EPS
REAL F
REAL F0
REAL F1
REAL FRAC
REAL G
REAL H
INTEGER I
INTEGER IN
INTEGER LC
INTEGER LLC
INTEGER LLS
INTEGER LN
INTEGER LS
INTEGER M
INTEGER MAX
INTEGER MSWTCH
INTEGER N
INTEGER NST
INTEGER NSTOP
REAL ODCOS
REAL ODSIN
REAL P
REAL PI
REAL PVT1
REAL PVT2
REAL S
REAL S1
REAL S1SQ
REAL SUMCOS
REAL SUMSIN
REAL T
REAL T1
REAL T2
REAL TEMP1
REAL TEMP2
REAL THA
REAL TMAX
REAL TP
REAL TSQ
REAL XI
REAL XM
PI = 3.1415926535898E+00
XM = M
C
C F1 = COS ( M * PI ) TEMPORARY.
C
F1 = 1.0E+00 - 2.0E+00 * ( M - ( M / 2 ) * 2 )
F0 = F ( 0.0E+00 )
F1 = F ( 1.0E+00 ) * F1
C
C 'CIR' WILL BE USED THROUGHOUT THESE COMMENTS TO STAND FOR 'SIN' OR
C 'COS' WHENEVER THOSE TWO SYMBOLS MAY OCCUR.
C NOW DEFINE SUMCIR OF THE ENDPOINTS.
C
SUMCOS = ( F1 + F0 ) * 0.5E+00
SUMSIN = 0.0E+00
B1 = 2.0E+00 / 3.0E+00
C
C TMAX IS THE SWITCH-OVER POINT IN THE ANGLE T.
C OUR ANALYSIS INDICATES THAT TMAX = 1/6 IS THE BEST FOR THE ILIAC II,
C WHICH HAS A 44 BIT FLOATING POINT MANTISSA.
C
TMAX = 0.166E+00
C
C N IS THE NUMBER OF ITERATIONS.
C
N = LOG ( 2.0E+00 * XM ) / 0.693E+00
C
C BOTH TMAX AND N MAY BE CHANGED IF THE MACHINE FOR WHICH THIS
C ROUTINE IS INTENDED HAS GREATER OR LESS ACCURACY THAN ILIAC II.
C IF N IS CHANGED, THEN THE CORRESPONDING CHANGES MUST BE MADE
C IN THE ASSIGNMENTS OF H AND NSTOP.
C
H = 1.0E+00 / REAL( 2**N )
C
C H = 1 / 2**N.
C
NSTOP = 2**N - 1
C
C NSTOP = 2**N - 1.
C
T = H * XM
TP = T * PI
NST = 1
C ***** trh ******
C Initialized T2 to zero following Unassigned trap @ line259
T2 = 0.0E0
C ***** end trh ******
ASSIGN 67 TO MSWTCH
C
C LLC AND LLS ARE USED BY THE ROUTINE IN COMPUTED-GO-TO STATEMENTS.
C AS SOON AS LLS AND LLC HAVE BEEN DEFINED, WE CAN USE LS AND LC
C AS RETURN PARAMETERS (SEE ABOVE).
C
IF ( LS ) 1, 1, 2
1 LLS = 2
GO TO 3
2 LLS = 1
LS = MAX
3 IF ( LC ) 4, 4, 5
4 LLC = 2
GO TO 7
5 LLC = 1
LC = MAX
7 LN = 1
C
C ALL OF THE ABOVE IS EXECUTED ONLY ONCE PER CALL.
C NOW THE ITERATION BEGINS.
C
10 ODCOS = 0.0E+00
ODSIN = 0.0E+00
C
C BEGIN SUMMATION FOR ODCOS AND ODSIN.
C
DO 65 I = 1, NSTOP, NST
XI = I
THA = XI * T
C
C THA * PI IS THE ANGLE USED IN THIS I-TH TERM.
C
C CIR(I*T*PI) IS CALCULATED HERE USING THE IDENTITY
C CIR ( INTEGER MULTIPLE OF PI + FRACTIONAL MULTIPLE OF PI )
C = COS ( INTEGER * PI ) * CIR ( FRAC * PI )
C = ( + 0R - ) * CIR ( FRAC * PI ).
C
FRAC = THA
IN = THA
THA = IN
FRAC = ( FRAC - THA ) * PI
C
C THA IS A FLOATING POINT INTEGER, FRAC IS THE FRACTIONAL PART * PI.
C
COSIP = 1 - 2 * ( IN - 2 * ( IN / 2 ) )
TEMP1 = COSIP * F ( XI * H )
C
C TEMP1 = COS ( INTEGER PART ) * F ( I * H ).
C
GO TO ( 50, 55 ) LLS
50 ODSIN = TEMP1 * SIN ( FRAC ) + ODSIN
55 GO TO ( 60, 65 ) LLC
60 ODCOS = TEMP1 * COS ( FRAC ) + ODCOS
65 CONTINUE
GO TO MSWTCH, ( 67, 70 )
67 NST = 2
C
C NOW HAVE MADE UP FOR THE FIRST 4 ITERATION STEPS, SO RESET THESE
C TWO NUMBERS TO LOOK LIKE THE GENERAL CASE.
C
NSTOP = 2**N
C
C NSTOP = 2**N, IN CASE YOU CHANGE STARTING VALUE OF N.
C
ASSIGN 70 TO MSWTCH
GO TO 92
70 TSQ = TP * TP
IF ( T - TMAX ) 74, 74, 75
C
C 74 IS THE POWER SERIES FOR SMALL T, 75 IS THE CLOSED FORM USED WITH
C LARGER VALUES OF T.
C THE POWER SERIES (WITH 'TN' = TP**N)
C A = (2/45)*T3 - (2/315)*T5 + (2/4725)*T7
C B = (2/3) + (2/15)*T2 - (4/105)*T4 + (2/567)*T6 - (4/22275)*T8
C G = (4/3) - (2/15)*T2 - (1/210)*T4 + (1/11340)*T6
C THE NEXT TERM IN G IS TOO SMALL. IT IS (1/997920)*T8.
C
74 A = TP * TSQ * ( 1.0E+00 - TSQ * ( 1.0E+00 - TSQ / 15.0E+00 )
& / 7.0E+00 ) / 22.5E+00
B2 = B1 * TSQ * 0.2E+00
B3 = B2 * TSQ * 2.0E+00 / 7.0E+00
B4 = B3 * TSQ / 10.8E+00
B5 = B4 * TSQ * 14.0E+00 / 275.0E+00
B = B1 + B2 - B3 + B4 - B5
G = 2.0E+00 * B1 - B2 + B3 / 8.0E+00 - B4 / 40.0E+00
C
C G = 2*B1 - B2 + B3/8 - B4/40 + 5*B5/896. IF YOU WANT THE T8
C TERM INCLUDED IN G.
C
GO TO 80
C
C CLOSED FORM OF THE COEFFICIENTS, WHERE AGAIN 'TN' MEANS TP**N.
C A = 1/TP + COS(TP)*SIN(TP)/T2 - 2*(SIN(TP))**2/T3
C B = 2*((1+(COS(TP))**2)/T2 - 2*SIN(TP)*COS(TP)/T3).
C G = 4*( SIN(TP)/T3 - COS(TP)/T2 ).
C
75 IN = T
TEMP1 = 1 - 2 * ( IN - 2 * ( IN / 2 ) )
TEMP2 = IN
C
C TEMP1 IS COS ( INTEGER PART OF TP), TEMP2 IS FRACTIONAL PART OF TP.
C
TEMP2 = ( T - TEMP2 ) * PI
S1 = TEMP1 * SIN ( TEMP2 )
C
C S1 = SIN(TP).
C
C1 = TEMP1 * COS ( TEMP2 )
C
C C1 = COS(TP).
C
P = S1 * C1
S1SQ = S1 * S1
A = ((( -2.0E+00 * S1SQ / TP ) + P ) / TP + 1.0E+00 ) / TP
B = 2.0E+00 * ( ( -2.0E+00 * P / TP ) + 2.0E+00 - S1SQ ) / TSQ
G = 4.0E+00 * ( S1 / TP - C1 ) / TSQ
80 GO TO ( 81, 85 ) LLS
C
C HAVE CALCULATED THE COEFFICIENTS. NOW READY FOR THE INTEGRATION
C FORMULAS.
C
81 T2 = H * ( A * ( F0 - F1 ) + B * SUMSIN + G * ODSIN )
C
C ENDT1 IS A SUBROUTINE WHICH CHECKS FOR THE CONVERGENCE OF THE
C ITERATIONS. ENDT1 REQUIRES THE PRESENT VALUE TO AGREE WITH THE
C PREVIOUS VALUE TO WITHIN EPS2, WHERE
C EPS2 = ( 1.0 + ABS ( PRESENT VALUE ) * EPS.
C EPS IS SUPPLIED BY THE USER.
C
CALL ENDT1 ( PVT2, T2, EPS, S, LLS, LN )
GO TO ( 85, 84 ), LLS
84 LS = N
85 GO TO ( 86, 90 ), LLC
C
C THIS IS THE COSINE INTEGRAL.
C
86 T1 = H * ( B * SUMCOS + G * ODCOS )
CALL ENDT1 ( PVT1, T1, EPS, C, LLC, LN )
GO TO ( 90, 89 ), LLC
89 LC = N
90 LN = 2
C
C NOW TEST TO SEE IF DONE.
C
IF ( LLC + LLS - 3 ) 92, 92, 100
92 N = N + 1
C
C THIS IS THE BEGINNING OF THE ITERATION.
C
IF ( N - MAX ) 95, 95, 100
95 H = 0.5E+00 * H
T = 0.5E+00 * T
TP = 0.5E+00 * TP
NSTOP = 2 * NSTOP
SUMSIN = SUMSIN + ODSIN
SUMCOS = SUMCOS + ODCOS
GO TO 10
100 S = T2
C = T1
RETURN
END
SUBROUTINE ENDT1 ( PREVQT, QUANT, EPS, VALUE, L1, L2 )
IMPLICIT NONE
REAL EPS
INTEGER L1
INTEGER L2
REAL PREVQT
REAL QUANT
REAL REPS
REAL VALUE
GO TO ( 29, 20 ), L2
20 REPS = EPS * ( 1.0E+00 + ABS ( QUANT ) )
23 IF ( ABS ( PREVQT - QUANT ) - REPS ) 25, 25, 29
25 VALUE = QUANT
L1 = 2
GO TO 30
29 PREVQT = QUANT
30 RETURN
END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
# End of shell archive
exit 0
.