MINICERN/mathlib/gen/d/riwiad.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:16  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE RIWIAD(F)
  11 #include "gen/imp64.inc"
  12 C
  13 C********** COMMUNICATION TO CALLING PROGRAM
  14 C
  15       COMMON/PARAMS/ACC,NDIM,NSUB,ITER
  16       COMMON/ANSWER/VALUE,ERROR
  17       COMMON/STORE/XA(11),XB(11),XC(11),XD(11),MA(11),MB(11),MC(11)
  18       COMMON/STORE1/R(10000),LR
  19       COMMON/OPTION/IPRINT,ICONV,IRESET
  20       COMMON/RANDOM/NSHOTS
  21       COMMON/INTERN/FACTOR,ALFA,BETA,GAMMA,DELTA,LEVEL,NMIN
  22 C*NS  INTEGER XLOC, FNAME
  23       EXTERNAL F
  24 C
  25 #if defined(CERNLIB_DOUBLE)
  26       REAL Q
  27 C     MODIFICATIONS FOR DOUBLE PRECISION
  28       REAL(N)=DBLE(N)
  29 C............
  30 #endif
  31       CALL RIWIBD
  32 C***************CONSTANTS AND COUNTERS
  33       ZERO=1.0E-12
  34       ONE=1.0-ZERO
  35       ITC=0
  36       NCALL=0
  37       SWI=0.0
  38       SRV=0.0
  39       SWISQ=0.0
  40       ANMIN=NMIN
  41       ANDIM=NDIM
  42       ANSUB=NSUB
  43       RNDIM=1.0/ANDIM
  44       ANSHT1=NSHOTS-1
  45       RNSHTS=1.0/NSHOTS
  46       RNSHT1=1.0/(NSHOTS-1)
  47 C
  48 C********** INITIALISE INTERVALS
  49 C
  50       IF(MA(1).EQ.0) GOTO 201
  51       IF(IRESET.EQ.0) GOTO 202
  52   201 K=INT(ANSUB**RNDIM+ZERO)
  53       IF(K.LT.NMIN) GOTO 207
  54       A=1.0/K
  55       L=NDIM*K
  56       DO 204 I=1,L
  57   204 R(I)=A
  58       DO 203 I=1,NDIM
  59   203 MA(I)=K
  60   202 CONTINUE
  61 C
  62 C********** PRINT OUT INITIAL VALUES
  63 C
  64       IF(IPRINT.EQ.0) GOTO 1
  65       WRITE(6,2)
  66     2 FORMAT('1........ RIEMANN INTEGRATION WITH INTERVAL ADJUSTMENT',
  67      *'........')
  68       WRITE(6,3)      ACC,NDIM,NSUB,ITER
  69 3      FORMAT('0******** INPUT PARAM ********'//
  70      +,' NAME OF INTEGRAND...',13X,'F   '/
  71      +,' RELATIVE ACCURACY... ',12X,F7.5/
  72      +,' DIMENSION OF INTEGRAL... ',I15/
  73      +,' MAXIMUM NUMBER OF SUBVOLUMES... ',I8/
  74      +,' MAXIMUM NUMBER OF ITERATIONS... ',I8)
  75       IF(IPRINT.EQ.2) GOTO 1
  76       WRITE(6,4)NSHOTS,LEVEL,FACTOR,IPRINT,
  77      1ICONV,IRESET,ALFA,BETA,GAMMA,DELTA,NMIN
  78     4 FORMAT('0'/'0******** INTERNAL PARAMETERS ********'//
  79      +,' NUMBER OF CALLS PER SUBVOLUME... ',I7/
  80      +,' LEVEL OF CONFIDENCE (PER CENT)...',I7/
  81      +,' CORRESPONDING FACTOR...',F17.2/
  82      +,' PRINT OPTION...',I25/
  83      +,' CONVERGENCE CRITERION...',I16/
  84      +,' RESET OPTION...',I25/
  85      +,' RATE OF CONVERGENCE...',2X,4F4.2/
  86      +,' MINIMUM NUMBER OF INTERVALS...',I10)
  87 C
  88 C********** ENTRY AND REENTRY FOR ITERATIONS
  89 C
  90     1 ITC=ITC+1
  91       J=1
  92       NID114=0
  93       DO 6 I=1,NDIM
  94       MB(I)=NID114
  95       K=MA(I)
  96       J=J*K
  97     6 NID114=NID114+K
  98       NCALL=NCALL+NSHOTS*J
  99       LA=0
 100       LB=NID114
 101       LC=2*NID114
 102       LD=3*NID114
 103       LE=4*NID114
 104       IF(LE.GT.LR) GOTO 107
 105 C
 106 C********** INITIALISE INTEGRATION LOOP
 107 C
 108 C*UL7 L=0
 109       L=0
 110       V=1.0
 111       DO 110 I=1,NDIM
 112       MC(I)=1
 113       K=L+1
 114       L=L+MA(I)
 115       A=R(LA+K)
 116       XA(I)=A
 117       XB(I)=0.0
 118       XC(I)=V
 119       V=V*A
 120       X=0.0
 121       DO 8 J=K,L
 122       R(LB+J)=X
 123       X=X+R(LA+J)
 124       R(LC+J)=0.0
 125     8 R(LD+J)=0.0
 126       R(LA+L)=R(LA+L)+1.0-X
 127   110 CONTINUE
 128 C
 129 C********** INTEGRATION LOOP
 130 C
 131    15 SF1=0.0
 132       SF2=0.0
 133       DO 9 N=1,NSHOTS
 134       DO 10 I=1,NDIM
 135    10 XD(I)=XB(I)+XA(I)*RNDM(I)
 136       Y=F(XD)
 137       SF1=SF1+Y
 138     9 SF2=SF2+Y**2
 139       SF1=SF1*RNSHTS
 140       SF2=SF2*RNSHTS
 141       C=V*SF1
 142       D=V**2*(SF2-SF1**2)*RNSHT1
 143       DO 11 I=1,NDIM
 144       J=MB(I)+MC(I)
 145       R(LC+J)=R(LC+J)+C
 146    11 R(LD+J)=R(LD+J)+D
 147 C
 148 C********** LOOP THE LOOP
 149 C
 150       I=NDIM
 151    13 J=MC(I)
 152       K=MB(I)+1
 153       IF(J.NE.MA(I)) GOTO 12
 154       MC(I)=1
 155       XA(I)=R(LA+K)
 156       XB(I)=0.0
 157       I=I-1
 158       IF(I) 14,14,13
 159    12 MC(I)=J+1
 160       V=XC(I)
 161       K=K+J
 162       A=R(LA+K)
 163       XA(I)=A
 164       XB(I)=R(LB+K)
 165       V=V*A
 166    16 IF(I.EQ.NDIM) GOTO 15
 167       I=I+1
 168       XC(I)=V
 169       V=V*XA(I)
 170       GOTO 16
 171 C
 172 C********** CALCULATE INTEGRAL AND ERROR
 173 C
 174    14 J=MA(1)
 175       SI=0.0
 176       SV=0.0
 177       DO 17 I=1,J
 178       SI=SI+R(LC+I)
 179    17 SV=SV+R(LD+I)
 180       IF(SV.EQ.0.0) GOTO 400
 181       RV=1.0/SV
 182       WI=SI*RV
 183       SRV=SRV+RV
 184       SWI=SWI+WI
 185       VALUE=SWI/SRV
 186       FSV=FACTOR*SQRT(SV)
 187       ERROR=FACTOR/SQRT(SRV)
 188       SWISQ=SWISQ+RV*SI**2
 189       CHISQ=SWISQ-VALUE**2*SRV
 190 C
 191 C********** REDUCE TO RELATIVE QUANTITIES
 192 C
 193       RI=1.0/SI
 194       DO 18 J=1,NID114
 195       R(LC+J)=R(LC+J)*RI
 196    18 R(LD+J)=R(LD+J)*RV
 197 C
 198 C********** PRINT OUT THE RESULTS OF THIS ITERATION
 199 C
 200       IF(IPRINT.EQ.0) GOTO 25
 201       WRITE(6,20)ITC,SI,FSV,VALUE,ERROR,CHISQ
 202      1,NCALL,(MA(I),I=1,NDIM)
 203  20   FORMAT('0'//////' ******** ITERATION NUMBER',I3,' ********'//
 204      +,' VALUE OF INTEGRAL...',E13.5/
 205      +,' ESTIMATED ERROR.....',E13.5//
 206      +,' ACCUMULATED VALUE...',E13.5/
 207      +,' ACCUMULATED ERROR...',E13.5//
 208      +,' CHI SQUARE..........',F13.2//
 209      +,' NUMBER OF CALLS TO F',I13/
 210      +,' INTERVAL STRUCTURE..',5X,(25I4)//)
 211       IF(IPRINT.EQ.2) GOTO 25
 212       DO 19 I=1,NDIM
 213       WRITE(6,21)I
 214    21 FORMAT('0'/'0........ AXIS NUMBER',I3)
 215       K=MB(I)
 216       L=MA(I)+K
 217       K=K+1
 218       WRITE(6,22)(R(LA+J),J=K,L)
 219    22 FORMAT('0INTERVAL CHOICE'/(10E13.2))
 220       WRITE(6,23)(R(LC+J),J=K,L)
 221    23 FORMAT('0RELATIVE CONTRIBUTION TO INTEGRAL'/(10E13.2))
 222       WRITE(6,24)(R(LD+J),J=K,L)
 223    24 FORMAT('0RELATIVE CONTRIBUTION TO ERROR'/(10E13.2))
 224    19 CONTINUE
 225 C
 226 C********** CHECK IF ACCURACY REACHED OR END OF ITERATIONS
 227 C
 228    25 ACC1=ERROR/ABS(VALUE)
 229       IF(ACC1.GT.ACC) GOTO 26
 230    28 IF(IPRINT.NE.0) WRITE(6,27) 'WAS',ITC,VALUE,ERROR
 231       RETURN
 232   26  IF(ITC.EQ.ITER) THEN
 233         IF(IPRINT.NE.0) WRITE(6,27)'NOT',ITC,VALUE,ERROR
 234         RETURN
 235       ENDIF
 236 C
 237 C********** RECOMPUTE INTERVAL SIZES
 238 C
 239       L=0
 240       V=1.0
 241       DO 29 I=1,NDIM
 242       K=L+1
 243       M=MA(I)
 244       L=L+M
 245       X=0.0
 246       Y=0.0
 247       Z=0.0
 248       DO 30 J=K,L
 249       A=R(LA+J)
 250       C=R(LC+J)
 251       D=R(LD+J)
 252       GOTO (301,302,303,304),ICONV
 253   301 CONTINUE
 254       P=(1.0-A+ZERO)**ALFA
 255       Q=(-LOG(D+ZERO)+ZERO)**BETA
 256       GOTO 399
 257   302 CONTINUE
 258       P=A**ALFA
 259       Q=(D+ZERO)**(-BETA)
 260       GOTO 399
 261   303 CONTINUE
 262       P=((-LOG(D+ZERO)+ZERO)/(1.0-D+ZERO))**ALFA
 263       Q=((-LOG(A+ZERO)+ZERO)/(1.0-A+ZERO))**(-BETA)
 264       GOTO 399
 265   304 CONTINUE
 266       P=(1.0-LOG(D+ZERO))**ALFA
 267       Q=(1.0-LOG(A+ZERO))**(-BETA)
 268       GOTO 399
 269   399 CONTINUE
 270       B=A*P*Q
 271       R(LB+J)=B
 272       X=X+B
 273       Y=Y+D/A
 274    30 Z=Z+C**2/A
 275       B=ANSHT1 *Y*SV+(Z-1.0)*SI**2
 276       Y=B**GAMMA*M**DELTA
 277       V=V*Y
 278       XA(I)=Y
 279       X=1.0/X
 280       MC(I)=0
 281       DO 31 J=K,L
 282    31 R(LB+J)=R(LB+J)*X
 283    29 CONTINUE
 284 C
 285 C********** RECOMPUTE INTERVAL NUMBERS
 286 C
 287       L=0
 288    34 Q=(ANSUB/(ANMIN**L*V))**(1.0/(NDIM-L))
 289       V=1.0
 290       K=0
 291       DO 32 I=1,NDIM
 292       M=MC(I)
 293       IF(M.EQ.NMIN) GOTO 32
 294       A=XA(I)*Q
 295       M=INT(A+0.5)
 296       IF(M.LE.NMIN) GOTO 33
 297       XA(I)=A
 298       V=V*A
 299       MC(I)=M
 300       GOTO 32
 301    33 MC(I)=NMIN
 302       K=K+1
 303    32 CONTINUE
 304       L=L+K
 305       IF(K.NE.0) GOTO 34
 306 C
 307 C********** COMBINE SIZES AND NUMBERS INTO NEW INTERVALS
 308 C
 309       NID114=0
 310       DO 35 I=1,NDIM
 311    35 NID114=NID114+MC(I)
 312       IF(5*NID114.GT.LR) GOTO 107
 313       L=0
 314       L1=0
 315       DO 36 I=1,NDIM
 316       K=L+1
 317       K1=L1+1
 318       M=MC(I)
 319       M1=MA(I)
 320       L=L+M
 321       L1=L1+M1
 322       A=REAL(M1)/REAL(M)
 323       A=A*ONE
 324       C=0.0
 325       DO 38 J=K,L
 326       B=C
 327       C=C+A
 328       I1=INT(B)
 329       I2=INT(C)
 330       X=0.0
 331       DO 39 I3=I1,I2
 332    39 X=X+R(LB+K1+I3)
 333    38 R(LE+J)=X-(I2+1-C)*R(LB+K1+I2)-(B-I1)*R(LB+K1+I1)
 334    36 CONTINUE
 335       DO 40 L=1,NID114
 336    40 R(LA+L)=R(LE+L)
 337       DO 41 I=1,NDIM
 338    41 MA(I)=MC(I)
 339       GOTO 1
 340 C
 341 C********** NOT STORAGE ENOUGH
 342 C
 343   107 X1=ANDIM*ANSUB**RNDIM
 344       X2=ANMIN*(ANDIM-1.0)+ANSUB/ANMIN**(ANDIM-1.0)
 345       X1=5.0*X1
 346       X2=5.0*X2
 347       X3=SQRT(X1*X2)
 348       I1=INT(X1)
 349       I2=INT(X2)
 350       I3=INT(X3)
 351       WRITE(6,108)LR,I1,I2,I3
 352   108 FORMAT('0'/' STORAGE LENGTH =',I5,'  IS TOO SMALL'/
 353      +,' MINIMUM NECESSARY -',I7/' MAXIMUM POSSIBLE - ',I7/
 354      +,' PROBABLE VALUE - ',I9)
 355       STOP 1
 356 C
 357 C********** TOO SMALL NUMBER OF SUBVOLUMES
 358 C
 359   207 L=NMIN**NDIM
 360       WRITE(6,206)NSUB,L
 361   206 FORMAT('0 NSUB TOO SMALL...',I7/
 362      +,' MINIMUM NECESSARY...',I7)
 363       STOP 2
 364 C
 365 C ********** ZERO VARIANCE ABORT
 366 C
 367   400 WRITE(6,401)
 368   401 FORMAT('0FUNCTION HAS ZERO VARIANCE')
 369       STOP 3
 370    27 FORMAT('0'/'0********'/
 371      +,' THE DESIRED ACCURACY ',A,' OBTAINED AFTER',I4,' ITERATIONS'//
 372      + ' FINAL VALUE...',E13.5/
 373      + ' FINAL ERROR...',E13.5/'0........')
 374       END