program drsilm c>> 2001-05-22 DRSILM Krogh Minor change for making .f90 version. c>> 1994-10-19 DRSILM Krogh Changes to use M77CON c>> 1992-03-04 DRSILM Krogh Initial Code. c--S replaces "?": DR?ILM, ?ILUPM c Demonstration driver for SILUPM. c Given a table of sin(xy), for x = 0, .08, .16, ..., and y = 0, .12, .2 c ... interpolates for (x,y) = (2.7, -.1), (.93, .05), (.765, .87) and c obtains an error estimate. c integer ND1, ND2, NDEG1, NDEG2, NUMTAB, NDIM, NX, NDIMT, NEOPT parameter (ND1=50, ND2=40, NDEG1=8, NDEG2=10, NUMTAB=ND1*ND2, 1 NDIM=2, NX=3, NDIMT=2*NDIM+1, NEOPT=2*(2*NDIM+NDEG1)+1) integer I, J, IOPT(4), NTAB(NDIMT), NDEG(NDIM), LUP(NDIM) real X1(NX), X2(NX), YT(NUMTAB), YT2(ND2, ND1), 1 XT(2*NDIM), X(NDIM), Y, EOPT(NEOPT), H1, H2, ANS parameter (H1 = .08E0, H2 = .12E0) equivalence (YT, YT2) data X1 / 2.7E0, .93E0, .765E0 / data X2 / -.1E0, .05E0, .87E0 / c Set IOPT to get an error estimate. data IOPT / 0, NEOPT, 1, 0 / data NTAB(1), NTAB(2), NTAB(3) / ND1, ND2, 0 / data NDEG, LUP / NDEG1, NDEG2, 3, 3 / c c Compute the XT and YT tables XT(1) = 0.E0 XT(2) = H1 XT(3) = 0.E0 XT(4) = H2 do 20 I = 1, ND1 do 10 J = 1, ND2 YT2(J, I) = sin(real((I-1)*(J-1)) * H1 * H2) 10 continue 20 continue print *, 1 'IOP(1) X1 X2 Y Est. Error True Error' do 50 I = 1, NX X(1) = X1(I) X(2) = X2(I) ANS = sin(X(1)*X(2)) call SILUPM (NDIM,X,Y,NTAB,XT,YT,NDEG,LUP,IOPT,EOPT) print '(I5, 2F9.4, F12.8, 1P,E10.2, E11.2)', 1 IOPT(1), X1(I), X2(I), Y, EOPT(1), Y - ANS 50 continue stop end