c program DRSWCOM2
c>> 1996-05-28 DRSWCOM2 Krogh Removed implicit statement.
c>> 1996-05-28 DRSWCOM2 Krogh Added external statement.
c>> 1995-10-24 DRSWCOM2 Krogh Changed misleading message.
c>> 1994-11-02 DRSWCOM2 Krogh Changes to use M77CON
c>> 1992-05-15 DRSWCOM2 CLL Removed "stop '... finished'"
c>> 1992-03-24 CLL Add parameter MXCASE.
c>> 1992-03-12 CLL
c>> 1987-10-28 Original time stamp
c Demo driver for the SWCOMP package. This code was adapted from the
c test driver to reduce the number of tests and the amount of output.
c Here we have NCASES = 1 and NN(1) = 6, whereas the test driver had
c NCASES = 7 and NN(1) = 0.
c The package SWCOMP computes values and derivatives using the
c approach described by Wengert.
c C. L. Lawson, JPL, 1971. Revised for Fortran77, Jan 1987.
c CLL 8/18/87. Added SWASIN & SWACOS to package.
c CLL 8/31/87. Added SWRCHN and changed order of args in SWCHN.
c CLL 9/1/87. Added SWSINH, SWCOSH, SWTANH.
c 1992-05-15 CLL Removed "stop '... finished'" to simplify
c comparison of output from different systems.
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c--S replaces "?": DR?WCOM2, ?WCOMP
c--& ?COPY, ?MXDIF, ?PASCL, ?VECP, ?WACOS, ?WASIN
c--& ?WATAN, ?WATN2, ?WCHN, ?WCOS, ?WCOSH, ?WDIF, ?WDIF1, ?WEXP
c--& ?WLOG, ?WPRO, ?WPRO1, ?WPWRI, ?WQUO, ?WQUO1, ?WRCHN, ?WSET
c--& ?WSIN, ?WSINH, ?WSQRT, ?WSUM, ?WSUM1, ?WTAN, ?WTANH
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
integer I, IA, ICASE, ICOUNT, II, J, MXCASE, N, NCASES, NMAX,
1 NDIM, NP1, NW, NA, NX
parameter(MXCASE = 7)
parameter(NCASES = 1, NMAX = 6, NDIM = NMAX+1, NW = 20)
parameter(NA = 4, NX = 10)
external SMXDIF
real SMXDIF
real S(295),W(NDIM, NW), TEMP(NDIM), X(NDIM,NX)
real SAVE1(NDIM), SAVE2(NDIM), XEXP(NDIM), XLOG(NDIM)
real A(NA), FIX3(NA)
real ZERO, ONE, TWO, HALF, FOUR, PI, C9, C8, C7
parameter(ZERO = 0.0e0, ONE = 1.0e0, TWO = 2.0e0, HALF = 0.50e0)
parameter(FOUR = 4.0e0, C9 = 0.9e0, C8 = 0.8e0, C7 = 0.7e0)
logical DETAIL, AS(NA), AC(NA)
integer NN(MXCASE), IPWR
data NN/6,1,2,3,4,5,6/
data A / -2.0e0, -1.0e0, 0.75e0, 2.5e0/
data FIX3 / 4*0.0e0 /
data AS / .false., .true., .true., .false. /
data AC / .false., .false., .true., .true. /
data DETAIL / .false. /
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
print'(a)',
* ' DRSWCOM2.. Demo driver for the whole SWCOMP package.',
* ' Will print the numerical error in various calculations.'
PI = ATAN(ONE) * FOUR
FIX3(1) = PI
FIX3(4) = -PI
c
if(DETAIL) then
do 20 I=1,6
do 10 J=1,40
S(J)=ZERO
10 continue
call SPASCL(I-1,S)
II=((I-1)*I)/2 + 1
call SVECP (S,II, '0DPASCL')
20 continue
endif
c Set X(*,1) and X(*,2)
X(1,1) = C8
X(1,2) = C7
do 30 I = 2,NDIM
X(I,1) = C8 * X(I-1,1)
X(I,2) = -C7 * X(I-1,2)
30 continue
c Loop through different values of N
do 500 ICASE= 1, NCASES
N=NN(ICASE)
write(*,'(1x/1x,a,i3/1X)') '>>>>>> Tests with Nderiv = ', N
NP1=N+1
c Set W(*,1) and W(*,2)
W(1,1) = C9
W(1,2) = C9
W(2,1) = C9*C9
W(2,2) = -C9
do 40 I = 3,NDIM
W(I,1) = -C9*W(I-1,1)
W(I,2) = ZERO
40 continue
c Test SWSET
call SWSET (N, C9, -C9, W(1,3))
write(*,'(1x,a,g11.3)')
* 'Error in SET =',
* SMXDIF(NP1,W(1,2),1, W(1,3),1)
c
c Test SWSUM, SWSUM1, SWDIF, SWDIF1, SWPRO1
c
call SWSUM1(N, FOUR, X(1,1), X(1,3))
call SWSUM (N, X(1,3), X(1,2), X(1,3))
call SWDIF (N, X(1,3), X(1,1), X(1,3))
call SWDIF1(N, FOUR, X(1,3), X(1,3))
call SWPRO1(N, -ONE, X(1,3), X(1,3))
write(*,'(1x,a,g11.3)')
* 'Error in -(4-(((4+x)+y))-x) - y =',
* SMXDIF(NP1,X(1,3),1, X(1,2),1)
c
c Test SWQUO1, SWPRO1, and SWPWRI
c
call SWQUO1(N,TWO, X(1,1), X(1,3))
call SWPWRI(N,-1, X(1,3), X(1,4))
call SWPRO1(N,TWO, X(1,4), X(1,5))
write(*,'(1x,a,g11.3)')
* 'Error in 2*((2/x)**-1) - x =',
* SMXDIF(NP1,X(1,5),1, X(1,1),1)
c Test SWPRO and SWQUO
call SWPRO(N, X(1,1), X(1,2), X(1,3))
call SWQUO(N, X(1,3), X(1,2), X(1,4))
write(*,'(1x,a,g11.3)')
* 'Error in (X1*X2)/X2 - X1 =',
* SMXDIF(NP1,X(1,4),1, X(1,1),1)
if(DETAIL) then
print*,'X1, X2, X3 = X1 * X2, X4 = X3/X2'
do 50 J = 1,4
write(*,'(1x,I3,4g15.7/(4x,4g15.7))') J,(X(I,J),I=1,NP1)
50 continue
endif
call SWQUO(N, X(1,3), X(1,2), X(1,3))
write(*,'(1x,a,g11.3)')
* 'Error in (X1*X2)/X2 - X1 =',
* SMXDIF(NP1,X(1,3),1, X(1,1),1)
if(DETAIL) then
print*,' '
print*,'X3 = X3/X2'
write(*,'(1x,I3,4g15.7)') 3,(X(I,3),I=1,NDIM)
endif
c
c Test SWLOG and SWEXP, uses SWCHN
c
call SWLOG (N,W(1,1),W(1,2))
call SWEXP (N,W(1,2),W(1,3))
write(*,'(1x,a,g11.3)')
* 'Error in Exp(Log(2)) - 2 =',
* SMXDIF(NP1,W(1,3),1, W(1,1),1)
if(DETAIL) then
call SVECP(W(1,1),NP1, '0W1 = T = 2.')
call SVECP (W(1,2),NP1, '0W2 = LOG(W1)')
call SVECP (W(1,3),NP1, '0W3 = EXP(W2). SHOULD MATCH W1.')
call SVECP (W(1,4),NP1, '0W4 = 0.+W1. SCALAR ADD.')
endif
c
c Test SWCHN and SWRCHN
call SWSET(N, HALF, ONE, X(1,8))
call SWEXP(N, X(1,8), XEXP)
call SWSET(N, HALF, ONE, X(1,3))
call SWRCHN(N, XEXP, X(1,3))
if(DETAIL) then
print*,' '
print*,'X1, X2 = EXP(X1), X3 = Reversion to Log'
do 60 J = 1,3
write(*,'(1x,I3,4g15.7)') J,(X(I,J),I=1,NDIM)
60 continue
endif
c
call SWSET(N, XEXP(1), ONE, X(1,4))
call SWLOG(N, X(1,4), XLOG)
call SWSET(N, X(1,4), ONE, X(1,6))
call SWRCHN(N, XLOG, X(1,6))
if(DETAIL) then
print*,' '
print*,'X4, X5 = LOG(X4), X6 = Reversion to Exp'
do 70 J = 4,6
write(*,'(1x,I3,4g15.7)') J,(X(I,J),I=1,NDIM)
70 continue
endif
write(*,'(1x,a,g11.3)')
* 'Error in Log(x) - Rev. of Exp(x) =',
* SMXDIF(NP1,XLOG,1, X(1,3),1)
write(*,'(1x,a,g11.3)')
* 'Error in Exp(x) - Rev. of Log(x) =',
* SMXDIF(NP1,XEXP,1, X(1,6),1)
c
call SCOPY(NP1, X(1,1),1, X(1,7),1)
call SCOPY(NP1, X(1,7),1, SAVE1,1)
call SWCHN (N, XEXP, X(1,7))
call SCOPY(NP1, X(1,7),1, SAVE2,1)
call SWRCHN(N, XEXP, X(1,7))
write(*,'(1x,a,g11.3)')
* 'Error in X7;CHN(XEXP,X7);RCHN(XEXP,X7) =',
* SMXDIF(NP1,SAVE1,1, X(1,7),1)
if(DETAIL) then
print*,' '
print*,'X7 ='
J = 7
write(*,'(1x,I3,4g15.7)') J,(SAVE1(I),I=1,NDIM)
print*,'call SWCHN ( N, X2, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(SAVE2(I),I=1,NDIM)
print*,'call SWRCHN ( N, X2, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(X(I,J),I=1,NDIM)
endif
c
c
call SCOPY(NP1, X(1,1),1, X(1,7),1)
call SCOPY(NP1, X(1,7),1, SAVE1,1)
call SWCHN (N, XEXP, X(1,7))
call SCOPY(NP1, X(1,7),1, SAVE2,1)
call SWCHN(N, XLOG, X(1,7))
write(*,'(1x,a,g11.3)')
* 'Error in X7;CHN(XEXP,X7);CHN(XLOG,X7) =',
* SMXDIF(NP1,SAVE1,1, X(1,7),1)
if(DETAIL) then
print*,' '
print*,'X7 ='
J = 7
write(*,'(1x,I3,4g15.7)') J,(SAVE1(I),I=1,NDIM)
print*,'call SWCHN ( N, X2, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(SAVE2(I),I=1,NDIM)
print*,'call SWCHN ( N, X5, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(X(I,J),I=1,NDIM)
endif
c
call SCOPY(NP1, X(1,1),1, X(1,7),1)
call SCOPY(NP1, X(1,7),1, SAVE1,1)
call SWRCHN (N, XEXP, X(1,7))
call SCOPY(NP1, X(1,7),1, SAVE2,1)
call SWRCHN(N, XLOG, X(1,7))
write(*,'(1x,a,g11.3)')
* 'Error in X7;RCHN(XEXP,X7);RCHN(XLOG,X7)=',
* SMXDIF(NP1,SAVE1,1, X(1,7),1)
if(DETAIL) then
print*,' '
print*,'X7 ='
J = 7
write(*,'(1x,I3,4g15.7)') J,(SAVE1(I),I=1,NDIM)
print*,'call SWRCHN ( N, X2, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(SAVE2(I),I=1,NDIM)
print*,'call SWRCHN ( N, X5, X7). X7 ='
write(*,'(1x,I3,4g15.7)') J,(X(I,J),I=1,NDIM)
endif
c Test SWPRO and SWSQRT
c
call SCOPY (NP1, W(1,1),1, W(1,4),1)
do 80 I=1,3
call SWPRO (N,W(1,4),W(1,4),W(1,4))
if(I .eq. 2) call SCOPY(NP1, W(1,4),1, TEMP,1)
if(DETAIL) call SVECP (W(1,4),NP1, '0W4 = W4*W4')
80 continue
c
call SWSQRT (N, W(1,4), W(1,5) )
write(*,'(1x,a,g11.3)')
* 'Error in Sqrt(x*x) - x =',
* SMXDIF(NP1,TEMP,1, W(1,5),1)
if(DETAIL) call SVECP (W(1,5), NP1, '0W5=SQRT(W4)')
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Test of SWPWRI( ,0, , )
c
call SWPWRI(N, 0, X(1,1), X(1,3))
call SWSET (N, ONE, ZERO, X(1,4))
write(*,'(1x,a,g11.3)')
* 'Error in x**0 - 1 =',
* SMXDIF(NP1,X(1,3),1, X(1,4),1)
c
c Loop through tests of SWPWRI
do 100 IPWR = 1,3
write(*,'(1x/1x,a,i3)') 'Test SWPWRI using I = ',IPWR
call SWPWRI(N, IPWR, X(1,1), X(1,3))
call SCOPY(NP1, X(1,1),1, X(1,4),1)
do 90 ICOUNT = 2,IPWR
call SWPRO(N, X(1,1), X(1,4), X(1,4))
90 continue
write(*,'(1x,a,g11.3)')
* 'Error in x**I - x * x * ... * x =',
* SMXDIF(NP1,X(1,3),1, X(1,4),1)
call SWPWRI(N, -IPWR, X(1,1), X(1,5))
call SWPRO(N, X(1,5), X(1,4), X(1,6))
call SWSET (N, ONE, ZERO, X(1,7))
write(*,'(1x,a,g11.3)')
* 'Error in x**(-I) * x**I - 1 =',
* SMXDIF(NP1,X(1,6),1, X(1,7),1)
100 continue
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Loop through quadrants to test trig functions.
c
do 110 IA=1,NA
write(*,'(1x/1x,a,f11.4)') 'Using angle argument,',A(IA)
call SCOPY (NP1, W(1,1),1, W(1,5),1)
W(1,5) = A(IA)
c Test SWSIN and SWASIN
call SWSIN (N,W(1,5),W(1,6))
if(AS(IA)) then
call SWASIN (N,W(1,6),W(1,13))
write(*,'(1x,a,g11.3)')
* 'Error in Asin(Sin(x)) - x =',
* SMXDIF(NP1,W(1,5),1, W(1,13),1)
endif
c Test SWCOS and SWACOS
call SWCOS (N,W(1,5),W(1,7))
if(AC(IA)) then
call SWACOS (N,W(1,7),W(1,14))
write(*,'(1x,a,g11.3)')
* 'Error in Acos(Cos(x)) - x =',
* SMXDIF(NP1,W(1,5),1, W(1,14),1)
endif
c Test SWTAN
call SWTAN(N, W(1,5), W(1,8))
call SWQUO (N,W(1,6),W(1,7),W(1,9))
write(*,'(1x,a,g11.3)')
* 'Error in Tan(x) - Sin(x)/Cos(x) =',
* SMXDIF(NP1,W(1,8),1, W(1,9),1)
c Test SWATN2, SWATAN
call SWATN2(N,W(1,6),W(1,7),W(1,9))
call SWATAN (N,W(1,8),W(1,10))
call SCOPY(NP1,W(1,9),1, TEMP,1)
TEMP(1) = TEMP(1) + FIX3(IA)
write(*,'(1x,a,g11.3)')
* 'Error in FIX + Atan2(y,x) - Atan(y/x) =',
* SMXDIF(NP1,TEMP,1, W(1,10),1)
call SWSET(N,ONE,ZERO,W(1,11) )
call SWATN2(N,W(1,8), W(1,11), W(1,12))
write(*,'(1x,a,g11.3)')
* 'Error in Atan2(y/x, 1) - Atan(y/x) =',
* SMXDIF(NP1,W(1,12),1, W(1,10),1)
if(DETAIL) then
call SVECP (W(1,5),NP1, '0W5 = ANGLE(IA)')
call SVECP (W(1,6),NP1, '0W6 = SIN(W5)')
call SVECP (W(1,13),NP1,
* '0W13 = ASIN(W6). Compare with W5.')
call SVECP (W(1,7),NP1, '0W7 = COS(W5)')
call SVECP (W(1,14),NP1,
* '0W13 = Acos(W7). Compare with W5.')
call SVECP (W(1,8),NP1, '0W8 = W6/W7')
call SVECP (W(1,9),NP1,
* '0W9 = ATAN2(W6,W7). SHOULD MATCH W5')
call SVECP (W(1,10),NP1,
* '0W10 = ATAN(W8). SHOULD MATCH W5 OR DIFFER BY 3.14159265')
call SVECP(W(1,12),NP1, '0W12=ATAN2(W8,1.)' )
endif
110 continue
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c Test SWSINH, SWCOSH, TANH
c First check identity: cosh**2 = sinh**2 + 1
c
call SWSINH(N, X(1,1), X(1,3))
call SWPRO(N, X(1,3), X(1,3), X(1,6))
call SWSUM1(N, ONE, X(1,6), X(1,8))
call SWCOSH(N, X(1,1), X(1,4))
call SWPRO(N, X(1,4), X(1,4), X(1,7))
write(*,'(1x,a,g11.3)')
* 'Error in cosh**2 - sinh**2 - 1 =',
* SMXDIF(NP1,X(1,7),1, X(1,8),1)
c
c Check identity: Tanh = Sinh/Cosh
c
call SWQUO(N, X(1,3), X(1,4), X(1,9))
call SWTANH(N, X(1,1), X(1,5))
write(*,'(1x,a,g11.3)')
* 'Error in tanh - sinh/cosh =',
* SMXDIF(NP1,X(1,5),1, X(1,9),1)
500 continue
end
c ==================================================================
real function SMXDIF(N,X,INCX, Y,INCY)
c
c Compute max norm of difference between the N-vectors x and y.
c The vectors are stored with a storage increment of INCX and INCY
c between successive components.
c C. L. Lawson, JPL, Sept 1987.
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
integer I, INCX, INCY, IX, IY, N
real X(*), Y(*), ZERO, TEMP
parameter(ZERO = 0.0e0)
c -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
TEMP = ZERO
IX = 1 - INCX
IY = 1 - INCY
do 10 I = 1,N
IX = IX + INCX
IY = IY + INCY
TEMP = max(TEMP, abs(X(IX) - Y(IY)))
10 continue
SMXDIF = TEMP
return
end