Eliminating warnings issued by gfortran -Wall
[openmx:openmx.git] / src / sadmvn.f
1       SUBROUTINE RANMVN( N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
2      &                   ABSEPS, RELEPS, ERROR, VALUE, INFORM )
3 *
4 *     A subroutine for computing multivariate normal probabilities.
5 *     This subroutine uses the Monte-Carlo algorithm given in the paper
6 *     "Numerical Computation of Multivariate Normal Probabilities", in
7 *     J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
8 *          Alan Genz
9 *          Department of Mathematics
10 *          Washington State University
11 *          Pullman, WA 99164-3113
12 *          Email : alangenz@wsu.edu
13 *
14 *  Parameters
15 *
16 *     N      INTEGER, the number of variables.
17 *     LOWER  REAL, array of lower integration limits.
18 *     UPPER  REAL, array of upper integration limits.
19 *     INFIN  INTEGER, array of integration limits flags:
20 *            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
21 *            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
22 *            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
23 *            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
24 *     CORREL REAL, array of correlation coefficients; the correlation
25 *            coefficient in row I column J of the correlation matrix
26 *            should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
27 *     MAXPTS INTEGER, maximum number of function values allowed. This
28 *            parameter can be used to limit the time taken. A
29 *            sensible strategy is to start with MAXPTS = 1000*N, and then
30 *            increase MAXPTS if ERROR is too large.
31 *     ABSEPS REAL absolute error tolerance.
32 *     RELEPS REAL relative error tolerance.
33 *     ERROR  REAL estimated absolute error, with 99% confidence level.
34 *     VALUE  REAL estimated value for the integral
35 *     INFORM INTEGER, termination status parameter:
36 *            if INFORM = 0, normal completion with ERROR < EPS;
37 *            if INFORM = 1, completion with ERROR > EPS and MAXPTS
38 *                           function vaules used; increase MAXPTS to
39 *                           decrease ERROR;
40 *            if INFORM = 2, N > 100 or N < 1.
41 *
42       EXTERNAL MVNFNC
43       INTEGER N, INFIN(*), MAXPTS, MPT, INFORM, INFIS, IVLS
44       DOUBLE PRECISION
45      &     CORREL(*), LOWER(*), UPPER(*), MVNFNC,
46      &     ABSEPS, RELEPS, ERROR, VALUE, D, E, EPS, MVNNIT
47       IF ( N .GT. 100 .OR. N .LT. 1 ) THEN
48          INFORM = 2
49          VALUE = 0
50          ERROR = 1
51          RETURN
52       ENDIF
53       INFORM = MVNNIT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E)
54       IF ( N-INFIS .EQ. 0 ) THEN
55          VALUE = 1
56          ERROR = 0
57       ELSE IF ( N-INFIS .EQ. 1 ) THEN
58          VALUE = E - D
59          ERROR = 2E-16
60       ELSE
61 *
62 *        Call then Monte-Carlo integration subroutine
63 *
64          MPT = 25 + 10*N
65          CALL RCRUDE(N-INFIS-1, MPT, MVNFNC, ERROR, VALUE, 0)
66          IVLS = MPT
67  10      EPS = MAX( ABSEPS, RELEPS*ABS(VALUE) )
68          IF ( ERROR .GT. EPS .AND. IVLS .LT. MAXPTS ) THEN
69             MPT = MAX( MIN( INT(MPT*(ERROR/(EPS))**2),
70      &                      MAXPTS-IVLS ), 10 )
71             CALL RCRUDE(N-INFIS-1, MPT, MVNFNC, ERROR, VALUE, 1)
72             IVLS = IVLS + MPT
73             GO TO 10
74          ENDIF
75          IF ( ERROR. GT. EPS .AND. IVLS .GE. MAXPTS ) INFORM = 1
76       ENDIF
77       END
78       SUBROUTINE RCRUDE(NDIM, MAXPTS, FUNCTN, ABSEST, FINEST, IR)
79 *
80 *     Crude Monte-Carlo Algorithm with simple antithetic variates
81 *      and weighted results on restart
82 *
83       EXTERNAL FUNCTN
84       INTEGER NDIM, MAXPTS, M, K, IR, NPTS
85       DOUBLE PRECISION FINEST, ABSEST, X(100), FUN, FUNCTN, UNI,
86      &     VARSQR, VAREST, VARPRD, FINDIF, FINVAL
87       SAVE VAREST
88       IF ( IR .LE. 0 ) THEN
89          VAREST = 0
90          FINEST = 0
91       ENDIF
92       FINVAL = 0
93       VARSQR = 0
94       NPTS = MAXPTS/2
95       DO M = 1,NPTS
96          DO K = 1,NDIM
97             X(K) = UNI()
98          END DO
99          FUN = FUNCTN(NDIM, X)
100          DO K = 1,NDIM
101             X(K) = 1 - X(K)
102          END DO
103          FUN = ( FUNCTN(NDIM, X) + FUN )/2
104          FINDIF = ( FUN - FINVAL )/M
105          VARSQR = ( M - 2 )*VARSQR/M + FINDIF**2
106          FINVAL = FINVAL + FINDIF
107       END DO
108       VARPRD = VAREST*VARSQR
109       FINEST = FINEST + ( FINVAL - FINEST )/(1 + VARPRD)
110       IF ( VARSQR .GT. 0 ) VAREST = (1 + VARPRD)/VARSQR
111       ABSEST = 3*SQRT( VARSQR/( 1 + VARPRD ) )
112       END
113       SUBROUTINE SPHMVN(N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
114      &     ABSEPS, RELEPS, ERROR, VALUE, INFORM)
115 *
116 *     A subroutine for computing multivariate normal probabilities.
117 *     This subroutine uses a Mont-Carlo algorithm given in the paper
118 *       "Three Digit Accurate Multiple Normal Probabilities",
119 *          pp. 369-380, Numer. Math. 35(1980), by I. Deak
120 *
121 *
122 *  Parameters
123 *
124 *     N      INTEGER, the number of variables.
125 *     LOWER  REAL, array of lower integration limits.
126 *     UPPER  REAL, array of upper integration limits.
127 *     INFIN  INTEGER, array of integration limits flags:
128 *            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
129 *            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
130 *            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
131 *            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
132 *     CORREL REAL, array of correlation coefficients; the correlation
133 *            coefficient in row I column J of the correlation matrix
134 *            should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
135 *     MAXPTS INTEGER, maximum number of function values allowed. This
136 *            parameter can be used to limit the time. A sensible
137 *            strategy is to start with MAXPTS = 1000*N, and then
138 *            increase MAXPTS if ERROR is too large.
139 *     ABSEPS REAL absolute error tolerance.
140 *     RELEPS REAL relative error tolerance.
141 *     ERROR  REAL, estimated absolute error, with 99% confidence level.
142 *     VALUE  REAL, estimated value for the integral
143 *     INFORM INTEGER, termination status parameter:
144 *            if INFORM = 0, normal completion with ERROR < EPS;
145 *            if INFORM = 1, completion with ERROR > EPS and MAXPTS
146 *                           function vaules used; increase MAXPTS to
147 *                           decrease ERROR;
148 *            if INFORM = 2, N > 100.
149 *
150       EXTERNAL SPNRML
151       INTEGER N, INFIS, INFIN(*), MAXPTS, MPT, INFORM, NS, IVLS
152       DOUBLE PRECISION CORREL(*), LOWER(*), UPPER(*),
153      &     ABSEPS, RELEPS, ERROR, VALUE, D, E, EPS, SPNRNT
154       IF ( N .GT. 100 ) THEN
155          INFORM = 2
156          VALUE = 0
157          ERROR = 1
158          RETURN
159       ENDIF
160       INFORM = SPNRNT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E, NS)
161       IF ( N-INFIS .EQ. 0 ) THEN
162          VALUE = 1
163          ERROR = 0
164       ELSE IF ( N-INFIS .EQ. 1 ) THEN
165          VALUE = E - D
166          ERROR = 2E-16
167       ELSE
168 *
169 *        Call then Monte-Carlo integration subroutine
170 *
171          MPT = 25 + NS/N**3
172          CALL SCRUDE( N-INFIS, MPT, ERROR, VALUE, 0 )
173          IVLS = MPT*NS
174  10      EPS = MAX( ABSEPS, RELEPS*ABS(VALUE) )
175          IF ( ERROR .GT. EPS .AND. IVLS .LT. MAXPTS ) THEN
176             MPT = MAX( MIN( INT(MPT*(ERROR/(EPS))**2),
177      &                      ( MAXPTS - IVLS )/NS ), 10 )
178             CALL SCRUDE( N-INFIS, MPT, ERROR, VALUE, 1 )
179             IVLS = IVLS + MPT*NS
180             GO TO 10
181          ENDIF
182          IF ( ERROR. GT. EPS .AND. IVLS .GE. MAXPTS ) INFORM = 1
183       ENDIF
184       END
185       DOUBLE PRECISION FUNCTION SPNRML(N)
186 *
187 *     Integrand subroutine
188 *
189       DOUBLE PRECISION LOWER(*), UPPER(*), CORREL(*), D, E, ZERO
190       INTEGER N, INFIN(*), INFIS
191       INTEGER NL, IJ, I, J, K, NS, NSO, ND
192       PARAMETER ( NL = 100, ND = 3, ZERO = 0 )
193       DOUBLE PRECISION A(NL), B(NL), U(NL,NL), Y(NL)
194       INTEGER INFI(NL), IS(NL), IC(NL)
195       DOUBLE PRECISION RS, TMP, BT, RNOR, SPHLIM, SPNRNT
196       SAVE A, B, INFI, U
197 *
198 *    First generate U = COV*(random orthogonal matrix)
199 *
200       DO K = N-1, 1, -1
201          TMP = 0
202          DO J = K, N
203             Y(J) = RNOR()
204             TMP = TMP + Y(J)**2
205          END DO
206          TMP = -SQRT(TMP)
207          BT = 1/( TMP*( Y(K) + TMP ) )
208          Y(K) = Y(K) + TMP
209          DO I = 1, N
210             TMP = 0
211             DO J = K, N
212                TMP = TMP + U(I,J)*Y(J)
213             END DO
214             TMP = BT*TMP
215             DO J = K, N
216                U(I,J) = U(I,J) - TMP*Y(J)
217             END DO
218          END DO
219       END DO
220 *
221 *     Compute integrand average
222 *
223       RS = SQRT( DBLE(ND) )
224       DO I = 1,ND
225          IC(I) = I
226       END DO
227       IC(ND+1) = N+1
228       SPNRML = 0
229       NS = 0
230  10   DO I = 1,ND
231          IS(I) = -1
232       END DO
233  20   DO I = 1, N
234          TMP = 0
235          DO J = 1,ND
236             TMP = TMP + IS(J)*U(I,IC(J))
237          END DO
238          Y(I) = TMP/RS
239       END DO
240       NS = NS + 1
241       SPNRML = SPNRML + ( SPHLIM( N, A, B, INFI, Y ) - SPNRML )/NS
242       DO I = 1, ND
243          IS(I) = IS(I) + 2
244          IF ( IS(I) .LT. 2 ) GO TO 20
245          IS(I) = -1
246       END DO
247       DO I = 1, ND
248          IC(I) = IC(I) + 1
249          IF ( IC(I) .LT. IC(I+1)  ) GO TO 10
250          IC(I) = I
251       END DO
252       SPNRML = SPNRML/2
253       RETURN
254       ENTRY SPNRNT( N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E, NSO )
255       SPNRNT = 0
256 *
257 *     Initialisation
258 *
259       IJ = 0
260       INFIS = 0
261       DO I = 1, N
262          INFI(I) = INFIN(I)
263          IF ( INFI(I) .LT. 0 ) THEN
264             INFIS = INFIS + 1
265          ELSE
266             A(I) = 0
267             B(I) = 0
268             IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I)
269             IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I)
270          ENDIF
271          DO J = 1, I-1
272             IJ = IJ + 1
273             U(I,J) = CORREL(IJ)
274             U(J,I) = 0
275          END DO
276          U(I,I) = 1
277       END DO
278       NSO = 1
279       DO I = 1,ND
280          NSO = 2*NSO*( N - INFIS - I + 1 )/I
281       END DO
282 *
283 *     First move any doubly infinite limits to innermost positions
284 *
285       IF ( INFIS .LT. N ) THEN
286          outer: DO I = N, N-INFIS+1, -1
287             IF ( INFI(I) .GE. 0 ) THEN
288                DO J = 1,I-1
289                   IF ( INFI(J) .LT. 0 ) THEN
290                      DO K = 1, J-1
291                         TMP = U(J,K)
292                         U(J,K) = U(I,K)
293                         U(I,K) = TMP
294                      END DO
295                      DO K = J+1, I-1
296                         TMP = U(I,K)
297                         U(I,K) = U(K,J)
298                         U(K,J) = TMP
299                      END DO
300                      DO K = I+1, N
301                         TMP = U(K,J)
302                         U(K,J) = U(K,I)
303                         U(K,I) = TMP
304                      END DO
305                      TMP = A(J)
306                      A(J) = A(I)
307                      A(I) = TMP
308                      TMP = B(J)
309                      B(J) = B(I)
310                      B(I) = TMP
311                      TMP = INFI(J)
312                      INFI(J) = INFI(I)
313                      INFI(I) = TMP
314                      CYCLE outer
315                   ENDIF
316                END DO
317             ENDIF
318          END DO outer
319       ENDIF
320 *
321 *     Determine Cholesky decomposition
322 *
323       DO J = 1, N-INFIS
324          DO I = J, N-INFIS
325             TMP = U(I,J)
326             DO K = 1, J-1
327                TMP = TMP - U(I,K)*U(J,K)
328             END DO
329             IF ( I .EQ. J ) THEN
330                U(J,J) = SQRT( MAX( TMP, ZERO ) )
331             ELSE IF ( U(I,I) .GT. 0 ) THEN
332                U(I,J) = TMP/U(J,J)
333             ELSE
334                U(I,J) = 0
335             END IF
336          END DO
337       END DO
338       DO I = 1, N-INFIS
339          IF ( U(I,I) .GT. 0 ) THEN
340             IF ( INFI(I) .NE. 0 ) A(I) = A(I)/U(I,I)
341             IF ( INFI(I) .NE. 1 ) B(I) = B(I)/U(I,I)
342             DO J = 1,I
343                U(I,J) = U(I,J)/U(I,I)
344             END DO
345          ENDIF
346       END DO
347       CALL LIMITS( A(1), B(1), INFI(1), D, E )
348       END
349       DOUBLE PRECISION FUNCTION SPHLIM( N, A, B, INFI, Y )
350       DOUBLE PRECISION A(*), B(*), Y(*), CMN, CMX, SPHINC
351       INTEGER INFI(*), I, N
352       CMN = -10*N
353       CMX =  10*N
354       DO I = 1,N
355          IF ( Y(I) .GT. 0 ) THEN
356             IF ( INFI(I) .NE. 1 ) CMX = MIN( CMX, B(I)/Y(I) )
357             IF ( INFI(I) .NE. 0 ) CMN = MAX( CMN, A(I)/Y(I) )
358          ELSE
359             IF ( INFI(I) .NE. 1 ) CMN = MAX( CMN, B(I)/Y(I) )
360             IF ( INFI(I) .NE. 0 ) CMX = MIN( CMX, A(I)/Y(I) )
361          ENDIF
362       END DO
363       IF ( CMN .LT. CMX ) THEN
364          IF ( CMN .GE. 0 .AND. CMX .GE. 0 ) THEN
365             SPHLIM = SPHINC( N,  CMX ) - SPHINC( N,  CMN )
366          ELSEIF ( CMN .LT. 0 .AND. CMX .GE. 0 ) THEN
367             SPHLIM = SPHINC( N, -CMN ) + SPHINC( N,  CMX )
368          ELSE
369             SPHLIM = SPHINC( N, -CMN ) - SPHINC( N, -CMX )
370          ENDIF
371       ELSE
372          SPHLIM = 0
373       ENDIF
374       END
375       SUBROUTINE SCRUDE( NDIM, MAXPTS, ABSEST, FINEST, IR )
376 *
377 *     Crude Monte-Carlo Algorithm for Deak method with
378 *      weighted results on restart
379 *
380       INTEGER NDIM, MAXPTS, M, K, IR, NPTS
381       DOUBLE PRECISION FINEST, ABSEST, SPNRML, UNI,
382      &     VARSQR, VAREST, VARPRD, FINDIF, FINVAL
383       SAVE VAREST
384       IF ( IR .LE. 0 ) THEN
385          VAREST = 0
386          FINEST = 0
387       ENDIF
388       FINVAL = 0
389       VARSQR = 0
390       DO M = 1,MAXPTS
391          FINDIF = ( SPNRML(NDIM) - FINVAL )/M
392          FINVAL = FINVAL + FINDIF
393          VARSQR = ( M - 2 )*VARSQR/M + FINDIF**2
394       END DO
395       VARPRD = VAREST*VARSQR
396       FINEST = FINEST + ( FINVAL - FINEST )/(1 + VARPRD)
397       IF ( VARSQR .GT. 0 ) VAREST = (1 + VARPRD)/VARSQR
398       ABSEST = 3*SQRT( VARSQR/( 1 + VARPRD ) )
399       END
400       DOUBLE PRECISION FUNCTION SPHINC( N, R )
401 *
402 *                   R
403 *     SPHINC =  K  I  exp(-t*t/2) t**(N-1) dt, for N > 1.
404 *                N  0
405 *
406       INTEGER I, N
407       DOUBLE PRECISION R, RR, RP, PF, ET, PHI
408       PARAMETER ( RP = 2.5066 28274 63100 04D0 )
409       IF ( R .GT. 0 ) THEN
410          RR = R*R
411          PF = 1
412          DO I = N-2, 2, -2
413             PF = 1 + RR*PF/I
414          END DO
415          IF ( MOD( N, 2 ) .EQ. 0 ) THEN
416             ET = LOG(PF) - RR/2
417             IF ( ET .GT. -40 ) THEN
418                SPHINC = 1 - EXP( ET )
419             ELSE
420                SPHINC = 1
421             END IF
422          ELSE
423             SPHINC = 1  - 2*PHI(-R)
424             ET = LOG(R*PF) - RR/2
425             IF ( ET .GT. -40 ) SPHINC = SPHINC - 2*EXP( ET )/RP
426          ENDIF
427       ELSE
428          SPHINC = 0
429       ENDIF
430       END
431       DOUBLE PRECISION FUNCTION RNOR()
432 *
433 *     RNOR generates normal random numbers with zero mean and unit
434 *     standard deviation, often denoted N(0,1),adapted from G. Marsaglia
435 *     and W. W. Tsang: "A Fast, Easily Implemented Method for Sampling
436 *     from Decreasing or Symmetric Unimodal Density Functions"
437 *      SIAM J. Sci. Stat. Comput. 5(1984), pp. 349-359.
438 *
439       INTEGER J, N, TN
440       DOUBLE PRECISION TWOPIS, AA, B, C, XDN
441       PARAMETER ( N = 64, TN = 2*N, TWOPIS = TN/2.506628274631000D0 )
442       PARAMETER ( XDN = 0.3601015713011893D0, B = 0.4878991777603940D0 )
443       PARAMETER (  AA =  12.37586029917064D0, C =  12.67705807886560D0 )
444       DOUBLE PRECISION XT, XX, Y, UNI
445       DOUBLE PRECISION X(0:N)
446       SAVE X
447       DATA ( X(J), J = 0, 31 ) /
448      &  0.3409450287039653D+00,  0.4573145918669259D+00,
449      &  0.5397792816116612D+00,  0.6062426796530441D+00,
450      &  0.6631690627645207D+00,  0.7136974590560222D+00,
451      &  0.7596124749339174D+00,  0.8020356003555283D+00,
452      &  0.8417226679789527D+00,  0.8792102232083114D+00,
453      &  0.9148948043867484D+00,  0.9490791137530882D+00,
454      &  0.9820004812398864D+00,  0.1013849238029940D+01,
455      &  0.1044781036740172D+01,  0.1074925382028552D+01,
456      &  0.1104391702268125D+01,  0.1133273776243940D+01,
457      &  0.1161653030133931D+01,  0.1189601040838737D+01,
458      &  0.1217181470700870D+01,  0.1244451587898246D+01,
459      &  0.1271463480572119D+01,  0.1298265041883197D+01,
460      &  0.1324900782180860D+01,  0.1351412509933371D+01,
461      &  0.1377839912870011D+01,  0.1404221063559975D+01,
462      &  0.1430592868502691D+01,  0.1456991476137671D+01,
463      &  0.1483452656603219D+01,  0.1510012164318519D+01 /
464       DATA ( X(J), J = 32, 64 ) /
465      &  0.1536706093359520D+01,  0.1563571235037691D+01,
466      &  0.1590645447014253D+01,  0.1617968043674446D+01,
467      &  0.1645580218369081D+01,  0.1673525509567038D+01,
468      &  0.1701850325062740D+01,  0.1730604541317782D+01,
469      &  0.1759842199038300D+01,  0.1789622321566574D+01,
470      &  0.1820009890130691D+01,  0.1851077020230275D+01,
471      &  0.1882904397592872D+01,  0.1915583051943031D+01,
472      &  0.1949216574916360D+01,  0.1983923928905685D+01,
473      &  0.2019843052906235D+01,  0.2057135559990095D+01,
474      &  0.2095992956249391D+01,  0.2136645022544389D+01,
475      &  0.2179371340398135D+01,  0.2224517507216017D+01,
476      &  0.2272518554850147D+01,  0.2323933820094302D+01,
477      &  0.2379500774082828D+01,  0.2440221797979943D+01,
478      &  0.2507511701865317D+01,  0.2583465835225429D+01,
479      &  0.2671391590320836D+01,4*0.2776994269662875D+01 /
480       Y = UNI()
481       J = MOD( INT( TN*UNI() ), N )
482       XT = X(J+1)
483       RNOR = ( Y + Y - 1 )*XT
484       IF ( ABS(RNOR) .GT. X(J) ) THEN
485          XX = B*( XT - ABS(RNOR) )/( XT - X(J) )
486          Y = UNI()
487          IF ( Y .GT. C - AA*EXP( -XX**2/2 ) ) THEN
488             RNOR = SIGN( XX, RNOR )
489          ELSE
490             IF ( EXP(-XT**2/2)+Y/(TWOPIS*XT).GT.EXP(-RNOR**2/2) ) THEN
491  10            XX = XDN*LOG( UNI() )
492                IF ( -2*LOG( UNI() ) .LE. XX**2 ) GO TO 10
493                RNOR = SIGN( X(N) - XX, RNOR )
494             END IF
495          END IF
496       END IF
497       END
498 *
499       DOUBLE PRECISION FUNCTION UNI()
500 *
501 *     Uniform (0, 1) random number generator
502 *
503 *     Reference:
504 *     L'Ecuyer, Pierre (1996),
505 *     "Combined Multiple Recursive Random Number Generators"
506 *     Operations Research 44, pp. 816-822.
507 *
508 *
509       INTEGER A12, A13, A21, A23, P12, P13, P21, P23
510       INTEGER Q12, Q13, Q21, Q23, R12, R13, R21, R23
511       INTEGER X10, X11, X12, X20, X21, X22, Z, M1, M2, H
512       DOUBLE PRECISION INVMP1
513       PARAMETER (  M1 = 2147483647,  M2 = 2145483479 )
514       PARAMETER ( A12 =   63308,    Q12 = 33921, R12 = 12979 )
515       PARAMETER ( A13 = -183326,    Q13 = 11714, R13 =  2883 )
516       PARAMETER ( A21 =   86098,    Q21 = 24919, R21 =  7417 )
517       PARAMETER ( A23 = -539608,    Q23 =  3976, R23 =  2071 )
518       PARAMETER ( INVMP1 = 4.656612873077392578125D-10 )
519 *                 INVMP1 = 1/(M1+1)
520       SAVE X10, X11, X12, X20, X21, X22
521       DATA       X10,      X11,      X12,      X20,      X21,      X22
522      &    / 11111111, 22222223, 33333335, 44444447, 55555559, 66666661 /
523 *
524 *     Component 1
525 *
526       H = X10/Q13
527       P13 = -A13*( X10 - H*Q13 ) - H*R13
528       H = X11/Q12
529       P12 =  A12*( X11 - H*Q12 ) - H*R12
530       IF ( P13 .LT. 0 ) P13 = P13 + M1
531       IF ( P12 .LT. 0 ) P12 = P12 + M1
532       X10 = X11
533       X11 = X12
534       X12 = P12 - P13
535       IF ( X12 .LT. 0 ) X12 = X12 + M1
536 *
537 *     Component 2
538 *
539       H = X20/Q23
540       P23 = -A23*( X20 - H*Q23 ) - H*R23
541       H = X22/Q21
542       P21 =  A21*( X22 - H*Q21 ) - H*R21
543       IF ( P23 .LT. 0 ) P23 = P23 + M2
544       IF ( P21 .LT. 0 ) P21 = P21 + M2
545       X20 = X21
546       X21 = X22
547       X22 = P21 - P23
548       IF ( X22 .LT. 0 ) X22 = X22 + M2
549 *
550 *     Combination
551 *
552       Z = X12 - X22
553       IF ( Z .LE. 0 ) Z = Z + M1
554       UNI = Z*INVMP1
555       END
556       SUBROUTINE SADMVN( N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
557      &                   ABSEPS, RELEPS, ERROR, VALUE, INFORM )
558 *
559 *     A subroutine for computing multivariate normal probabilities.
560 *     This subroutine uses an algorithm given in the paper
561 *     "Numerical Computation of Multivariate Normal Probabilities", in
562 *     J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
563 *          Alan Genz
564 *          Department of Mathematics
565 *          Washington State University
566 *          Pullman, WA 99164-3113
567 *          Email : alangenz@wsu.edu
568 *
569 *  Parameters
570 *
571 *     N      INTEGER, the number of variables.
572 *     LOWER  REAL, array of lower integration limits.
573 *     UPPER  REAL, array of upper integration limits.
574 *     INFIN  INTEGER, array of integration limits flags:
575 *            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
576 *            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
577 *            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
578 *            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
579 *     CORREL REAL, array of correlation coefficients; the correlation
580 *            coefficient in row I column J of the correlation matrix
581 *            should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
582 *     MAXPTS INTEGER, maximum number of function values allowed. This
583 *            parameter can be used to limit the time taken. A
584 *            sensible strategy is to start with MAXPTS = 1000*N, and then
585 *            increase MAXPTS if ERROR is too large.
586 *     ABSEPS REAL absolute error tolerance.
587 *     RELEPS REAL relative error tolerance.
588 *     ERROR  REAL estimated absolute error, with 99% confidence level.
589 *     VALUE  REAL estimated value for the integral
590 *     INFORM INTEGER, termination status parameter:
591 *            if INFORM = 0, normal completion with ERROR < EPS;
592 *            if INFORM = 1, completion with ERROR > EPS and MAXPTS
593 *                           function vaules used; increase MAXPTS to
594 *                           decrease ERROR;
595 *            if INFORM = 2, N > 20 or N < 1.
596 *
597       EXTERNAL MVNFNC
598       INTEGER N, NL, M, INFIN(*), LENWRK, MAXPTS, INFORM, INFIS,
599      &     RULCLS, TOTCLS, NEWCLS, MAXCLS
600       DOUBLE PRECISION
601      &     CORREL(*), LOWER(*), UPPER(*), ABSEPS, RELEPS, ERROR, VALUE,
602      &     OLDVAL, D, E, MVNNIT, MVNFNC
603       PARAMETER ( NL = 20 )
604       PARAMETER ( LENWRK = 20*NL**2 )
605       DOUBLE PRECISION WORK(LENWRK)
606       IF ( N .GT. 20 .OR. N .LT. 1 ) THEN
607          INFORM = 2
608          VALUE = 0
609          ERROR = 1
610          RETURN
611       ENDIF
612       INFORM = MVNNIT( N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E )
613       M = N - INFIS
614       IF ( M .EQ. 0 ) THEN
615          VALUE = 1
616          ERROR = 0
617       ELSE IF ( M .EQ. 1 ) THEN
618          VALUE = E - D
619          ERROR = 2E-16
620       ELSE
621 *
622 *        Call the subregion adaptive integration subroutine
623 *
624          M = M - 1
625          RULCLS = 1
626          CALL ADAPT( M, RULCLS, 0, MVNFNC, ABSEPS, RELEPS,
627      &               LENWRK, WORK, ERROR, VALUE, INFORM )
628          MAXCLS = MIN( 10*RULCLS, MAXPTS )
629          TOTCLS = 0
630          CALL ADAPT(M, TOTCLS, MAXCLS, MVNFNC, ABSEPS, RELEPS,
631      &        LENWRK, WORK, ERROR, VALUE, INFORM)
632          IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
633  10         OLDVAL = VALUE
634             MAXCLS = MAX( 2*RULCLS, MIN( 3*MAXCLS/2, MAXPTS - TOTCLS ) )
635             NEWCLS = -1
636             CALL ADAPT(M, NEWCLS, MAXCLS, MVNFNC, ABSEPS, RELEPS,
637      &           LENWRK, WORK, ERROR, VALUE, INFORM)
638             TOTCLS = TOTCLS + NEWCLS
639             ERROR = ABS(VALUE-OLDVAL) + SQRT(RULCLS*ERROR**2/TOTCLS)
640             IF ( ERROR .GT. MAX( ABSEPS, RELEPS*ABS(VALUE) ) ) THEN
641                IF ( MAXPTS - TOTCLS .GT. 2*RULCLS ) GO TO 10
642             ELSE
643                INFORM = 0
644             END IF
645          ENDIF
646       ENDIF
647       END
648 *
649       SUBROUTINE KROMVN( N, LOWER, UPPER, INFIN, CORREL, MAXPTS,
650      &                   ABSEPS, RELEPS, ERROR, VALUE, INFORM )
651 *
652 *     A subroutine for computing multivariate normal probabilities.
653 *     This subroutine uses an algorithm given in the paper
654 *     "Numerical Computation of Multivariate Normal Probabilities", in
655 *     J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
656 *          Alan Genz
657 *          Department of Mathematics
658 *          Washington State University
659 *          Pullman, WA 99164-3113
660 *          Email : AlanGenz@wsu.edu
661 *
662 *  Parameters
663 *
664 *     N      INTEGER, the number of variables.
665 *     LOWER  REAL, array of lower integration limits.
666 *     UPPER  REAL, array of upper integration limits.
667 *     INFIN  INTEGER, array of integration limits flags:
668 *            if INFIN(I) < 0, Ith limits are (-infinity, infinity);
669 *            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
670 *            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
671 *            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
672 *     CORREL REAL, array of correlation coefficients; the correlation
673 *            coefficient in row I column J of the correlation matrix
674 *            should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I.
675 *     MAXPTS INTEGER, maximum number of function values allowed. This
676 *            parameter can be used to limit the time. A sensible
677 *            strategy is to start with MAXPTS = 1000*N, and then
678 *            increase MAXPTS if ERROR is too large.
679 *     ABSEPS REAL absolute error tolerance.
680 *     RELEPS REAL relative error tolerance.
681 *     ERROR  REAL estimated absolute error, with 99% confidence level.
682 *     VALUE  REAL estimated value for the integral
683 *     INFORM INTEGER, termination status parameter:
684 *            if INFORM = 0, normal completion with ERROR < EPS;
685 *            if INFORM = 1, completion with ERROR > EPS and MAXPTS
686 *                           function vaules used; increase MAXPTS to
687 *                           decrease ERROR;
688 *            if INFORM = 2, N > 100 or N < 1.
689 *
690       EXTERNAL MVNFNC
691       INTEGER N, INFIN(*), MAXPTS, INFORM, INFIS, IVLS
692       DOUBLE PRECISION CORREL(*), LOWER(*), UPPER(*), RELEPS, ABSEPS,
693      &       ERROR, VALUE, E, D, MVNNIT, MVNFNC
694       IF ( N .GT. 100 .OR. N .LT. 1 ) THEN
695          INFORM = 2
696          VALUE = 0
697          ERROR = 1
698       ELSE
699          INFORM = MVNNIT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E)
700          IF ( N-INFIS .EQ. 0 ) THEN
701             VALUE = 1
702             ERROR = 0
703          ELSE IF ( N-INFIS .EQ. 1 ) THEN
704             VALUE = E - D
705             ERROR = 2E-16
706          ELSE
707 *
708 *        Call the lattice rule integration subroutine
709 *
710             IVLS = 0
711             CALL KROBOV( N-INFIS-1, IVLS, MAXPTS, MVNFNC,
712      &                   ABSEPS, RELEPS, ERROR, VALUE, INFORM )
713          ENDIF
714       ENDIF
715       END
716       SUBROUTINE KROBOV( NDIM, MINVLS, MAXVLS, FUNCTN, ABSEPS, RELEPS,
717      &                   ABSERR, FINEST, INFORM )
718 *
719 *  Automatic Multidimensional Integration Subroutine
720 *
721 *         AUTHOR: Alan Genz
722 *                 Department of Mathematics
723 *                 Washington State University
724 *                 Pulman, WA 99164-3113
725 *                 Email: AlanGenz@wsu.edu
726 *
727 *         Last Change: 4/15/98
728 *
729 *  KROBOV computes an approximation to the integral
730 *
731 *      1  1     1
732 *     I  I ... I       F(X)  dx(NDIM)...dx(2)dx(1)
733 *      0  0     0
734 *
735 *
736 *  KROBOV uses randomized Korobov rules. The primary references are
737 *  "Randomization of Number Theoretic Methods for Multiple Integration"
738 *   R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13, pp. 904-14,
739 *  and
740 *   "Optimal Parameters for Multidimensional Integration",
741 *    P. Keast, SIAM J Numer Anal, 10, pp.831-838.
742 *
743 ***************  Parameters ********************************************
744 ****** Input parameters
745 *  NDIM    Number of variables, must exceed 1, but not exceed 40
746 *  MINVLS  Integer minimum number of function evaluations allowed.
747 *          MINVLS must not exceed MAXVLS.  If MINVLS < 0 then the
748 *          routine assumes a previous call has been made with
749 *          the same integrand and continues that calculation.
750 *  MAXVLS  Integer maximum number of function evaluations allowed.
751 *  FUNCTN  EXTERNALly declared user defined function to be integrated.
752 *          It must have parameters (NDIM,Z), where Z is a real array
753 *          of dimension NDIM.
754 *  ABSEPS  Required absolute accuracy.
755 *  RELEPS  Required relative accuracy.
756 ****** Output parameters
757 *  MINVLS  Actual number of function evaluations used.
758 *  ABSERR  Estimated absolute accuracy of FINEST.
759 *  FINEST  Estimated value of integral.
760 *  INFORM  INFORM = 0 for normal exit, when
761 *                     ABSERR <= MAX(ABSEPS, RELEPS*ABS(FINEST))
762 *                  and
763 *                     INTVLS <= MAXCLS.
764 *          INFORM = 1 If MAXVLS was too small to obtain the required
765 *          accuracy. In this case a value FINEST is returned with
766 *          estimated absolute accuracy ABSERR.
767 ************************************************************************
768       EXTERNAL FUNCTN
769       INTEGER NDIM, MINVLS, MAXVLS, INFORM, NP, PLIM, NLIM,
770      &        SAMPLS, I, INTVLS, MINSMP
771       PARAMETER ( PLIM = 20, NLIM = 100, MINSMP = 6 )
772       INTEGER C(PLIM,NLIM), P(PLIM)
773       DOUBLE PRECISION FUNCTN, ABSEPS, RELEPS, FINEST, ABSERR, DIFINT,
774      &       FINVAL, VARSQR, VAREST, VARPRD, VALUE
775       DOUBLE PRECISION ALPHA(NLIM), X(NLIM), VK(NLIM), ONE
776       PARAMETER ( ONE = 1 )
777       SAVE P, C, SAMPLS, NP, VAREST
778       INFORM = 1
779       INTVLS = 0
780       IF ( MINVLS .GE. 0 ) THEN
781          FINEST = 0
782          VAREST = 0
783          SAMPLS = MINSMP
784          DO I = 1, PLIM
785             NP = I
786             IF ( MINVLS .LT. 2*SAMPLS*P(I) ) GO TO 10
787          END DO
788          SAMPLS = MAX( MINSMP, MINVLS/( 2*P(NP) ) )
789       ENDIF
790  10   VK(1) = ONE/P(NP)
791       DO I = 2, NDIM
792          VK(I) = MOD( C(NP,NDIM-1)*VK(I-1), ONE )
793       END DO
794       FINVAL = 0
795       VARSQR = 0
796       DO I = 1, SAMPLS
797          CALL KROSUM( NDIM, VALUE, P(NP), VK, FUNCTN, ALPHA, X )
798          DIFINT = ( VALUE - FINVAL )/I
799          FINVAL = FINVAL + DIFINT
800          VARSQR = ( I - 2 )*VARSQR/I + DIFINT**2
801       END DO
802       INTVLS = INTVLS + 2*SAMPLS*P(NP)
803       VARPRD = VAREST*VARSQR
804       FINEST = FINEST + ( FINVAL - FINEST )/( 1 + VARPRD )
805       IF ( VARSQR .GT. 0 ) VAREST = ( 1 + VARPRD )/VARSQR
806       ABSERR = 3*SQRT( VARSQR/( 1 + VARPRD ) )
807       IF ( ABSERR .GT. MAX( ABSEPS, ABS(FINEST)*RELEPS ) ) THEN
808          IF ( NP .LT. PLIM ) THEN
809             NP = NP + 1
810          ELSE
811             SAMPLS = MIN( 3*SAMPLS/2, ( MAXVLS - INTVLS )/( 2*P(NP) ) )
812             SAMPLS = MAX( MINSMP, SAMPLS )
813          ENDIF
814          IF ( INTVLS + 2*SAMPLS*P(NP) .LE. MAXVLS ) GO TO 10
815       ELSE
816          INFORM = 0
817       ENDIF
818       MINVLS = INTVLS
819       DATA P( 1), ( C( 1,I), I = 1, 99 ) /    113,
820      &     42,    54,    55,    32,    13,    26,    26,    13,    26,
821      &     14,    13,    26,    35,     2,     2,     2,     2,    56,
822      &     28,     7,     7,    28,     4,    49,     4,    40,    48,
823      &      5,    35,    27,    16,    16,     2,     2,     7,    28,
824      &      4,    49,     4,    56,     8,     2,     2,    56,     7,
825      &     16,    28,     7,     7,    28,     4,    49,     4,    37,
826      &     55,    21,    33,    40,    16,    16,    28,     7,    16,
827      &     28,     4,    49,     4,    56,    35,     2,     2,     2,
828      &     16,    16,    28,     4,    16,    28,     4,    49,     4,
829      &     40,    40,     5,    42,    27,    16,    16,    28,     4,
830      &     16,    28,     4,    49,     4,     8,     8,     2,     2/
831       DATA P( 2), ( C( 2,I), I = 1, 99 ) /    173,
832      &     64,    34,    57,     9,    72,    86,    16,    75,    75,
833      &     70,    42,     2,    86,    62,    62,    30,    30,     5,
834      &     42,    70,    70,    70,    53,    70,    70,    53,    42,
835      &     62,    53,    53,    53,    69,    75,     5,    53,    86,
836      &      2,     5,    30,    75,    59,     2,    69,     5,     5,
837      &     63,    62,     5,    69,    30,    44,    30,    86,    86,
838      &      2,    69,     5,     5,     2,     2,    61,    69,    17,
839      &      2,     2,     2,    53,    69,     2,     2,    86,    69,
840      &     13,     2,     2,    37,    43,    65,     2,     2,    30,
841      &     86,    45,    16,    32,    18,    86,    86,    86,     9,
842      &     63,    63,    11,    76,    76,    76,    63,    60,    70/
843       DATA P( 3), ( C( 3,I), I = 1, 99 ) /    263,
844      &    111,    67,    98,    36,    48,   110,     2,   131,     2,
845      &      2,   124,   124,    48,     2,     2,   124,   124,    70,
846      &     70,    48,   126,    48,   126,    56,    65,    48,    48,
847      &     70,     2,    92,   124,    92,   126,   131,   124,    70,
848      &     70,    70,    20,   105,    70,     2,     2,    27,   108,
849      &     27,    39,     2,   131,   131,    92,    92,    48,     2,
850      &    126,    20,   126,     2,     2,   131,    38,   117,     2,
851      &    131,    68,    58,    38,    90,    38,   108,    38,     2,
852      &    131,   131,   131,    68,    14,    94,   131,   131,   131,
853      &    108,    18,   131,    56,    85,   117,   117,     9,   131,
854      &    131,    55,    92,    92,    92,   131,   131,    48,    48/
855       DATA P( 4), ( C( 4,I), I = 1, 99 ) /    397,
856      &    151,   168,    46,   197,    69,    64,     2,   198,   191,
857      &    134,   134,   167,   124,    16,   124,   124,   124,   124,
858      &    141,   134,   128,     2,     2,    32,    32,    32,    31,
859      &     31,    64,    64,    99,     4,     4,   167,   124,   124,
860      &    124,   124,   124,   124,   107,    85,    79,    85,   111,
861      &     85,   128,    31,    31,    31,    31,    64,   167,     4,
862      &    107,   167,   124,   124,   124,   124,   124,   124,   107,
863      &    183,     2,     2,     2,    62,    32,    31,    31,    31,
864      &     31,    31,   167,     4,   107,   167,   124,   124,   124,
865      &    124,   124,   124,   107,   142,   184,   184,    65,    65,
866      &    183,    31,    31,    31,    31,    31,   167,     4,   107/
867       DATA P( 5), ( C( 5,I), I = 1, 99 ) /    593,
868      &    229,    40,   268,    42,   153,   294,    71,     2,   130,
869      &    199,   199,   199,   149,   199,   149,   153,   130,   149,
870      &    149,    15,   119,   294,    31,    82,   260,   122,   209,
871      &    209,   122,   296,   130,   130,   260,   260,    30,   206,
872      &     94,   209,    94,   122,   209,   209,   122,   122,   209,
873      &    130,     2,   130,   130,    38,    38,    79,    82,    94,
874      &     82,   122,   122,   209,   209,   122,   122,   168,   220,
875      &     62,    60,   168,   282,   282,    82,   209,   122,    94,
876      &    209,   122,   122,   122,   122,   258,   148,   286,   256,
877      &    256,    62,    62,    82,   122,    82,    82,   122,   122,
878      &    122,   209,   122,    15,    79,    79,    79,    79,   168/
879       DATA P( 6), ( C( 6,I), I = 1, 99 ) /    907,
880      &    264,   402,   406,   147,   452,   153,   224,     2,     2,
881      &    224,   224,   449,   101,   182,   449,   101,   451,   181,
882      &    181,   101,   101,   377,    85,   453,   453,   453,    85,
883      &    197,   451,     2,     2,   101,   449,   449,   449,   173,
884      &    173,     2,   453,   453,     2,   426,    66,   367,   426,
885      &    101,   453,     2,    32,    32,    32,   101,     2,     2,
886      &    453,   223,   147,   449,   290,     2,   453,     2,    83,
887      &    223,   101,   453,     2,    83,    83,   147,     2,   453,
888      &    147,   147,   147,   147,   147,   147,   147,   453,   153,
889      &    153,   147,     2,   224,   290,   320,   453,   147,   431,
890      &    383,   290,   290,     2,   162,   162,   147,     2,   162/
891       DATA P( 7), ( C( 7,I), I = 1, 99 ) /   1361,
892      &    505,   220,   195,   410,   199,   248,   460,   471,     2,
893      &    331,   662,   547,   209,   547,   547,   209,     2,   680,
894      &    680,   629,   370,   574,    63,    63,   259,   268,   259,
895      &    547,   209,   209,   209,   547,   547,   209,   209,   547,
896      &    547,   108,    63,    63,   108,    63,    63,   108,   259,
897      &    268,   268,   547,   209,   209,   209,   209,   547,   209,
898      &    209,   209,   547,   108,    63,    63,    63,   405,   285,
899      &    234,   259,   259,   259,   259,   209,   209,   209,   209,
900      &    209,   209,   209,   209,   547,   289,   289,   234,   285,
901      &    316,     2,   410,   259,   259,   259,   268,   209,   209,
902      &    209,   209,   547,   547,   209,   209,   209,   285,   316/
903       DATA P( 8), ( C( 8,I), I = 1, 99 ) /   2053,
904      &    468,   635,   849,   687,   948,    37,  1014,   513,     2,
905      &      2,     2,     2,     2,  1026,     2,     2,  1026,   201,
906      &    201,     2,  1026,   413,  1026,  1026,     2,     2,   703,
907      &    703,     2,     2,   393,   393,   678,   413,  1026,     2,
908      &      2,  1026,  1026,     2,   405,   953,     2,  1026,   123,
909      &    123,   953,   953,   123,   405,   794,   123,   647,   613,
910      &   1026,   647,   768,   953,   405,   953,   405,   918,   918,
911      &    123,   953,   953,   918,   953,   536,   405,    70,   124,
912      &   1005,   529,   207,   405,   405,   953,   953,   123,   918,
913      &    918,   953,   405,   918,   953,   468,   405,   794,   794,
914      &    647,   613,   548,   405,   953,   405,   953,   123,   918/
915       DATA P( 9), ( C( 9,I), I = 1, 99 ) /   3079,
916      &   1189,  1423,   287,   186,   341,    77,   733,   733,  1116,
917      &      2,  1539,     2,     2,     2,     2,     2,  1116,   847,
918      &   1174,     2,   827,   713,   910,   944,   139,  1174,  1174,
919      &   1539,  1397,  1397,  1174,   370,    33,  1210,     2,   370,
920      &   1423,   370,   370,  1423,  1423,  1423,   434,  1423,   901,
921      &    139,  1174,   427,   427,   200,  1247,   114,   114,  1441,
922      &    139,   728,  1116,  1174,   139,   113,   113,   113,  1406,
923      &   1247,   200,   200,   200,   200,  1247,  1247,    27,   427,
924      &    427,  1122,  1122,   696,   696,   427,  1539,   435,  1122,
925      &    758,  1247,  1247,  1247,   200,   200,   200,  1247,   114,
926      &     27,   118,   118,   113,   118,   453,   453,  1084,  1406/
927       DATA P(10), ( C(10,I), I = 1, 99 ) /   4621,
928      &   1764,  1349,  1859,   693,    78,   438,   531,    68,  2234,
929      &   2310,  2310,  2310,     2,  2310,  2310,  2102,  2102,   178,
930      &    314,   921,  1074,  1074,  1074,  2147,   314,  1869,   178,
931      &    178,  1324,  1324,   510,  2309,  1541,  1541,  1541,  1541,
932      &    342,  1324,  1324,  1324,  1324,   510,   570,   570,  2197,
933      &    173,  1202,   998,  1324,  1324,   178,  1324,  1324,  1541,
934      &   1541,  1541,   342,  1541,   886,   178,  1324,  1324,  1324,
935      &    510,   784,   784,   501,   652,  1541,  1541,  1324,   178,
936      &   1324,   178,  1324,  1541,   342,  1541,  2144,   784,  2132,
937      &   1324,  1324,  1324,  1324,   510,   652,  1804,  1541,  1541,
938      &   1541,  2132,  1324,  1324,  1324,   178,   510,  1541,   652/
939       DATA P(11), ( C(11,I), I = 1, 99 ) /   6947,
940      &   2872,  1238,   387,  2135,   235,  1565,   221,  1515,  2950,
941      &    486,  3473,     2,  2950,   982,  2950,  3122,  2950,  3172,
942      &   2091,  2091,     9,  3449,  3122,  2846,  3122,  3122,  1947,
943      &   2846,  3122,   772,  1387,  2895,  1387,     3,     3,     3,
944      &   1320,  1320,  2963,  2963,  1320,  1320,  2380,   108,  1284,
945      &    702,  1429,   907,  3220,  3125,  1320,  2963,  1320,  1320,
946      &   2963,  1320,  1639,  3168,  1660,  2895,  2895,  2895,  2895,
947      &   1639,  1297,  1639,   404,  3168,  2963,  2943,  2943,   550,
948      &   1387,  1387,  2895,  2895,  2895,  1387,  2895,  1387,  2895,
949      &   1320,  1320,  2963,  1320,  1320,  1320,  2963,  1320,     2,
950      &   3473,     2,  3473,   772,  2550,     9,  1320,  2963,  1320/
951       DATA P(12), ( C(12,I), I = 1, 99 ) /  10427,
952      &   4309,  2339,  4154,  4480,  4967,   630,  5212,  2592,  4715,
953      &   1808,  1808,  5213,     2,   216,  4014,  3499,  3499,  4204,
954      &   2701,  2701,  5213,  4157,  1209,  4157,  4460,   335,  4460,
955      &   1533,  4575,  4013,  4460,  1881,  2701,  4030,  4030,  1881,
956      &   4030,  1738,   249,   335,    57,  2561,  2561,  2561,  1533,
957      &   1533,  1533,  4013,  4013,  4013,  4013,  4013,  1533,   856,
958      &    856,   468,   468,   468,  2561,   468,  2022,  2022,  2434,
959      &    138,  4605,  1100,  2561,  2561,    57,    57,  3249,   468,
960      &    468,   468,    57,   468,  1738,   313,   856,     6,  3877,
961      &    468,   557,   468,    57,   468,  4605,  2022,     2,  4605,
962      &    138,  1100,    57,  2561,    57,    57,  2022,  5213,  3249/
963       DATA P(13), ( C(13,I), I = 1, 99 ) /  15641,
964      &   6610,  1658,  3022,  2603,  5211,   265,  4985,     3,  4971,
965      &   2127,  1877,  1877,     2,  2925,  3175,  3878,  1940,  1940,
966      &   1940,  5117,  5117,  5771,  5117,  5117,  5117,  5117,  5117,
967      &   5771,  5771,  5117,  3658,  3658,  3658,  3658,  3658,  3658,
968      &   5255,  2925,  2619,  1714,  4100,  6718,  6718,  4100,  2322,
969      &    842,  4100,  6718,  5119,  4728,  5255,  5771,  5771,  5771,
970      &   5117,  5771,  5117,  5117,  5117,  5117,  5117,  5117,  5771,
971      &   5771,  1868,  4483,  4728,  3658,  5255,  3658,  5255,  3658,
972      &   3658,  5255,  5255,  3658,  6718,  6718,   842,  2322,  6718,
973      &   4100,  6718,  4100,  4100,  5117,  5771,  5771,  5117,  5771,
974      &   5771,  5771,  5771,  5117,  5117,  5117,  5771,  5771,  1868/
975       DATA P(14), ( C(14,I), I = 1, 99 ) /  23473,
976      &   9861,  7101,  6257,  7878, 11170, 11638,  7542,  2592,  2591,
977      &   6074,  1428,  8925, 11736,  8925,  5623,  5623,  1535,  6759,
978      &   9953,  9953, 11459,  9953,  7615,  7615, 11377, 11377,  2762,
979      &  11734, 11459,  6892,  1535,  6759,  4695,  1535,  6892,     2,
980      &      2,  6892,  6892,  4177,  4177,  6339,  6950,  1226,  1226,
981      &   1226,  4177,  6892,  6890,  3640,  3640,  1226, 10590, 10590,
982      &   6950,  6950,  6950,  1226,  6950,  6950,  7586,  7586,  7565,
983      &   7565,  3640,  3640,  6950,  7565,  6950,  3599,  3599,  3599,
984      &   2441,  4885,  4885,  4885,  7565,  7565,  1226,  1226,  1226,
985      &   6950,  7586,  1346,  2441,  6339,  3640,  6950, 10590,  6339,
986      &   6950,  6950,  6950,  1226,  1226,  6950,   836,  6891,  7565/
987       DATA P(15), ( C(15,I), I = 1, 99 ) /  35221,
988      &  13482,  5629,  6068, 11974,  4732, 14946, 12097, 17609, 11740,
989      &  15170, 10478, 10478, 17610,     2,     2,  7064,  7064,  7064,
990      &   5665,  1771,  2947,  4453, 12323, 17610, 14809, 14809,  5665,
991      &   5665,  2947,  2947,  2947,  2947, 12323, 12323,  4453,  4453,
992      &   2026, 11772,  2026, 11665, 12323, 12323,  3582,  2940,  2940,
993      &   6654,  4449,  9254, 11470,   304,   304, 11470,   304, 11470,
994      &   6156,  9254, 11772,  6654, 11772,  6156, 11470, 11470, 11772,
995      &  11772, 11772, 11470, 11470,   304, 11470, 11470,   304, 11470,
996      &    304, 11470,   304,   304,   304,  6654, 11508,   304,   304,
997      &   6156,  3582, 11470, 11470, 11470, 17274,  6654,  6654,  6744,
998      &   6711,  6654,  6156,  3370,  6654, 12134,  3370,  6654,  3582/
999       DATA P(16), ( C(16,I), I = 1, 99 ) /  52837,
1000      &  13482,  5629,  6068, 11974,  4732, 14946, 12097, 17609, 11740,
1001      &  15170, 10478, 10478, 17610,     2,     2,  7064,  7064,  7064,
1002      &   5665,  1771,  2947,  4453, 12323, 17610, 14809, 14809,  5665,
1003      &   5665,  2947,  2947,  2947,  2947, 12323, 12323,  4453,  4453,
1004      &   2026, 11772,  2026, 11665, 12323, 12323,  3582,  2940,  2940,
1005      &   6654,  4449,  9254, 11470,   304,   304, 11470,   304, 11470,
1006      &   6156,  9254, 11772,  6654, 11772,  6156, 11470, 11470, 11772,
1007      &  11772, 11772, 11470, 11470,   304, 11470, 11470,   304, 11470,
1008      &    304, 11470,   304,   304,   304,  6654, 11508,   304,   304,
1009      &   6156,  3582, 11470, 11470, 11470, 17274,  6654,  6654,  6744,
1010      &   6711,  6654,  6156,  3370,  6654, 12134,  3370,  6654,  3582/
1011       DATA P(17), ( C(17,I), I = 1, 99 ) /  79259,
1012      &  34566, 38838, 23965, 17279, 35325, 33471,   330, 36050, 26419,
1013      &   3012, 38428, 36430, 36430, 36755, 39629,  5749,  5749, 36755,
1014      &   5749, 14353, 14353, 14353, 32395, 32395, 32395, 32395, 32396,
1015      &  32396, 32396, 32396, 27739, 14353, 36430, 36430, 36430, 15727,
1016      &  38428, 28987, 28987, 27739, 38428, 27739, 18786, 14353, 15727,
1017      &  28987, 19151, 19757, 19757, 19757, 14353, 22876, 19151, 24737,
1018      &  24737,  4412, 30567, 30537, 19757, 30537, 19757, 30537, 30537,
1019      &   4412, 24737, 28987, 19757, 19757, 19757, 30537, 30537, 33186,
1020      &   4010,  4010,  4010, 17307, 15217, 32789, 37709,  4010,  4010,
1021      &   4010, 33186, 33186,  4010, 11057, 39388, 33186,  1122, 15089,
1022      &  39629,     2,     2, 23899, 16466, 16466, 17038,  9477,  9260/
1023       DATA P(18), ( C(18,I), I = 1, 99 ) / 118891,
1024      &  31929, 40295,  2610,  5177, 17271, 23770,  9140,   952, 39631,
1025      &      3, 11424, 49719, 38267, 25172,     2,     2, 59445,     2,
1026      &  59445, 38267, 44358, 14673, 53892, 14674, 14673, 14674, 41368,
1027      &  17875, 17875, 30190, 20444, 55869, 15644, 25499, 15644, 20983,
1028      &  44358, 15644, 15644,   485, 41428,   485,   485,   485, 41428,
1029      &  53798, 50230, 53798, 50253, 50253, 35677, 35677, 17474,  7592,
1030      &   4098, 17474,   485, 41428,   485, 41428,   485, 41428,   485,
1031      &  41428, 41428, 41428, 41428, 41428,  9020, 22816,  4098,  4098,
1032      &   4098,  7592, 42517,   485, 50006, 50006, 22816, 22816,  9020,
1033      &    485, 41428, 41428, 41428, 41428, 50006,   485, 41428, 41428,
1034      &  41428, 41428, 22816, 41428, 41428,   485,   485,   485,  9020/
1035       DATA P(19), ( C(19,I), I = 1, 99 ) / 178349,
1036      &  73726, 16352, 16297, 74268, 60788,  8555,  1077, 25486, 86595,
1037      &  59450, 19958, 62205, 62205,  4825,  4825, 89174, 89174, 62205,
1038      &  19958, 62205, 19958, 27626, 63080, 62205, 62205, 62205, 19958,
1039      &   8914, 83856, 30760, 47774, 47774, 19958, 62205, 39865, 39865,
1040      &  74988, 75715, 75715, 74988, 34522, 74988, 74988, 25101, 44621,
1041      &  44621, 44621, 25101, 25101, 25101, 44621, 47768, 41547, 44621,
1042      &  10273, 74988, 74988, 74988, 74988, 74988, 74988, 34522, 34522,
1043      &  67796, 67796, 30208,     2, 67062, 18500, 29251, 29251,     2,
1044      &  67796, 67062, 38649, 59302,  6225, 67062,  6475,  6225, 46772,
1045      &  38649, 67062, 46772, 46772, 67062, 46772, 25372, 67062,  6475,
1046      &  25372, 67062, 67062, 67062,  6225, 67062, 67062, 68247, 80676/
1047       DATA P(20), ( C(20,I), I = 1, 99 )/ 267523,
1048      & 103650, 50089, 70223, 41805, 74847,112775, 40889, 64866, 44053,
1049      &   1754,129471, 13630, 53467, 53467, 61378,133761,     2,133761,
1050      &      2,133761,133761, 65531, 65531, 65531, 38080,133761,133761,
1051      & 131061,  5431, 65531, 78250, 11397, 38841, 38841,107233,107233,
1052      & 111286, 19065, 38841, 19065, 19065, 16099,127638, 82411, 96659,
1053      &  96659, 82411, 96659, 82411, 51986,101677, 39264, 39264,101677,
1054      &  39264, 39264, 47996, 96659, 82411, 47996, 10971, 10004, 82411,
1055      &  96659, 82411, 82411, 82411, 96659, 96659, 96659, 82411, 96659,
1056      &  51986,110913, 51986, 51986,110913, 82411, 54713, 54713, 22360,
1057      & 117652, 22360, 78250, 78250, 91996, 22360, 91996, 97781, 91996,
1058      &  97781, 91996, 97781, 97781, 91996, 97781, 97781, 36249, 39779/
1059       END
1060 *
1061       SUBROUTINE KROSUM( NDIM, SUMKRO, PRIME, VK, FUNCTN, ALPHA, X )
1062       EXTERNAL FUNCTN
1063       INTEGER NDIM, PRIME, K, J
1064       DOUBLE PRECISION SUMKRO, VK(*), FUNCTN, ALPHA(*), X(*), ONE, UNI
1065       PARAMETER ( ONE = 1 )
1066       SUMKRO = 0
1067       DO J = 1, NDIM
1068          ALPHA(J) = UNI()
1069       END DO
1070       DO K = 1, PRIME
1071          DO J = 1, NDIM
1072             X(J) = MOD( K*VK(J) + ALPHA(J), ONE )
1073             X(J) = ABS( 2*X(J) - 1 )
1074          END DO
1075          SUMKRO = SUMKRO + ( FUNCTN(NDIM,X) - SUMKRO )/( 2*K - 1 )
1076          DO J = 1, NDIM
1077             X(J) = 1 - X(J)
1078          END DO
1079          SUMKRO = SUMKRO + ( FUNCTN(NDIM,X) - SUMKRO )/( 2*K )
1080       END DO
1081       END
1082 *
1083       DOUBLE PRECISION FUNCTION MVNFNC(N, W)
1084 *
1085 *     Integrand subroutine
1086 *
1087       INTEGER N, INFIN(*), INFIS
1088       DOUBLE PRECISION W(*), LOWER(*), UPPER(*), CORREL(*), ONE
1089       INTEGER NL, IJ, I, J
1090       PARAMETER ( NL = 100, ONE = 1 )
1091       DOUBLE PRECISION COV((NL*(NL+1))/2), A(NL), B(NL), Y(NL), BVN
1092       INTEGER INFI(NL)
1093       DOUBLE PRECISION PROD, D1, E1, DI, EI, SUM, PHINV, D, E, MVNNIT
1094       SAVE D1, E1, A, B, INFI, COV
1095       DI = D1
1096       EI = E1
1097       PROD = E1 - D1
1098       IJ = 1
1099       DO I = 1,N
1100          Y(I) = PHINV( DI + W(I)*(EI-DI) )
1101          SUM = 0
1102          DO J = 1,I
1103             IJ = IJ + 1
1104             SUM = SUM + COV(IJ)*Y(J)
1105          END DO
1106          IJ = IJ + 1
1107          IF ( COV(IJ) .GT. 0 ) THEN
1108             CALL LIMITS( A(I+1)-SUM, B(I+1)-SUM, INFI(I+1), DI, EI )
1109          ELSE
1110             DI = ( 1 + SIGN( ONE, A(I+1)-SUM ) )/2
1111             EI = ( 1 + SIGN( ONE, B(I+1)-SUM ) )/2
1112          ENDIF
1113          PROD = PROD*(EI-DI)
1114       END DO
1115       MVNFNC = PROD
1116       RETURN
1117 *
1118 *     Entry point for intialization.
1119 *
1120       ENTRY MVNNIT(N, CORREL, LOWER, UPPER, INFIN, INFIS, D, E)
1121       MVNNIT = 0
1122 *
1123 *     Initialization and computation of covariance Cholesky factor.
1124 *
1125       CALL NCVSRT(N, LOWER,UPPER,CORREL,INFIN,Y, INFIS,A,B,INFI,COV,D,E)
1126       D1 = D
1127       E1 = E
1128       IF ( N - INFIS .EQ. 2 ) THEN
1129          D = SQRT( 1 + COV(2)**2 )
1130          A(2) = A(2)/D
1131          B(2) = B(2)/D
1132          E = BVN( A, B, INFI, COV(2)/D )
1133          D = 0
1134          INFIS = INFIS + 1
1135       END IF
1136       END
1137       SUBROUTINE LIMITS( A, B, INFIN, LOWER, UPPER )
1138       DOUBLE PRECISION A, B, LOWER, UPPER, PHI
1139       INTEGER INFIN
1140       LOWER = 0
1141       UPPER = 1
1142       IF ( INFIN .GE. 0 ) THEN
1143          IF ( INFIN .NE. 0 ) LOWER = PHI(A)
1144          IF ( INFIN .NE. 1 ) UPPER = PHI(B)
1145       ENDIF
1146       END
1147       SUBROUTINE NCVSRT( N, LOWER, UPPER, CORREL, INFIN, Y, INFIS,
1148      &                   A, B, INFI, COV, D, E )
1149 *
1150 *     Subroutine to sort integration limits.
1151 *
1152       INTEGER N, INFI(*), INFIN(*), INFIS
1153       DOUBLE PRECISION
1154      &     A(*), B(*), COV(*), LOWER(*), UPPER(*), CORREL(*), Y(*), D, E
1155       INTEGER I, J, K, IJ, II, JMIN
1156       DOUBLE PRECISION SUMSQ, ZERO
1157       PARAMETER ( ZERO = 0 )
1158       DOUBLE PRECISION AJ, BJ, SUM, SQTWPI
1159       DOUBLE PRECISION CVDIAG, AMIN, BMIN, DMIN, EMIN, YL, YU
1160       PARAMETER ( SQTWPI = 2.50662 82746 31000 50240 )
1161       IJ = 0
1162       II = 0
1163       INFIS = 0
1164       DO I = 1,N
1165          INFI(I) = INFIN(I)
1166          IF ( INFI(I) .LT. 0 ) THEN
1167             INFIS = INFIS + 1
1168          ELSE
1169             A(I) = 0
1170             B(I) = 0
1171             IF ( INFI(I) .NE. 0 ) A(I) = LOWER(I)
1172             IF ( INFI(I) .NE. 1 ) B(I) = UPPER(I)
1173          ENDIF
1174          DO J = 1,I-1
1175             IJ = IJ + 1
1176             II = II + 1
1177             COV(IJ) = CORREL(II)
1178          END DO
1179          IJ = IJ + 1
1180          COV(IJ) = 1
1181       END DO
1182 *
1183 *     First move any doubly infinite limits to innermost positions
1184 *
1185       IF ( INFIS .LT. N ) THEN
1186          outer: DO I = N,N-INFIS+1,-1
1187             IF ( INFI(I) .GE. 0 ) THEN
1188                inner: DO J = 1,I-1
1189                   IF ( INFI(J) .LT. 0 ) THEN
1190                      CALL RCSWAP(J, I, A, B, INFI, N, COV)
1191                      CYCLE outer
1192                   ENDIF
1193                END DO inner
1194             ENDIF
1195          END DO outer
1196 *
1197 *     Sort remaining limits and determine Cholesky decomposition
1198 *
1199          II = 0
1200          DO I = 1,N-INFIS
1201 *
1202 *     Determine the integration limits for variable with minimum
1203 *      expected probability and interchange that variable with Ith.
1204 *
1205             EMIN = 1
1206             DMIN = 0
1207             JMIN = I
1208             CVDIAG = 0
1209             IJ = II
1210             DO J = I, N-INFIS
1211                SUM = 0
1212                SUMSQ = 0
1213                DO K = 1, I-1
1214                   SUM = SUM + COV(IJ+K)*Y(K)
1215                   SUMSQ = SUMSQ + COV(IJ+K)**2
1216                END DO
1217                IJ = IJ + J
1218                SUMSQ = SQRT( MAX( COV(IJ)-SUMSQ, ZERO ) )
1219                IF ( SUMSQ .GT. 0 ) THEN
1220                   IF ( INFI(J) .NE. 0 ) AJ = ( A(J) - SUM )/SUMSQ
1221                   IF ( INFI(J) .NE. 1 ) BJ = ( B(J) - SUM )/SUMSQ
1222                   CALL LIMITS( AJ, BJ, INFI(J), D, E )
1223                   IF ( EMIN - DMIN .GE. E - D ) THEN
1224                      JMIN = J
1225                      IF ( INFI(J) .NE. 0 ) AMIN = AJ
1226                      IF ( INFI(J) .NE. 1 ) BMIN = BJ
1227                      DMIN = D
1228                      EMIN = E
1229                      CVDIAG = SUMSQ
1230                   ENDIF
1231                ENDIF
1232             END DO
1233             IF ( JMIN .NE. I) CALL RCSWAP(I, JMIN, A, B, INFI, N, COV)
1234 *
1235 *     Compute Ith column of Cholesky factor.
1236 *
1237             IJ = II + I
1238             COV(IJ) = CVDIAG
1239             DO J = I+1, N-INFIS
1240                IF ( CVDIAG .GT. 0 ) THEN
1241                   SUM = COV(IJ+I)
1242                   DO K = 1, I-1
1243                      SUM = SUM - COV(II+K)*COV(IJ+K)
1244                   END DO
1245                   COV(IJ+I) = SUM/CVDIAG
1246                ELSE
1247                   COV(IJ+I) = 0
1248                ENDIF
1249                IJ = IJ + J
1250             END DO
1251 *
1252 *     Compute expected value for Ith integration variable and
1253 *     scale Ith covariance matrix row and limits.
1254 *
1255             IF ( CVDIAG .GT. 0 ) THEN
1256                IF ( EMIN .GT. DMIN + 1D-8 ) THEN
1257                   YL = 0
1258                   YU = 0
1259                   IF ( INFI(I) .NE. 0 ) YL = -EXP( -AMIN**2/2 )/SQTWPI
1260                   IF ( INFI(I) .NE. 1 ) YU = -EXP( -BMIN**2/2 )/SQTWPI
1261                   Y(I) = ( YU - YL )/( EMIN - DMIN )
1262                ELSE
1263                   IF ( INFI(I) .EQ. 0 ) Y(I) = BMIN
1264                   IF ( INFI(I) .EQ. 1 ) Y(I) = AMIN
1265                   IF ( INFI(I) .EQ. 2 ) Y(I) = ( AMIN + BMIN )/2
1266                END IF
1267                DO J = 1,I
1268                   II = II + 1
1269                   COV(II) = COV(II)/CVDIAG
1270                END DO
1271                IF ( INFI(I) .NE. 0 ) A(I) = A(I)/CVDIAG
1272                IF ( INFI(I) .NE. 1 ) B(I) = B(I)/CVDIAG
1273             ELSE
1274                Y(I) = 0
1275                II = II + I
1276             ENDIF
1277          END DO
1278          CALL LIMITS( A(1), B(1), INFI(1), D, E)
1279       ENDIF
1280       END
1281       DOUBLE PRECISION FUNCTION CONDIT( N, SYMIN )
1282 *
1283 *     Computes condition number of symmetric matix in situ
1284 *
1285       INTEGER NL, N
1286       PARAMETER ( NL = 100 )
1287       DOUBLE PRECISION DET, SYMIN(*), SUM, ROWMX, ROWMXI,
1288      & SYM(NL*(NL+1)/2)
1289       INTEGER II, IJ, I, J, IM
1290       ROWMX = 0
1291       IJ = 0
1292       DO I = 1,N
1293          SUM = 0
1294          IM = (I-2)*(I-1)/2
1295          DO J = 1,I-1
1296             IM = IM + 1
1297             SUM = SUM + ABS(SYMIN(IM))
1298             IJ = IJ + 1
1299             SYM(IJ) = SYMIN(IM)
1300          END DO
1301          SUM = SUM + 1
1302          IJ = IJ + 1
1303          SYM(IJ) = 1
1304          IM = IM + I
1305          DO J = I,N-1
1306             SUM = SUM + ABS(SYMIN(IM))
1307             IM = IM + J
1308          END DO
1309          ROWMX = MAX( SUM, ROWMX )
1310       END DO
1311       CALL SYMINV2(N, SYM, DET)
1312       ROWMXI = 0
1313       II = 0
1314       DO I = 1,N
1315          SUM = 0
1316          IJ = II
1317          DO J = 1,I
1318             IJ = IJ + 1
1319             SUM = SUM + ABS(SYM(IJ))
1320          END DO
1321          DO J = I,N-1
1322             IJ = IJ + J
1323             SUM = SUM + ABS(SYM(IJ))
1324          END DO
1325          ROWMXI = MAX( SUM, ROWMXI )
1326          II = II + I
1327       END DO
1328       CONDIT = ROWMX*ROWMXI
1329       END
1330       SUBROUTINE SYMINV2(N, LOWINV, DET)
1331 *
1332 *     Computes lower symmetric inverse and determinant in situ
1333 *
1334       INTEGER I, II, N
1335       DOUBLE PRECISION LOWINV(*), DET
1336       CALL CHOLSK(N, LOWINV)
1337       DET = 1
1338       II = 0
1339       DO I = 1,N
1340          II = II + I
1341          DET = DET*LOWINV(II)
1342       END DO
1343       DET = DET*DET
1344       CALL CHOLNV(N, LOWINV)
1345       CALL CHOLPI(N, LOWINV)
1346       END
1347       SUBROUTINE CHOLSK(N, CHOFAC)
1348 *
1349 *     Computes Choleski factor in situ
1350 *
1351       INTEGER I, II, J, JJ, K, N
1352       DOUBLE PRECISION CHOFAC(*), T
1353       DOUBLE PRECISION S, ZERO
1354       PARAMETER ( ZERO = 0 )
1355       JJ = 0
1356       DO J = 1,N
1357          II = JJ
1358          DO I = J,N
1359             S = CHOFAC(II+J)
1360             DO K = 1,J-1
1361                S = S - CHOFAC(II+K)*CHOFAC(JJ+K)
1362             END DO
1363             IF ( I .EQ. J ) THEN
1364                T = SQRT( MAX( S, ZERO ) )
1365                CHOFAC(II+J) = T
1366             ELSE
1367                CHOFAC(II+J) = S/T
1368             ENDIF
1369             II = II + I
1370          END DO
1371          JJ = JJ + J
1372       END DO
1373       END
1374       SUBROUTINE CHOLNV(N, CHOINV)
1375 *
1376 *     Inverts a lower triangular matrix in situ
1377 *
1378       INTEGER I, II, J, JJ, K, KK, N
1379       DOUBLE PRECISION CHOINV(*), T
1380       DOUBLE PRECISION S
1381       II = 0
1382       DO I = 1,N
1383          T = 1/CHOINV(II+I)
1384          JJ = 0
1385          DO J = 1,I-1
1386             S = 0
1387             JJ = JJ + J
1388             KK = JJ
1389             DO K = J,I-1
1390                S = S + CHOINV(II+K)*CHOINV(KK)
1391                KK = KK + K
1392             END DO
1393             CHOINV(II+J) = -S*T
1394          END DO
1395          II = II + I
1396          CHOINV(II) = T
1397       END DO
1398       END
1399       SUBROUTINE CHOLPI(N, CHOPDI)
1400 *
1401 *     Multiplies Choleski inverse factors in situ
1402 *
1403       INTEGER I, II, J, JJ, K, KK, N
1404       DOUBLE PRECISION CHOPDI(*)
1405       DOUBLE PRECISION S
1406       II = 0
1407       DO I = 1,N
1408          DO J = 1,I
1409             S = 0
1410             JJ = II + I
1411             KK = II + J
1412             DO K = I,N
1413                S = S + CHOPDI(KK)*CHOPDI(JJ)
1414                JJ = JJ + K
1415                KK = KK + K
1416             END DO
1417             CHOPDI(II+J) = S
1418          END DO
1419          II = II + I
1420       END DO
1421       END
1422       SUBROUTINE RCSWAP(P, Q, A, B, INFIN, N, C)
1423 *
1424 *     Swaps rows and columns P and Q in situ.
1425 *
1426       DOUBLE PRECISION A(*), B(*), C(*), T
1427       INTEGER INFIN(*), P, Q, N, I, J, II, JJ
1428       T = A(P)
1429       A(P) = A(Q)
1430       A(Q) = T
1431       T = B(P)
1432       B(P) = B(Q)
1433       B(Q) = T
1434       J = INFIN(P)
1435       INFIN(P) = INFIN(Q)
1436       INFIN(Q) = J
1437       JJ = (P*(P-1))/2
1438       II = (Q*(Q-1))/2
1439       T = C(JJ+P)
1440       C(JJ+P) = C(II+Q)
1441       C(II+Q) = T
1442       DO J = 1, P-1
1443          T = C(JJ+J)
1444          C(JJ+J) = C(II+J)
1445          C(II+J) = T
1446       END DO
1447       JJ = JJ + P
1448       DO I = P+1, Q-1
1449          T = C(JJ+P)
1450          C(JJ+P) = C(II+I)
1451          C(II+I) = T
1452          JJ = JJ + I
1453       END DO
1454       II = II + Q
1455       DO I = Q+1, N
1456          T = C(II+P)
1457          C(II+P) = C(II+Q)
1458          C(II+Q) = T
1459          II = II + I
1460       END DO
1461       END
1462       DOUBLE PRECISION FUNCTION PHI(Z)
1463 *
1464 *     Normal distribution probabilities accurate to 1.e-15.
1465 *     Z = no. of standard deviations from the mean.
1466 *
1467 *     Based upon algorithm 5666 for the error function, from:
1468 *     Hart, J.F. et al, 'Computer Approximations', Wiley 1968
1469 *
1470 *     Programmer: Alan Miller
1471 *
1472 *     Latest revision - 30 March 1986
1473 *
1474       DOUBLE PRECISION P0, P1, P2, P3, P4, P5, P6,
1475      &     Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7,
1476      &     Z, P, EXPNTL, CUTOFF, ROOTPI, ZABS
1477       PARAMETER(
1478      &     P0 = 220.20 68679 12376 1D0,
1479      &     P1 = 221.21 35961 69931 1D0,
1480      &     P2 = 112.07 92914 97870 9D0,
1481      &     P3 = 33.912 86607 83830 0D0,
1482      &     P4 = 6.3739 62203 53165 0D0,
1483      &     P5 = .70038 30644 43688 1D0,
1484      &     P6 = .035262 49659 98910 9D0)
1485       PARAMETER(
1486      &     Q0 = 440.41 37358 24752 2D0,
1487      &     Q1 = 793.82 65125 19948 4D0,
1488      &     Q2 = 637.33 36333 78831 1D0,
1489      &     Q3 = 296.56 42487 79673 7D0,
1490      &     Q4 = 86.780 73220 29460 8D0,
1491      &     Q5 = 16.064 17757 92069 5D0,
1492      &     Q6 = 1.7556 67163 18264 2D0,
1493      &     Q7 = .088388 34764 83184 4D0)
1494       PARAMETER(ROOTPI = 2.5066 28274 63100 1D0)
1495       PARAMETER(CUTOFF = 7.0710 67811 86547 5D0)
1496 *
1497       ZABS = ABS(Z)
1498 *
1499 *     |Z| > 37
1500 *
1501       IF (ZABS .GT. 37) THEN
1502          P = 0
1503       ELSE
1504 *
1505 *     |Z| <= 37
1506 *
1507          EXPNTL = EXP(-ZABS**2/2)
1508 *
1509 *     |Z| < CUTOFF = 10/SQRT(2)
1510 *
1511          IF (ZABS .LT. CUTOFF) THEN
1512             P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS
1513      &           + P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS
1514      &           + Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS
1515      &           + Q0)
1516 *
1517 *     |Z| >= CUTOFF.
1518 *
1519          ELSE
1520             P = EXPNTL/(ZABS + 1/(ZABS + 2/(ZABS + 3/(ZABS + 4/
1521      &           (ZABS + 0.65D0)))))/ROOTPI
1522          END IF
1523       END IF
1524       IF (Z .GT. 0) P = 1 - P
1525       PHI = P
1526       END
1527       DOUBLE PRECISION FUNCTION PHINV(P)
1528 *
1529 *       ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
1530 *
1531 *       Produces the normal deviate Z corresponding to a given lower
1532 *       tail area of P.
1533 *
1534 *       The hash sums below are the sums of the mantissas of the
1535 *       coefficients.   They are included for use in checking
1536 *       transcription.
1537 *
1538       DOUBLE PRECISION SPLIT1, SPLIT2, CONST1, CONST2,
1539      &     A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, B4, B5, B6, B7,
1540      &     C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, D6, D7,
1541      &     E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, F4, F5, F6, F7,
1542      &     P, Q, R
1543       PARAMETER (SPLIT1 = 0.425, SPLIT2 = 5,
1544      &     CONST1 = 0.180625D0, CONST2 = 1.6D0)
1545 *
1546 *     Coefficients for P close to 0.5
1547 *
1548       PARAMETER (
1549      &     A0 = 3.38713 28727 96366 6080D0,
1550      &     A1 = 1.33141 66789 17843 7745D+2,
1551      &     A2 = 1.97159 09503 06551 4427D+3,
1552      &     A3 = 1.37316 93765 50946 1125D+4,
1553      &     A4 = 4.59219 53931 54987 1457D+4,
1554      &     A5 = 6.72657 70927 00870 0853D+4,
1555      &     A6 = 3.34305 75583 58812 8105D+4,
1556      &     A7 = 2.50908 09287 30122 6727D+3,
1557      &     B1 = 4.23133 30701 60091 1252D+1,
1558      &     B2 = 6.87187 00749 20579 0830D+2,
1559      &     B3 = 5.39419 60214 24751 1077D+3,
1560      &     B4 = 2.12137 94301 58659 5867D+4,
1561      &     B5 = 3.93078 95800 09271 0610D+4,
1562      &     B6 = 2.87290 85735 72194 2674D+4,
1563      &     B7 = 5.22649 52788 52854 5610D+3)
1564 *     HASH SUM AB    55.88319 28806 14901 4439
1565 *
1566 *     Coefficients for P not close to 0, 0.5 or 1.
1567 *
1568       PARAMETER (
1569      &     C0 = 1.42343 71107 49683 57734D0,
1570      &     C1 = 4.63033 78461 56545 29590D0,
1571      &     C2 = 5.76949 72214 60691 40550D0,
1572      &     C3 = 3.64784 83247 63204 60504D0,
1573      &     C4 = 1.27045 82524 52368 38258D0,
1574      &     C5 = 2.41780 72517 74506 11770D-1,
1575      &     C6 = 2.27238 44989 26918 45833D-2,
1576      &     C7 = 7.74545 01427 83414 07640D-4,
1577      &     D1 = 2.05319 16266 37758 82187D0,
1578      &     D2 = 1.67638 48301 83803 84940D0,
1579      &     D3 = 6.89767 33498 51000 04550D-1,
1580      &     D4 = 1.48103 97642 74800 74590D-1,
1581      &     D5 = 1.51986 66563 61645 71966D-2,
1582      &     D6 = 5.47593 80849 95344 94600D-4,
1583      &     D7 = 1.05075 00716 44416 84324D-9)
1584 *     HASH SUM CD    49.33206 50330 16102 89036
1585 *
1586 *       Coefficients for P near 0 or 1.
1587 *
1588       PARAMETER (
1589      &     E0 = 6.65790 46435 01103 77720D0,
1590      &     E1 = 5.46378 49111 64114 36990D0,
1591      &     E2 = 1.78482 65399 17291 33580D0,
1592      &     E3 = 2.96560 57182 85048 91230D-1,
1593      &     E4 = 2.65321 89526 57612 30930D-2,
1594      &     E5 = 1.24266 09473 88078 43860D-3,
1595      &     E6 = 2.71155 55687 43487 57815D-5,
1596      &     E7 = 2.01033 43992 92288 13265D-7,
1597      &     F1 = 5.99832 20655 58879 37690D-1,
1598      &     F2 = 1.36929 88092 27358 05310D-1,
1599      &     F3 = 1.48753 61290 85061 48525D-2,
1600      &     F4 = 7.86869 13114 56132 59100D-4,
1601      &     F5 = 1.84631 83175 10054 68180D-5,
1602      &     F6 = 1.42151 17583 16445 88870D-7,
1603      &     F7 = 2.04426 31033 89939 78564D-15)
1604 *     HASH SUM EF    47.52583 31754 92896 71629
1605 *
1606       Q = ( 2*P - 1 )/2
1607       IF ( ABS(Q) .LE. SPLIT1 ) THEN
1608          R = CONST1 - Q*Q
1609          PHINV = Q*(((((((A7*R + A6)*R + A5)*R + A4)*R + A3)
1610      &        *R + A2)*R + A1)*R + A0) /
1611      &        (((((((B7*R + B6)*R + B5)*R + B4)*R + B3)
1612      &        *R + B2)*R + B1)*R + 1)
1613       ELSE
1614          R = MIN( P, 1 - P )
1615          IF (R .GT. 0) THEN
1616             R = SQRT( -LOG(R) )
1617             IF ( R .LE. SPLIT2 ) THEN
1618                R = R - CONST2
1619                PHINV = (((((((C7*R + C6)*R + C5)*R + C4)*R + C3)
1620      &              *R + C2)*R + C1)*R + C0) /
1621      &              (((((((D7*R + D6)*R + D5)*R + D4)*R + D3)
1622      &              *R + D2)*R + D1)*R + 1)
1623             ELSE
1624                R = R - SPLIT2
1625                PHINV = (((((((E7*R + E6)*R + E5)*R + E4)*R + E3)
1626      &              *R + E2)*R + E1)*R + E0) /
1627      &              (((((((F7*R + F6)*R + F5)*R + F4)*R + F3)
1628      &              *R + F2)*R + F1)*R + 1)
1629             END IF
1630          ELSE
1631             PHINV = 9
1632          END IF
1633          IF ( Q .LT. 0 ) PHINV = - PHINV
1634       END IF
1635       END
1636       DOUBLE PRECISION FUNCTION BVN ( LOWER, UPPER, INFIN, CORREL )
1637 *
1638 *     A function for computing bivariate normal probabilities.
1639 *
1640 *  Parameters
1641 *
1642 *     LOWER  REAL, array of lower integration limits.
1643 *     UPPER  REAL, array of upper integration limits.
1644 *     INFIN  INTEGER, array of integration limits flags:
1645 *            if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)];
1646 *            if INFIN(I) = 1, Ith limits are [LOWER(I), infinity);
1647 *            if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)].
1648 *     CORREL REAL, correlation coefficient.
1649 *
1650       DOUBLE PRECISION LOWER(*), UPPER(*), CORREL, BVNU
1651       INTEGER INFIN(*)
1652       IF ( INFIN(1) .EQ. 2  .AND. INFIN(2) .EQ. 2 ) THEN
1653          BVN =  BVNU ( LOWER(1), LOWER(2), CORREL )
1654      +        - BVNU ( UPPER(1), LOWER(2), CORREL )
1655      +        - BVNU ( LOWER(1), UPPER(2), CORREL )
1656      +        + BVNU ( UPPER(1), UPPER(2), CORREL )
1657       ELSE IF ( INFIN(1) .EQ. 2  .AND. INFIN(2) .EQ. 1 ) THEN
1658          BVN =  BVNU ( LOWER(1), LOWER(2), CORREL )
1659      +        - BVNU ( UPPER(1), LOWER(2), CORREL )
1660       ELSE IF ( INFIN(1) .EQ. 1  .AND. INFIN(2) .EQ. 2 ) THEN
1661          BVN =  BVNU ( LOWER(1), LOWER(2), CORREL )
1662      +        - BVNU ( LOWER(1), UPPER(2), CORREL )
1663       ELSE IF ( INFIN(1) .EQ. 2  .AND. INFIN(2) .EQ. 0 ) THEN
1664          BVN =  BVNU ( -UPPER(1), -UPPER(2), CORREL )
1665      +        - BVNU ( -LOWER(1), -UPPER(2), CORREL )
1666       ELSE IF ( INFIN(1) .EQ. 0  .AND. INFIN(2) .EQ. 2 ) THEN
1667          BVN =  BVNU ( -UPPER(1), -UPPER(2), CORREL )
1668      +        - BVNU ( -UPPER(1), -LOWER(2), CORREL )
1669       ELSE IF ( INFIN(1) .EQ. 1  .AND. INFIN(2) .EQ. 0 ) THEN
1670          BVN =  BVNU ( LOWER(1), -UPPER(2), -CORREL )
1671       ELSE IF ( INFIN(1) .EQ. 0  .AND. INFIN(2) .EQ. 1 ) THEN
1672          BVN =  BVNU ( -UPPER(1), LOWER(2), -CORREL )
1673       ELSE IF ( INFIN(1) .EQ. 1  .AND. INFIN(2) .EQ. 1 ) THEN
1674          BVN =  BVNU ( LOWER(1), LOWER(2), CORREL )
1675       ELSE IF ( INFIN(1) .EQ. 0  .AND. INFIN(2) .EQ. 0 ) THEN
1676          BVN =  BVNU ( -UPPER(1), -UPPER(2), CORREL )
1677       END IF
1678       END
1679       DOUBLE PRECISION FUNCTION BVNU( SH, SK, R )
1680 *
1681 *     A function for computing bivariate normal probabilities.
1682 *
1683 *       Yihong Ge
1684 *       Department of Computer Science and Electrical Engineering
1685 *       Washington State University
1686 *       Pullman, WA 99164-2752
1687 *       Email : yge@eecs.wsu.edu
1688 *     and
1689 *       Alan Genz
1690 *       Department of Mathematics
1691 *       Washington State University
1692 *       Pullman, WA 99164-3113
1693 *       Email : alangenz@wsu.edu
1694 *
1695 * BVN - calculate the probability that X is larger than SH and Y is
1696 *       larger than SK.
1697 *
1698 * Parameters
1699 *
1700 *   SH  REAL, integration limit
1701 *   SK  REAL, integration limit
1702 *   R   REAL, correlation coefficient
1703 *   LG  INTEGER, number of Gauss Rule Points and Weights
1704 *
1705       DOUBLE PRECISION BVN, SH, SK, R, ZERO, TWOPI
1706       INTEGER I, LG, NG
1707       PARAMETER ( ZERO = 0, TWOPI = 6.2831 85307 179586 )
1708       DOUBLE PRECISION X(10,3), W(10,3), AS, A, B, C, D, RS, XS
1709       DOUBLE PRECISION PHI, SN, ASR, H, K, BS, HS, HK
1710 *     Gauss Legendre Points and Weights, N =  6
1711       DATA ( W(I,1), X(I,1), I = 1,3) /
1712      &  0.1713244923791705D+00,-0.9324695142031522D+00,
1713      &  0.3607615730481384D+00,-0.6612093864662647D+00,
1714      &  0.4679139345726904D+00,-0.2386191860831970D+00/
1715 *     Gauss Legendre Points and Weights, N = 12
1716       DATA ( W(I,2), X(I,2), I = 1,6) /
1717      &  0.4717533638651177D-01,-0.9815606342467191D+00,
1718      &  0.1069393259953183D+00,-0.9041172563704750D+00,
1719      &  0.1600783285433464D+00,-0.7699026741943050D+00,
1720      &  0.2031674267230659D+00,-0.5873179542866171D+00,
1721      &  0.2334925365383547D+00,-0.3678314989981802D+00,
1722      &  0.2491470458134029D+00,-0.1252334085114692D+00/
1723 *     Gauss Legendre Points and Weights, N = 20
1724       DATA ( W(I,3), X(I,3), I = 1,10) /
1725      &  0.1761400713915212D-01,-0.9931285991850949D+00,
1726      &  0.4060142980038694D-01,-0.9639719272779138D+00,
1727      &  0.6267204833410906D-01,-0.9122344282513259D+00,
1728      &  0.8327674157670475D-01,-0.8391169718222188D+00,
1729      &  0.1019301198172404D+00,-0.7463319064601508D+00,
1730      &  0.1181945319615184D+00,-0.6360536807265150D+00,
1731      &  0.1316886384491766D+00,-0.5108670019508271D+00,
1732      &  0.1420961093183821D+00,-0.3737060887154196D+00,
1733      &  0.1491729864726037D+00,-0.2277858511416451D+00,
1734      &  0.1527533871307259D+00,-0.7652652113349733D-01/
1735       SAVE X, W
1736       IF ( ABS(R) .LT. 0.3 ) THEN
1737          NG = 1
1738          LG = 3
1739       ELSE IF ( ABS(R) .LT. 0.75 ) THEN
1740          NG = 2
1741          LG = 6
1742       ELSE
1743          NG = 3
1744          LG = 10
1745       ENDIF
1746       H = SH
1747       K = SK
1748       HK = H*K
1749       BVN = 0
1750       IF ( ABS(R) .LT. 0.925 ) THEN
1751          HS = ( H*H + K*K )/2
1752          ASR = ASIN(R)
1753          DO 10  I = 1, LG
1754             SN = SIN(ASR*( X(I,NG)+1 )/2)
1755             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
1756             SN = SIN(ASR*(-X(I,NG)+1 )/2)
1757             BVN = BVN + W(I,NG)*EXP( ( SN*HK - HS )/( 1 - SN*SN ) )
1758  10      CONTINUE
1759          BVN = BVN*ASR/(2*TWOPI) + PHI(-H)*PHI(-K)
1760       ELSE
1761          IF ( R .LT. 0 ) THEN
1762             K = -K
1763             HK = -HK
1764          ENDIF
1765          IF ( ABS(R) .LT. 1 ) THEN
1766             AS = ( 1 - R )*( 1 + R )
1767             A = SQRT(AS)
1768             BS = ( H - K )**2
1769             C = ( 4 - HK )/8
1770             D = ( 12 - HK )/16
1771             BVN = A*EXP( -(BS/AS + HK)/2 )
1772      +             *( 1 - C*(BS - AS)*(1 - D*BS/5)/3 + C*D*AS*AS/5 )
1773             IF ( HK .GT. -160 ) THEN
1774                B = SQRT(BS)
1775                BVN = BVN - EXP(-HK/2)*SQRT(TWOPI)*PHI(-B/A)*B
1776      +                    *( 1 - C*BS*( 1 - D*BS/5 )/3 )
1777             ENDIF
1778             A = A/2
1779             DO 20 I = 1, LG
1780                XS = ( A*(X(I,NG)+1) )**2
1781                RS = SQRT( 1 - XS )
1782                BVN = BVN + A*W(I,NG)*
1783      +              ( EXP( -BS/(2*XS) - HK/(1+RS) )/RS
1784      +              - EXP( -(BS/XS+HK)/2 )*( 1 + C*XS*( 1 + D*XS ) ) )
1785                XS = AS*(-X(I,NG)+1)**2/4
1786                RS = SQRT( 1 - XS )
1787                BVN = BVN + A*W(I,NG)*EXP( -(BS/XS + HK)/2 )
1788      +                    *( EXP( -HK*(1-RS)/(2*(1+RS)) )/RS
1789      +                       - ( 1 + C*XS*( 1 + D*XS ) ) )
1790  20         CONTINUE
1791             BVN = -BVN/TWOPI
1792          ENDIF
1793          IF ( R .GT. 0 ) BVN =  BVN + PHI( -MAX( H, K ) )
1794          IF ( R .LT. 0 ) BVN = -BVN + MAX( ZERO, PHI(-H) - PHI(-K) )
1795       ENDIF
1796       BVNU = BVN
1797       END
1798 *
1799       SUBROUTINE ADAPT(NDIM, MINCLS, MAXCLS, FUNCTN,
1800      &     ABSREQ, RELREQ, LENWRK, WORK, ABSEST, FINEST, INFORM)
1801 *
1802 *   Adaptive Multidimensional Integration Subroutine
1803 *
1804 *   Author: Alan Genz
1805 *           Department of Mathematics
1806 *           Washington State University
1807 *           Pullman, WA 99164-3113 USA
1808 *
1809 *  This subroutine computes an approximation to the integral
1810 *
1811 *      1 1     1
1812 *     I I ... I       FUNCTN(NDIM,X)  dx(NDIM)...dx(2)dx(1)
1813 *      0 0     0
1814 *
1815 ***************  Parameters for ADAPT  ********************************
1816 *
1817 ****** Input Parameters
1818 *
1819 *  NDIM    Integer number of integration variables.
1820 *  MINCLS  Integer minimum number of FUNCTN calls to be allowed; MINCLS
1821 *          must not exceed MAXCLS. If MINCLS < 0, then ADAPT assumes
1822 *          that a previous call of ADAPT has been made with the same
1823 *          integrand and continues that calculation.
1824 *  MAXCLS  Integer maximum number of FUNCTN calls to be used; MAXCLS
1825 *          must be >= RULCLS, the number of function calls required for
1826 *          one application of the basic integration rule.
1827 *           IF ( NDIM .EQ. 1 ) THEN
1828 *              RULCLS = 11
1829 *           ELSE IF ( NDIM .LT. 15 ) THEN
1830 *              RULCLS = 2**NDIM + 2*NDIM*(NDIM+3) + 1
1831 *           ELSE
1832 *              RULCLS = 1 + NDIM*(24-NDIM*(6-NDIM*4))/3
1833 *           ENDIF
1834 *  FUNCTN  Externally declared real user defined integrand. Its
1835 *          parameters must be (NDIM, Z), where Z is a real array of
1836 *          length NDIM.
1837 *  ABSREQ  Real required absolute accuracy.
1838 *  RELREQ  Real required relative accuracy.
1839 *  LENWRK  Integer length of real array WORK (working storage); ADAPT
1840 *          needs LENWRK >= 16*NDIM + 27. For maximum efficiency LENWRK
1841 *          should be about 2*NDIM*MAXCLS/RULCLS if MAXCLS FUNCTN
1842 *          calls are needed. If LENWRK is significantly less than this,
1843 *          ADAPT may be less efficient.
1844 *
1845 ****** Output Parameters
1846 *
1847 *  MINCLS  Actual number of FUNCTN calls used by ADAPT.
1848 *  WORK    Real array (length LENWRK) of working storage. This contains
1849 *          information that is needed for additional calls of ADAPT
1850 *          using the same integrand (input MINCLS < 0).
1851 *  ABSEST  Real estimated absolute accuracy.
1852 *  FINEST  Real estimated value of integral.
1853 *  INFORM  INFORM = 0 for normal exit, when ABSEST <= ABSREQ or
1854 *                     ABSEST <= |FINEST|*RELREQ with MINCLS <= MAXCLS.
1855 *          INFORM = 1 if MAXCLS was too small for ADAPT to obtain the
1856 *                     result FINEST to within the requested accuracy.
1857 *          INFORM = 2 if MINCLS > MAXCLS, LENWRK < 16*NDIM + 27 or
1858 *                     RULCLS > MAXCLS.
1859 *
1860 ************************************************************************
1861 *
1862 *     Begin driver routine. This routine partitions the working storage
1863 *      array and then calls the main subroutine ADBASE.
1864 *
1865       EXTERNAL FUNCTN
1866       INTEGER NDIM, MINCLS, MAXCLS, LENWRK, INFORM
1867       DOUBLE PRECISION
1868      &     FUNCTN, ABSREQ, RELREQ, WORK(LENWRK), ABSEST, FINEST
1869       INTEGER SBRGNS, MXRGNS, RULCLS, LENRUL,
1870      & INERRS, INVALS, INPTRS, INLWRS, INUPRS, INMSHS, INPNTS, INWGTS,
1871      & INLOWR, INUPPR, INWDTH, INMESH, INWORK
1872       IF ( NDIM .EQ. 1 ) THEN
1873          LENRUL = 5
1874          RULCLS = 9
1875       ELSE IF ( NDIM .LT. 12 ) THEN
1876          LENRUL = 6
1877          RULCLS = 2**NDIM + 2*NDIM*(NDIM+2) + 1
1878       ELSE
1879          LENRUL = 6
1880          RULCLS = 1 + 2*NDIM*(1+2*NDIM)
1881       ENDIF
1882       IF ( LENWRK .GE. LENRUL*(NDIM+4) + 10*NDIM + 3 .AND.
1883      &     RULCLS. LE. MAXCLS .AND. MINCLS .LE. MAXCLS ) THEN
1884         MXRGNS = ( LENWRK - LENRUL*(NDIM+4) - 7*NDIM )/( 3*NDIM + 3 )
1885         INERRS = 1
1886         INVALS = INERRS + MXRGNS
1887         INPTRS = INVALS + MXRGNS
1888         INLWRS = INPTRS + MXRGNS
1889         INUPRS = INLWRS + MXRGNS*NDIM
1890         INMSHS = INUPRS + MXRGNS*NDIM
1891         INWGTS = INMSHS + MXRGNS*NDIM
1892         INPNTS = INWGTS + LENRUL*4
1893         INLOWR = INPNTS + LENRUL*NDIM
1894         INUPPR = INLOWR + NDIM
1895         INWDTH = INUPPR + NDIM
1896         INMESH = INWDTH + NDIM
1897         INWORK = INMESH + NDIM
1898         IF ( MINCLS .LT. 0 ) SBRGNS = WORK(LENWRK)
1899         CALL ADBASE(NDIM, MINCLS, MAXCLS, FUNCTN, ABSREQ, RELREQ,
1900      &       ABSEST, FINEST, SBRGNS, MXRGNS, RULCLS, LENRUL,
1901      &       WORK(INERRS), WORK(INVALS), WORK(INPTRS), WORK(INLWRS),
1902      &       WORK(INUPRS), WORK(INMSHS), WORK(INWGTS), WORK(INPNTS),
1903      &       WORK(INLOWR), WORK(INUPPR), WORK(INWDTH), WORK(INMESH),
1904      &       WORK(INWORK), INFORM)
1905         WORK(LENWRK) = SBRGNS
1906        ELSE
1907         INFORM = 2
1908         MINCLS = RULCLS
1909       ENDIF
1910       END
1911       SUBROUTINE BSINIT(NDIM, W, LENRUL, G)
1912 *
1913 *     For initializing basic rule weights and symmetric sum parameters.
1914 *
1915       INTEGER NDIM, LENRUL, RULPTS(6), I, J, NUMNUL, SDIM
1916       PARAMETER ( NUMNUL = 4, SDIM = 12 )
1917       DOUBLE PRECISION W(LENRUL,4), G(NDIM,LENRUL)
1918       DOUBLE PRECISION LAM1, LAM2, LAM3, LAMP, RULCON
1919 *
1920 *     The following code determines rule parameters and weights for a
1921 *      degree 7 rule (W(1,1),...,W(5,1)), two degree 5 comparison rules
1922 *      (W(1,2),...,W(5,2) and W(1,3),...,W(5,3)) and a degree 3
1923 *      comparison rule (W(1,4),...W(5,4)).
1924 *
1925 *       If NDIM = 1, then LENRUL = 5 and total points = 9.
1926 *       If NDIM < SDIM, then LENRUL = 6 and
1927 *                      total points = 1+2*NDIM*(NDIM+2)+2**NDIM.
1928 *       If NDIM > = SDIM, then LENRUL = 6 and
1929 *                      total points = 1+2*NDIM*(1+2*NDIM).
1930 *
1931       DO I = 1,LENRUL
1932          DO J = 1,NDIM
1933             G(J,I) = 0
1934          END DO
1935          DO J = 1,NUMNUL
1936             W(I,J) = 0
1937          END DO
1938       END DO
1939       RULPTS(5) = 2*NDIM*(NDIM-1)
1940       RULPTS(4) = 2*NDIM
1941       RULPTS(3) = 2*NDIM
1942       RULPTS(2) = 2*NDIM
1943       RULPTS(1) = 1
1944       LAMP = 0.85
1945       LAM3 = 0.4707
1946       LAM2 = 4/(15 - 5/LAM3)
1947       W(5,1) = ( 3 - 5*LAM3 )/( 180*(LAM2-LAM3)*LAM2**2 )
1948       IF ( NDIM .LT. SDIM ) THEN
1949          LAM1 = 8*LAM3*(31*LAM3-15)/( (3*LAM3-1)*(5*LAM3-3)*35 )
1950          W(LENRUL,1) = 1/(3*LAM3)**3/2**NDIM
1951       ELSE
1952          LAM1 = ( LAM3*(15 - 21*LAM2) + 35*(NDIM-1)*(LAM2-LAM3)/9 )
1953      &       /  ( LAM3*(21 - 35*LAM2) + 35*(NDIM-1)*(LAM2/LAM3-1)/9 )
1954          W(6,1) = 1/(4*(3*LAM3)**3)
1955       ENDIF
1956       W(3,1) = ( 15 - 21*(LAM3+LAM1) + 35*LAM3*LAM1 )
1957      &     /( 210*LAM2*(LAM2-LAM3)*(LAM2-LAM1) ) - 2*(NDIM-1)*W(5,1)
1958       W(2,1) = ( 15 - 21*(LAM3+LAM2) + 35*LAM3*LAM2 )
1959      &     /( 210*LAM1*(LAM1-LAM3)*(LAM1-LAM2) )
1960       IF ( NDIM .LT. SDIM ) THEN
1961          RULPTS(LENRUL) = 2**NDIM
1962          LAM3 = SQRT(LAM3)
1963          DO I = 1,NDIM
1964             G(I,LENRUL) = LAM3
1965          END DO
1966       ELSE
1967          W(6,1) = 1/(4*(3*LAM3)**3)
1968          RULPTS(6) = 2*NDIM*(NDIM-1)
1969          LAM3 = SQRT(LAM3)
1970          DO I = 1,2
1971             G(I,6) = LAM3
1972          END DO
1973       ENDIF
1974       IF ( NDIM .GT. 1 ) THEN
1975          W(5,2) = 1/(6*LAM2)**2
1976          W(5,3) = 1/(6*LAM2)**2
1977       ENDIF
1978       W(3,2) = ( 3 - 5*LAM1 )/( 30*LAM2*(LAM2-LAM1) )
1979      &     - 2*(NDIM-1)*W(5,2)
1980       W(2,2) = ( 3 - 5*LAM2 )/( 30*LAM1*(LAM1-LAM2) )
1981       W(4,3) = ( 3 - 5*LAM2 )/( 30*LAMP*(LAMP-LAM2) )
1982       W(3,3) = ( 3 - 5*LAMP )/( 30*LAM2*(LAM2-LAMP) )
1983      &     - 2*(NDIM-1)*W(5,3)
1984       W(2,4) = 1/(6*LAM1)
1985       LAMP = SQRT(LAMP)
1986       LAM2 = SQRT(LAM2)
1987       LAM1 = SQRT(LAM1)
1988       G(1,2) = LAM1
1989       G(1,3) = LAM2
1990       G(1,4) = LAMP
1991       IF ( NDIM .GT. 1 ) THEN
1992          G(1,5) = LAM2
1993          G(2,5) = LAM2
1994       ENDIF
1995       DO J = 1, NUMNUL
1996          W(1,J) = 1
1997          DO I = 2,LENRUL
1998             W(1,J) = W(1,J) - RULPTS(I)*W(I,J)
1999          END DO
2000       END DO
2001       RULCON = 2
2002       CALL RULNRM( LENRUL, NUMNUL, RULPTS, W, RULCON )
2003       END
2004       SUBROUTINE RULNRM( LENRUL, NUMNUL, RULPTS, W, RULCON )
2005       INTEGER LENRUL, NUMNUL, I, J, K, RULPTS(*)
2006       DOUBLE PRECISION ALPHA, NORMCF, NORMNL, W(LENRUL, *), RULCON
2007 *
2008 *     Compute orthonormalized null rules.
2009 *
2010       NORMCF = 0
2011       DO I = 1,LENRUL
2012          NORMCF = NORMCF + RULPTS(I)*W(I,1)*W(I,1)
2013       END DO
2014       DO K = 2,NUMNUL
2015          DO I = 1,LENRUL
2016             W(I,K) = W(I,K) - W(I,1)
2017          END DO
2018          DO J = 2,K-1
2019             ALPHA = 0
2020             DO I = 1,LENRUL
2021                ALPHA = ALPHA + RULPTS(I)*W(I,J)*W(I,K)
2022             END DO
2023             ALPHA = -ALPHA/NORMCF
2024             DO I = 1,LENRUL
2025                W(I,K) = W(I,K) + ALPHA*W(I,J)
2026             END DO
2027          END DO
2028          NORMNL = 0
2029          DO I = 1,LENRUL
2030             NORMNL = NORMNL + RULPTS(I)*W(I,K)*W(I,K)
2031          END DO
2032          ALPHA = SQRT(NORMCF/NORMNL)
2033          DO I = 1,LENRUL
2034             W(I,K) = ALPHA*W(I,K)
2035          END DO
2036       END DO
2037       DO J = 2, NUMNUL
2038          DO I = 1,LENRUL
2039             W(I,J) = W(I,J)/RULCON
2040          END DO
2041       END DO
2042       END
2043       SUBROUTINE ADBASE(NDIM, MINCLS, MAXCLS, FUNCTN, ABSREQ, RELREQ,
2044      &     ABSEST, FINEST, SBRGNS, MXRGNS, RULCLS, LENRUL,
2045      &     ERRORS, VALUES, PONTRS, LOWERS,
2046      &     UPPERS, MESHES, WEGHTS, POINTS,
2047      &     LOWER, UPPER, WIDTH, MESH, WORK, INFORM)
2048 *
2049 *        Main adaptive integration subroutine
2050 *
2051       EXTERNAL FUNCTN
2052       INTEGER I, J, NDIM, MINCLS, MAXCLS, SBRGNS, MXRGNS,
2053      &     RULCLS, LENRUL, INFORM, NWRGNS
2054       DOUBLE PRECISION FUNCTN, ABSREQ, RELREQ, ABSEST, FINEST,
2055      &     ERRORS(*), VALUES(*), PONTRS(*),
2056      &     LOWERS(NDIM,*), UPPERS(NDIM,*),
2057      &     MESHES(NDIM,*),WEGHTS(*), POINTS(*),
2058      &     LOWER(*), UPPER(*), WIDTH(*), MESH(*), WORK(*)
2059       INTEGER DIVAXN, TOP, RGNCLS, FUNCLS, DIFCLS
2060  
2061 *
2062 *     Initialization of subroutine
2063 *
2064       INFORM = 2
2065       FUNCLS = 0
2066       CALL BSINIT(NDIM, WEGHTS, LENRUL, POINTS)
2067       IF ( MINCLS .GE. 0) THEN
2068 *
2069 *       When MINCLS >= 0 determine initial subdivision of the
2070 *       integration region and apply basic rule to each subregion.
2071 *
2072          SBRGNS = 0
2073          DO I = 1,NDIM
2074             LOWER(I) = 0
2075             MESH(I) = 1
2076             WIDTH(I) = 1/(2*MESH(I))
2077             UPPER(I) = 1
2078          END DO
2079          DIVAXN = 0
2080          RGNCLS = RULCLS
2081          NWRGNS = 1
2082  10      CALL DIFFER(NDIM, LOWER, UPPER, WIDTH, WORK, WORK(NDIM+1),
2083      &        FUNCTN, DIVAXN, DIFCLS)
2084          FUNCLS = FUNCLS + DIFCLS
2085          IF ( FUNCLS + RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
2086      &        .LE. MINCLS ) THEN
2087             RGNCLS = RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
2088             NWRGNS = NWRGNS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
2089             MESH(DIVAXN) = MESH(DIVAXN) + 1
2090             WIDTH(DIVAXN) = 1/( 2*MESH(DIVAXN) )
2091             GO TO 10
2092          ENDIF
2093          IF ( NWRGNS .LE. MXRGNS ) THEN
2094             DO I = 1,NDIM
2095                UPPER(I) = LOWER(I) + 2*WIDTH(I)
2096                MESH(I) = 1
2097             END DO
2098          ENDIF
2099 *
2100 *     Apply basic rule to subregions and store results in heap.
2101 *
2102  20      SBRGNS = SBRGNS + 1
2103          CALL BASRUL(NDIM, LOWER, UPPER, WIDTH, FUNCTN,
2104      &        WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
2105      &        ERRORS(SBRGNS),VALUES(SBRGNS))
2106          CALL TRESTR(SBRGNS, SBRGNS, PONTRS, ERRORS)
2107          DO I = 1,NDIM
2108             LOWERS(I,SBRGNS) = LOWER(I)
2109             UPPERS(I,SBRGNS) = UPPER(I)
2110             MESHES(I,SBRGNS) = MESH(I)
2111          END DO
2112          DO I = 1,NDIM
2113             LOWER(I) = UPPER(I)
2114             UPPER(I) = LOWER(I) + 2*WIDTH(I)
2115             IF ( LOWER(I)+WIDTH(I) .LT. 1 )  GO TO 20
2116             LOWER(I) = 0
2117             UPPER(I) = LOWER(I) + 2*WIDTH(I)
2118          END DO
2119          FUNCLS = FUNCLS + SBRGNS*RULCLS
2120       ENDIF
2121 *
2122 *     Check for termination
2123 *
2124  30   FINEST = 0
2125       ABSEST = 0
2126       DO I = 1, SBRGNS
2127          FINEST = FINEST + VALUES(I)
2128          ABSEST = ABSEST + ERRORS(I)
2129       END DO
2130       IF ( ABSEST .GT. MAX( ABSREQ, RELREQ*ABS(FINEST) )
2131      &     .OR. FUNCLS .LT. MINCLS ) THEN
2132 *
2133 *     Prepare to apply basic rule in (parts of) subregion with
2134 *     largest error.
2135 *
2136          TOP = PONTRS(1)
2137          RGNCLS = RULCLS
2138          DO I = 1,NDIM
2139             LOWER(I) = LOWERS(I,TOP)
2140             UPPER(I) = UPPERS(I,TOP)
2141             MESH(I) = MESHES(I,TOP)
2142             WIDTH(I) = (UPPER(I)-LOWER(I))/(2*MESH(I))
2143             RGNCLS = RGNCLS*MESH(I)
2144          END DO
2145          CALL DIFFER(NDIM, LOWER, UPPER, WIDTH, WORK, WORK(NDIM+1),
2146      &        FUNCTN, DIVAXN, DIFCLS)
2147          FUNCLS = FUNCLS + DIFCLS
2148          RGNCLS = RGNCLS*(MESH(DIVAXN)+1)/MESH(DIVAXN)
2149          IF ( FUNCLS + RGNCLS .LE. MAXCLS ) THEN
2150             IF ( SBRGNS + 1 .LE. MXRGNS ) THEN
2151 *
2152 *     Prepare to subdivide into two pieces.
2153 *
2154                NWRGNS = 1
2155                WIDTH(DIVAXN) = WIDTH(DIVAXN)/2
2156             ELSE
2157                NWRGNS = 0
2158                WIDTH(DIVAXN) = WIDTH(DIVAXN)
2159      &                        *MESH(DIVAXN)/( MESH(DIVAXN) + 1 )
2160                MESHES(DIVAXN,TOP) = MESH(DIVAXN) + 1
2161             ENDIF
2162             IF ( NWRGNS .GT. 0 ) THEN
2163 *
2164 *     Only allow local subdivision when space is available.
2165 *
2166                DO J = SBRGNS+1,SBRGNS+NWRGNS
2167                   DO I = 1,NDIM
2168                      LOWERS(I,J) = LOWER(I)
2169                      UPPERS(I,J) = UPPER(I)
2170                      MESHES(I,J) = MESH(I)
2171                   END DO
2172                END DO
2173                UPPERS(DIVAXN,TOP) = LOWER(DIVAXN) + 2*WIDTH(DIVAXN)
2174                LOWERS(DIVAXN,SBRGNS+1) = UPPERS(DIVAXN,TOP)
2175             ENDIF
2176             FUNCLS = FUNCLS + RGNCLS
2177             CALL BASRUL(NDIM, LOWERS(1,TOP), UPPERS(1,TOP), WIDTH,
2178      &           FUNCTN, WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
2179      &           ERRORS(TOP), VALUES(TOP))
2180             CALL TRESTR(TOP, SBRGNS, PONTRS, ERRORS)
2181             DO I = SBRGNS+1, SBRGNS+NWRGNS
2182 *
2183 *     Apply basic rule and store results in heap.
2184 *
2185                CALL BASRUL(NDIM, LOWERS(1,I), UPPERS(1,I), WIDTH,
2186      &              FUNCTN, WEGHTS, LENRUL, POINTS, WORK, WORK(NDIM+1),
2187      &              ERRORS(I), VALUES(I))
2188                CALL TRESTR(I, I, PONTRS, ERRORS)
2189             END DO
2190             SBRGNS = SBRGNS + NWRGNS
2191             GO TO 30
2192          ELSE
2193             INFORM = 1
2194          ENDIF
2195       ELSE
2196          INFORM = 0
2197       ENDIF
2198       MINCLS = FUNCLS
2199       END
2200       SUBROUTINE BASRUL( NDIM, A, B, WIDTH, FUNCTN, W, LENRUL, G,
2201      &     CENTER, Z, RGNERT, BASEST )
2202 *
2203 *     For application of basic integration rule
2204 *
2205       EXTERNAL FUNCTN
2206       INTEGER I, LENRUL, NDIM
2207       DOUBLE PRECISION
2208      &     A(NDIM), B(NDIM), WIDTH(NDIM), FUNCTN, W(LENRUL,4),
2209      &     G(NDIM,LENRUL), CENTER(NDIM), Z(NDIM), RGNERT, BASEST
2210       DOUBLE PRECISION
2211      &     FULSUM, FSYMSM, RGNCMP, RGNVAL, RGNVOL, RGNCPT, RGNERR
2212 *
2213 *     Compute Volume and Center of Subregion
2214 *
2215       RGNVOL = 1
2216       DO I = 1,NDIM
2217          RGNVOL = 2*RGNVOL*WIDTH(I)
2218          CENTER(I) = A(I) + WIDTH(I)
2219       END DO
2220       BASEST = 0
2221       RGNERT = 0
2222 *
2223 *     Compute basic rule and error
2224 *
2225  10   RGNVAL = 0
2226       RGNERR = 0
2227       RGNCMP = 0
2228       RGNCPT = 0
2229       DO I = 1,LENRUL
2230          FSYMSM = FULSUM(NDIM, CENTER, WIDTH, Z, G(1,I), FUNCTN)
2231 *     Basic Rule
2232          RGNVAL = RGNVAL + W(I,1)*FSYMSM
2233 *     First comparison rule
2234          RGNERR = RGNERR + W(I,2)*FSYMSM
2235 *     Second comparison rule
2236          RGNCMP = RGNCMP + W(I,3)*FSYMSM
2237 *     Third Comparison rule
2238          RGNCPT = RGNCPT + W(I,4)*FSYMSM
2239       END DO
2240 *
2241 *     Error estimation
2242 *
2243       RGNERR = SQRT(RGNCMP**2 + RGNERR**2)
2244       RGNCMP = SQRT(RGNCPT**2 + RGNCMP**2)
2245       IF ( 4*RGNERR .LT. RGNCMP ) RGNERR = RGNERR/2
2246       IF ( 2*RGNERR .GT. RGNCMP ) RGNERR = MAX( RGNERR, RGNCMP )
2247       RGNERT = RGNERT +  RGNVOL*RGNERR
2248       BASEST = BASEST +  RGNVOL*RGNVAL
2249 *
2250 *     When subregion has more than one piece, determine next piece and
2251 *      loop back to apply basic rule.
2252 *
2253       DO I = 1,NDIM
2254          CENTER(I) = CENTER(I) + 2*WIDTH(I)
2255          IF ( CENTER(I) .LT. B(I) ) GO TO 10
2256          CENTER(I) = A(I) + WIDTH(I)
2257       END DO
2258       END
2259       DOUBLE PRECISION FUNCTION FULSUM(S, CENTER, HWIDTH, X, G, F)
2260 *
2261 ****  To compute fully symmetric basic rule sum
2262 *
2263       EXTERNAL F
2264       INTEGER S, IXCHNG, LXCHNG, I, L
2265       DOUBLE PRECISION CENTER(S), HWIDTH(S), X(S), G(S), F
2266       DOUBLE PRECISION INTSUM, GL, GI
2267       FULSUM = 0
2268 *
2269 *     Compute centrally symmetric sum for permutation of G
2270 *
2271  10   INTSUM = 0
2272       DO I = 1,S
2273          X(I) = CENTER(I) + G(I)*HWIDTH(I)
2274       END DO
2275  20   INTSUM = INTSUM + F(S,X)
2276       DO I = 1,S
2277          G(I) = -G(I)
2278          X(I) = CENTER(I) + G(I)*HWIDTH(I)
2279          IF ( G(I) .LT. 0 ) GO TO 20
2280       END DO
2281       FULSUM = FULSUM + INTSUM
2282 *
2283 *     Find next distinct permuation of G and loop back for next sum
2284 *
2285       DO I = 2,S
2286          IF ( G(I-1) .GT. G(I) ) THEN
2287             GI = G(I)
2288             IXCHNG = I - 1
2289             DO L = 1,(I-1)/2
2290                GL = G(L)
2291                G(L) = G(I-L)
2292                G(I-L) = GL
2293                IF (  GL  .LE. GI ) IXCHNG = IXCHNG - 1
2294                IF ( G(L) .GT. GI ) LXCHNG = L
2295             END DO
2296             IF ( G(IXCHNG) .LE. GI ) IXCHNG = LXCHNG
2297             G(I) = G(IXCHNG)
2298             G(IXCHNG) = GI
2299             GO TO 10
2300          ENDIF
2301       END DO
2302 *
2303 *     End loop for permutations of G and associated sums
2304 *
2305 *     Restore original order to G's
2306 *
2307       DO I = 1,S/2
2308          GI = G(I)
2309          G(I) = G(S+1-I)
2310          G(S+1-I) = GI
2311       END DO
2312       END
2313       SUBROUTINE DIFFER(NDIM, A, B, WIDTH, Z, DIF, FUNCTN,
2314      &     DIVAXN, DIFCLS)
2315 *
2316 *     Compute fourth differences and subdivision axes
2317 *
2318       EXTERNAL FUNCTN
2319       INTEGER I, NDIM, DIVAXN, DIFCLS
2320       DOUBLE PRECISION
2321      &     A(NDIM), B(NDIM), WIDTH(NDIM), Z(NDIM), DIF(NDIM), FUNCTN
2322       DOUBLE PRECISION FRTHDF, FUNCEN, WIDTHI
2323       DIFCLS = 0
2324       DIVAXN = MOD( DIVAXN, NDIM ) + 1
2325       IF ( NDIM .GT. 1 ) THEN
2326          DO I = 1,NDIM
2327             DIF(I) = 0
2328             Z(I) = A(I) + WIDTH(I)
2329          END DO
2330  10      FUNCEN = FUNCTN(NDIM, Z)
2331          DO I = 1,NDIM
2332             WIDTHI = WIDTH(I)/5
2333             FRTHDF = 6*FUNCEN
2334             Z(I) = Z(I) - 4*WIDTHI
2335             FRTHDF = FRTHDF + FUNCTN(NDIM,Z)
2336             Z(I) = Z(I) + 2*WIDTHI
2337             FRTHDF = FRTHDF - 4*FUNCTN(NDIM,Z)
2338             Z(I) = Z(I) + 4*WIDTHI
2339             FRTHDF = FRTHDF - 4*FUNCTN(NDIM,Z)
2340             Z(I) = Z(I) + 2*WIDTHI
2341             FRTHDF = FRTHDF + FUNCTN(NDIM,Z)
2342 *     Do not include differences below roundoff
2343             IF ( FUNCEN + FRTHDF/8 .NE. FUNCEN )
2344      &           DIF(I) = DIF(I) + ABS(FRTHDF)*WIDTH(I)
2345             Z(I) = Z(I) - 4*WIDTHI
2346          END DO
2347          DIFCLS = DIFCLS + 4*NDIM + 1
2348          DO I = 1,NDIM
2349             Z(I) = Z(I) + 2*WIDTH(I)
2350             IF ( Z(I) .LT. B(I) ) GO TO 10
2351             Z(I) = A(I) + WIDTH(I)
2352          END DO
2353          DO I = 1,NDIM
2354             IF ( DIF(DIVAXN) .LT. DIF(I) ) DIVAXN = I
2355          END DO
2356       ENDIF
2357       END
2358       SUBROUTINE TRESTR(POINTR, SBRGNS, PONTRS, RGNERS)
2359 ****BEGIN PROLOGUE TRESTR
2360 ****PURPOSE TRESTR maintains a heap for subregions.
2361 ****DESCRIPTION TRESTR maintains a heap for subregions.
2362 *            The subregions are ordered according to the size of the
2363 *            greatest error estimates of each subregion (RGNERS).
2364 *
2365 *   PARAMETERS
2366 *
2367 *     POINTR Integer.
2368 *            The index for the subregion to be inserted in the heap.
2369 *     SBRGNS Integer.
2370 *            Number of subregions in the heap.
2371 *     PONTRS Real array of dimension SBRGNS.
2372 *            Used to store the indices for the greatest estimated errors
2373 *            for each subregion.
2374 *     RGNERS Real array of dimension SBRGNS.
2375 *            Used to store the greatest estimated errors for each
2376 *            subregion.
2377 *
2378 ****ROUTINES CALLED NONE
2379 ****END PROLOGUE TRESTR
2380 *
2381 *   Global variables.
2382 *
2383       INTEGER POINTR, SBRGNS
2384       DOUBLE PRECISION PONTRS(*), RGNERS(*)
2385 *
2386 *   Local variables.
2387 *
2388 *   RGNERR Intermediate storage for the greatest error of a subregion.
2389 *   SUBRGN Position of child/parent subregion in the heap.
2390 *   SUBTMP Position of parent/child subregion in the heap.
2391 *
2392       INTEGER SUBRGN, SUBTMP
2393       DOUBLE PRECISION RGNERR
2394 *
2395 ****FIRST PROCESSING STATEMENT TRESTR
2396 *
2397       RGNERR = RGNERS(POINTR)
2398       IF ( POINTR .EQ. PONTRS(1)) THEN
2399 *
2400 *        Move the new subregion inserted at the top of the heap
2401 *        to its correct position in the heap.
2402 *
2403          SUBRGN = 1
2404  10      SUBTMP = 2*SUBRGN
2405          IF ( SUBTMP .LE. SBRGNS ) THEN
2406             IF ( SUBTMP .NE. SBRGNS ) THEN
2407 *
2408 *              Find maximum of left and right child.
2409 *
2410           IF ( RGNERS(INT(PONTRS(SUBTMP))) .LT.
2411      $         RGNERS(INT(PONTRS(SUBTMP+1))) ) SUBTMP = SUBTMP + 1
2412             ENDIF
2413 *
2414 *           Compare maximum child with parent.
2415 *           If parent is maximum, then done.
2416 *
2417             IF ( RGNERR .LT. RGNERS(INT(PONTRS(SUBTMP))) ) THEN
2418 *
2419 *              Move the pointer at position subtmp up the heap.
2420 *
2421                PONTRS(SUBRGN) = INT(PONTRS(SUBTMP))
2422                SUBRGN = SUBTMP
2423                GO TO 10
2424             ENDIF
2425          ENDIF
2426       ELSE
2427 *
2428 *        Insert new subregion in the heap.
2429 *
2430          SUBRGN = SBRGNS
2431  20      SUBTMP = SUBRGN/2
2432          IF ( SUBTMP .GE. 1 ) THEN
2433 *
2434 *           Compare child with parent. If parent is maximum, then done.
2435 *
2436             IF ( RGNERR .GT. RGNERS(INT(PONTRS(SUBTMP))) ) THEN
2437 *
2438 *              Move the pointer at position subtmp down the heap.
2439 *
2440                PONTRS(SUBRGN) = PONTRS(SUBTMP)
2441                SUBRGN = SUBTMP
2442                GO TO 20
2443             ENDIF
2444          ENDIF
2445       ENDIF
2446       PONTRS(SUBRGN) = POINTR
2447 *
2448 ****END TRESTR
2449 *
2450       END