MINICERN/mathlib/gen/d/radmul.F

   1 *
   2 * $Id
   3 *
   4 * $Log
   5 *
   6 #include "gen/pilot.h"
   7 #if defined(CERNLIB_DOUBLE)
   8       SUBROUTINE RADMUL
   9      1 (F,N,A,B,MINPTS,MAXPTS,EPS,WK,IWK,RESULT,RELERR,NFNEVL,IFAIL)
  10       CHARACTER NAME*(*)
  11       PARAMETER (NAME = 'RADMUL')
  12       CALL MTLPRT(NAME,'D120',
  13      +'not available on this machine - see documentation')
  14       RETURN
  15       END
  16
  17       SUBROUTINE DADMUL
  18      1 (F,N,A,B,MINPTS,MAXPTS,EPS,WK,IWK,RESULT,RELERR,NFNEVL,IFAIL)
  19 #include "gen/imp64.inc"
  20
  21 #else
  22       SUBROUTINE DADMUL
  23      1 (F,N,A,B,MINPTS,MAXPTS,EPS,WK,IWK,RESULT,RELERR,NFNEVL,IFAIL)
  24 #include "gen/imp128.inc"
  25       CHARACTER NAME*(*)
  26       PARAMETER (NAME = 'DADMUL')
  27       CALL MTLPRT(NAME,'D120',
  28      +'not available on this machine - see documentation')
  29       RETURN
  30       END
  31
  32       SUBROUTINE RADMUL
  33      1 (F,N,A,B,MINPTS,MAXPTS,EPS,WK,IWK,RESULT,RELERR,NFNEVL,IFAIL)
  34 #endif
  35
  36       LOGICAL LDV
  37
  38       DIMENSION A(*),B(*),WK(*)
  39       DIMENSION CTR(15),WTH(15),WTHL(15),Z(15)
  40       DIMENSION W(2:15,5),WP(2:15,3)
  41
  42       PARAMETER (R1 = 1, HF = R1/2)
  43
  44       PARAMETER (XL2 =  0.35856 85828 00318 073D0)
  45       PARAMETER (XL4 =  0.94868 32980 50513 796D0)
  46       PARAMETER (XL5 =  0.68824 72016 11685 289D0)
  47
  48       PARAMETER (W2 =  980*R1/6561, W4 = 200*R1/19683)
  49       PARAMETER (WP2 =  245*R1/486, WP4 = 25*R1/729)
  50
  51       DATA (W(N,1),W(N,3),N=2,15)
  52      1/-0.193872885230909911D+00,  0.518213686937966768D-01,
  53      2 -0.555606360818980835D+00,  0.314992633236803330D-01,
  54      3 -0.876695625666819078D+00,  0.111771579535639891D-01,
  55      4 -0.115714067977442459D+01, -0.914494741655235473D-02,
  56      5 -0.139694152314179743D+01, -0.294670527866686986D-01,
  57      6 -0.159609815576893754D+01, -0.497891581567850424D-01,
  58      7 -0.175461057765584494D+01, -0.701112635269013768D-01,
  59      8 -0.187247878880251983D+01, -0.904333688970177241D-01,
  60      9 -0.194970278920896201D+01, -0.110755474267134071D+00,
  61      A -0.198628257887517146D+01, -0.131077579637250419D+00,
  62      B -0.198221815780114818D+01, -0.151399685007366752D+00,
  63      C -0.193750952598689219D+01, -0.171721790377483099D+00,
  64      D -0.185215668343240347D+01, -0.192043895747599447D+00,
  65      E -0.172615963013768225D+01, -0.212366001117715794D+00/
  66
  67       DATA (W(N,5),W(N+1,5),N=2,14,2)
  68      1/ 0.871183254585174982D-01,  0.435591627292587508D-01,
  69      2  0.217795813646293754D-01,  0.108897906823146873D-01,
  70      3  0.544489534115734364D-02,  0.272244767057867193D-02,
  71      4  0.136122383528933596D-02,  0.680611917644667955D-03,
  72      5  0.340305958822333977D-03,  0.170152979411166995D-03,
  73      6  0.850764897055834977D-04,  0.425382448527917472D-04,
  74      7  0.212691224263958736D-04,  0.106345612131979372D-04/
  75
  76       DATA (WP(N,1),WP(N,3),N=2,15)
  77      1/-0.133196159122085045D+01,  0.445816186556927292D-01,
  78      2 -0.229218106995884763D+01, -0.240054869684499309D-01,
  79      3 -0.311522633744855959D+01, -0.925925925925925875D-01,
  80      4 -0.380109739368998611D+01, -0.161179698216735251D+00,
  81      5 -0.434979423868312742D+01, -0.229766803840877915D+00,
  82      6 -0.476131687242798352D+01, -0.298353909465020564D+00,
  83      7 -0.503566529492455417D+01, -0.366941015089163228D+00,
  84      8 -0.517283950617283939D+01, -0.435528120713305891D+00,
  85      9 -0.517283950617283939D+01, -0.504115226337448555D+00,
  86      A -0.503566529492455417D+01, -0.572702331961591218D+00,
  87      B -0.476131687242798352D+01, -0.641289437585733882D+00,
  88      C -0.434979423868312742D+01, -0.709876543209876532D+00,
  89      D -0.380109739368998611D+01, -0.778463648834019195D+00,
  90      E -0.311522633744855959D+01, -0.847050754458161859D+00/
  91
  92       RESULT=0
  93       ABSERR=0
  94       IFAIL=3
  95       IF(N .LT. 2 .OR. N .GT. 15) RETURN
  96       IF(MINPTS .GT. MAXPTS) RETURN
  97
  98       IFNCLS=0
  99       LDV=.FALSE.
 100       TWONDM=2**N
 101       IRGNST=2*N+3
 102       IRLCLS=2**N+2*N*(N+1)+1
 103       ISBRGN=IRGNST
 104       ISBRGS=IRGNST
 105       IF(MAXPTS .LT. IRLCLS) RETURN
 106       DO 10 J = 1,N
 107       CTR(J)=(B(J)+A(J))*HF
 108    10 WTH(J)=(B(J)-A(J))*HF
 109
 110    20 RGNVOL=TWONDM
 111       DO 30 J = 1,N
 112       RGNVOL=RGNVOL*WTH(J)
 113    30 Z(J)=CTR(J)
 114       SUM1=F(N,Z)
 115
 116       DIFMAX=0
 117       SUM2=0
 118       SUM3=0
 119       DO 40 J = 1,N
 120       Z(J)=CTR(J)-XL2*WTH(J)
 121       F2=F(N,Z)
 122       Z(J)=CTR(J)+XL2*WTH(J)
 123       F2=F2+F(N,Z)
 124       WTHL(J)=XL4*WTH(J)
 125       Z(J)=CTR(J)-WTHL(J)
 126       F3=F(N,Z)
 127       Z(J)=CTR(J)+WTHL(J)
 128       F3=F3+F(N,Z)
 129       SUM2=SUM2+F2
 130       SUM3=SUM3+F3
 131       DIF=ABS(7*F2-F3-12*SUM1)
 132       DIFMAX=MAX(DIF,DIFMAX)
 133       IF(DIFMAX .EQ. DIF) IDVAXN=J
 134    40 Z(J)=CTR(J)
 135
 136       SUM4=0
 137       DO 70 J = 2,N
 138       J1=J-1
 139       DO 60 K = J,N
 140       DO 50 L = 1,2
 141       WTHL(J1)=-WTHL(J1)
 142       Z(J1)=CTR(J1)+WTHL(J1)
 143       DO 50 M = 1,2
 144       WTHL(K)=-WTHL(K)
 145       Z(K)=CTR(K)+WTHL(K)
 146    50 SUM4=SUM4+F(N,Z)
 147    60 Z(K)=CTR(K)
 148    70 Z(J1)=CTR(J1)
 149
 150       SUM5=0
 151       DO 80 J = 1,N
 152       WTHL(J)=-XL5*WTH(J)
 153    80 Z(J)=CTR(J)+WTHL(J)
 154    90 SUM5=SUM5+F(N,Z)
 155       DO 100 J = 1,N
 156       WTHL(J)=-WTHL(J)
 157       Z(J)=CTR(J)+WTHL(J)
 158       IF(WTHL(J) .GT. 0) GO TO 90
 159   100 CONTINUE
 160
 161       RGNCMP=RGNVOL*(WP(N,1)*SUM1+WP2*SUM2+WP(N,3)*SUM3+WP4*SUM4)
 162       RGNVAL=W(N,1)*SUM1+W2*SUM2+W(N,3)*SUM3+W4*SUM4+W(N,5)*SUM5
 163       RGNVAL=RGNVOL*RGNVAL
 164       RGNERR=ABS(RGNVAL-RGNCMP)
 165       RESULT=RESULT+RGNVAL
 166       ABSERR=ABSERR+RGNERR
 167       IFNCLS=IFNCLS+IRLCLS
 168
 169       IF(LDV) THEN
 170   110  ISBTMP=2*ISBRGN
 171        IF(ISBTMP .GT. ISBRGS) GO TO 160
 172        IF(ISBTMP .LT. ISBRGS) THEN
 173         ISBTPP=ISBTMP+IRGNST
 174         IF(WK(ISBTMP) .LT. WK(ISBTPP)) ISBTMP=ISBTPP
 175        ENDIF
 176        IF(RGNERR .GE. WK(ISBTMP)) GO TO 160
 177        DO 130 K = 0,IRGNST-1
 178   130  WK(ISBRGN-K)=WK(ISBTMP-K)
 179        ISBRGN=ISBTMP
 180        GO TO 110
 181       ENDIF
 182   140 ISBTMP=(ISBRGN/(2*IRGNST))*IRGNST
 183       IF(ISBTMP .GE. IRGNST .AND. RGNERR .GT. WK(ISBTMP)) THEN
 184        DO 150 K = 0,IRGNST-1
 185   150  WK(ISBRGN-K)=WK(ISBTMP-K)
 186        ISBRGN=ISBTMP
 187        GO TO 140
 188       ENDIF
 189
 190   160 WK(ISBRGN)=RGNERR
 191       WK(ISBRGN-1)=RGNVAL
 192       WK(ISBRGN-2)=IDVAXN
 193       DO 170 J = 1,N
 194       ISBTMP=ISBRGN-2*J-2
 195       WK(ISBTMP+1)=CTR(J)
 196   170 WK(ISBTMP)=WTH(J)
 197       IF(LDV) THEN
 198        LDV=.FALSE.
 199        CTR(IDVAX0)=CTR(IDVAX0)+2*WTH(IDVAX0)
 200        ISBRGS=ISBRGS+IRGNST
 201        ISBRGN=ISBRGS
 202        GO TO 20
 203       ENDIF
 204       RELERR=ABSERR/ABS(RESULT)
 205       IF(ISBRGS+IRGNST .GT. IWK) IFAIL=2
 206       IF(IFNCLS+2*IRLCLS .GT. MAXPTS) IFAIL=1
 207       IF(RELERR .LT. EPS .AND. IFNCLS .GE. MINPTS) IFAIL=0
 208       IF(IFAIL .EQ. 3) THEN
 209        LDV=.TRUE.
 210        ISBRGN=IRGNST
 211        ABSERR=ABSERR-WK(ISBRGN)
 212        RESULT=RESULT-WK(ISBRGN-1)
 213        IDVAX0=WK(ISBRGN-2)
 214        DO 190 J = 1,N
 215        ISBTMP=ISBRGN-2*J-2
 216        CTR(J)=WK(ISBTMP+1)
 217   190  WTH(J)=WK(ISBTMP)
 218        WTH(IDVAX0)=HF*WTH(IDVAX0)
 219        CTR(IDVAX0)=CTR(IDVAX0)-WTH(IDVAX0)
 220        GO TO 20
 221       ENDIF
 222       NFNEVL=IFNCLS
 223       RETURN
 224       END