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