MODULE Mult_Normal
! Code converted using TO_F90 by Alan Miller
! Date: 2003-12-15  Time: 09:34:54

!     ALGORITHM AS285 APPL. STATIST. (1993) VOL.42, NO.3

!  FINDS THE PROBABILITY THAT A NORMALLY DISTRIBUTED RANDOM
!  N-VECTOR WITH MEAN 0 AND COVARIANCE COVAR FALLS
!  IN AREA ENCLOSED BY THE EXTERNAL USER-DEFINED FUNCTION F.

IMPLICIT NONE
INTEGER, SAVE  :: ix, iy, iz


CONTAINS


SUBROUTINE mulnor(covar, n, f, sdev, rhigh, init, iseed, t, work, prob,  &
    iter, ifault)

IMPLICIT NONE
REAL, INTENT(IN OUT)     :: covar(*)
INTEGER, INTENT(IN)      :: n
REAL, INTENT(IN)         :: sdev
REAL, INTENT(IN)         :: rhigh
INTEGER, INTENT(IN)      :: init
INTEGER, INTENT(IN OUT)  :: iseed
REAL, INTENT(IN OUT)     :: t(*)
REAL, INTENT(OUT)        :: work(n,*)
REAL, INTENT(OUT)        :: prob
INTEGER, INTENT(OUT)     :: iter
INTEGER, INTENT(OUT)     :: ifault

INTERFACE
  FUNCTION f(n, v) RESULT(fn_val)
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f
END INTERFACE

INTEGER  :: init1, k, ii, kk, q, nn, nullty, i, j
REAL     :: p, var, cox

REAL, PARAMETER     :: zero = 0.0
INTEGER, PARAMETER  :: maxitr = 10000

ifault = 0
init1 = init
IF (init < 25) init1 = 25
IF (init > 1000) init1 = 1000
cox = 1.0 + 2.0/init1

!  FIND CHOLESKY DECOMPOSITION OF COVARIANCE MATRIX

nn = n*(n+1)/2
CALL chol(covar, n, nn, t, nullty, ifault)
q = n - nullty
ifault = ifault*ifault
IF (ifault /= 0) RETURN

!  TRANSPOSE CHOLESKY FACTOR, OMITTING COLUMNS OF ZEROES.
!  RESULTING VECTOR IS PACKED FORM OF T'J,
!  WRITTEN ROWWISE BEGINNING WITH LAST ROW.

DO  i = 1,n
  work(1:i,i) = zero
END DO
ii = 0
DO  i = 1,n
  DO  j = 1,i
    ii = ii + 1
    work(i,j) = t(ii)
  END DO
END DO
j = 0
DO  i = n,1,-1
  IF (j == nullty) EXIT
  IF (work(i,i) == zero) THEN
    j = j+1
    DO  k = i+1,n
      DO  kk = i, n-j
        work(k,kk) = work(k,kk+1)
      END DO
    END DO
  END IF
END DO

k = 0
DO  j = n,1,-1
  DO  i = 1, MIN(j,q)
    k = k + 1
    t(k) = work(j,i)
  END DO
END DO

!  TAKE PILOT SAMPLE; ESTIMATE VARIANCE OF PROB FROM ORTHP.

prob = zero
var = zero
DO  k=1,init1
  CALL orthp(f, rhigh, n, q, t, work, p, iseed, ifault)
  IF (ifault /= 0) RETURN
  prob = prob + p
  var = var + p*p
END DO
prob = prob/init1
var = var/init1 - prob*prob

iter = INT(cox*var/(sdev*sdev)) + 1
IF (iter > maxitr) ifault = 5

DO  k=init1+1, iter
  CALL orthp(f, rhigh, n, q, t, work, p, iseed, ifault)
  IF (ifault /= 0) RETURN
  prob = prob + (p - prob)/k
END DO

RETURN
END SUBROUTINE mulnor



SUBROUTINE chol(a, n, nn, u, nullty, ifault)

IMPLICIT NONE
INTEGER, INTENT(IN)   :: nn
REAL, INTENT(IN)      :: a(nn)
INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN OUT)  :: u(nn)
INTEGER, INTENT(OUT)  :: nullty
INTEGER, INTENT(OUT)  :: ifault

REAL     :: eta2, x, w
INTEGER  :: j, k, ii, icol, l, kk, m, irow, i

REAL, PARAMETER  :: eta = 1.0E-05, zero = 0.0

! zabs(x) = ABS(x)
! zsqrt(x) = SQRT(x)

ifault = 1
IF (n <= 0) RETURN
ifault = 3
IF (nn /= n*(n+1)/2) RETURN
ifault = 2
nullty = 0
j = 1
k = 0
eta2 = eta*eta
ii = 0
DO  icol = 1,n
  ii = ii + icol
  x = eta2 * a(ii)
  l = 0
  kk = 0
  DO  irow = 1, icol
    kk = kk + irow
    k = k+1
    w = a(k)
    m = j
    DO  i = 1, irow
      l = l + 1
      IF (i == irow) EXIT
      w = w - u(l) * u(m)
      m = m + 1
    END DO
    IF (irow == icol) EXIT
    IF (u(l) == zero) GO TO 30
    u(k) = w / u(l)
    CYCLE
    30 IF (w * w > abs(x * a(kk))) RETURN
    u(k) = zero
  END DO
  IF (abs(w) <= abs(eta * a(k))) GO TO 60
  IF (w < zero) RETURN
  u(k) = sqrt(w)
  GO TO 70
  60 u(k) = zero
  nullty = nullty + 1
  70 j = j + icol
END DO
ifault = 0
RETURN
END SUBROUTINE chol



FUNCTION chiprb(x, df) RESULT(fn_val)

!  FINDS THE PROBABILITY THAT A CHI-SQUARED RANDOM VARIABLE
!  WITH DF DEGREES OF FREEDOM EXCEEDS X.
!  ADAPTED FROM HILL AND PIKE (1967) ALG. 299, CHI-SQUARED INTEGRAL,
!  COMM. ACM, VOL.10, 243-244.

IMPLICIT NONE
REAL, INTENT(IN)     :: x
INTEGER, INTENT(IN)  :: df
REAL                 :: fn_val

REAL     :: a, sqrta, y, c, e, z, s
INTEGER  :: odd, n2, i
!     REAL ERF

REAL, PARAMETER  :: rtpirp = 0.564189583547756
REAL, PARAMETER  :: big = 88.03, small = -85.2, zero = 0.0, one = 1.0, &
   half = 0.5

!  NOTE: IF THE FUNCTION "ERF" IS AVAILABLE, THE ABOVE STATEMENT
!  FUNCTION MAY BE REPLACED BY
!     ZERFP1(A) = ERF(A) + ONE

s = zero
IF (x > big) GO TO 50
s = one
IF (x < small) GO TO 50

a = half*x
y = EXP(-a)
odd = MOD(df,2)
IF (odd == 0) THEN
  s = y
  e = one
  z = zero
ELSE
  sqrta = SQRT(a)
  s = zerfp1(-sqrta)
  e = rtpirp/sqrta
  z = -half
END IF
IF (df < 3) GO TO 50

n2 = df/2 - 1 + odd
c = zero
DO  i = 1,n2
  z = z + 1
  e = e*(a/z)
  c = c + e
END DO
s = s + c*y

50 fn_val = s
RETURN
END FUNCTION chiprb



SUBROUTINE orthp(f, rhigh, n, q, t, z, prob, iseed, ifault)

!  FINDS THE PSEUDO-MONTE CARLO ESTIMATE OF PROB USING AN
!  ORTHONORMAL SYSTEM OF VECTORS.

IMPLICIT NONE
REAL, INTENT(IN)         :: rhigh
INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(IN)      :: q
REAL, INTENT(IN)         :: t(*)
REAL, INTENT(IN OUT)     :: z(n,*)
REAL, INTENT(OUT)        :: prob
INTEGER, INTENT(IN OUT)  :: iseed
INTEGER, INTENT(OUT)     :: ifault

INTERFACE
  FUNCTION f(n, v) RESULT(fn_val)
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f
END INTERFACE

INTEGER  :: i, j, k, jj, ii, qp1, qp2
REAL     :: s, a, b, r, fab

! EXTERNAL f
REAL, PARAMETER  :: zero = 0.0, rt2rcp = 0.707106781

prob = zero
qp1 = q+1
qp2 = q+2

!  PUT RANDOMLY ORIENTED ORTHONORMAL SYSTEM OF Q-VECTORS IN Z.

CALL rnortm(n,q, qp1, 1, 1, z(1,qp1), 1, z(1,qp1), iseed, z,  &
    z(1,qp1), fab, z(1,qp2), z(1,qp2), ifault)

!  REPLACE EACH COLUMN OF Z BY TRANS(CHOLESKY FACTOR)*J*COLUMN OF Z.

DO  i=1,q
  jj = 0
  DO  j = n,1,-1
    s = zero
    DO  k=1, MIN(j,q)
      jj = jj + 1
      s = s + t(jj)*z(k,i)
    END DO
    z(j,i) = s
  END DO
END DO

!  ESTIMATE THE PROBABILITY THAT N(0,COVAR) IS WITHIN THE BOUNDARY GIVEN BY F.

DO  i = 1,q
  a = zero
  b = rhigh
  CALL zerovc (f, r, a, b, n, z(1,i), z(1,qp2), ifault)
  prob = prob - chiprb(r*r, q)
END DO

DO  i = 1,q
  a = zero
  b = rhigh
  DO  k = 1,n
    z(k,qp1) = -z(k,i)
  END DO
  CALL zerovc (f, r, a, b, n, z(1,qp1), z(1,qp2), ifault)
  prob = prob - chiprb(r*r, q)
END DO

DO  i = 1,q-1
  DO  j = i+1,q
    DO  jj= -1,1,2
      DO  ii = -1,1,2
        DO  k=1,n
          z(k,qp1) = rt2rcp*(z(k,i)*ii + z(k,j)*jj)
        END DO
        a = zero
        b = rhigh
        CALL zerovc (f,r,a,b,n,z(1,qp1),z(1,qp2),ifault)
        prob = prob - chiprb(r*r,q)
      END DO
    END DO
  END DO
END DO
prob = prob/(2.0*q*q) + 1.0

RETURN
END SUBROUTINE orthp



SUBROUTINE rnortm(lda, n, np1, ibconf, nb, b, ndbi, dbi, iseed, a,  &
    chisq, fab, u, w, ifault)

!      ALGORITHM AS 127  APPL. STATIST.  (1978) VOL.27, NO.2

!      RNORTM GENERATES ORTHOGONAL MATRICES A FROM A DISTRIBUTION WITH
!      DENSITY FUNCTION DEFINED ON THE GROUP OF ORTHOGONAL MATRICES.
!      WHEN THE INPUT PARAMETER IBCONF = 1 THE MATRICES ARE
!      GENERATED FROM THE INVARIANT HAAR MEASURE.

IMPLICIT NONE
INTEGER, INTENT(IN)      :: lda
INTEGER, INTENT(IN)      :: n
INTEGER, INTENT(IN)      :: np1
INTEGER, INTENT(IN)      :: ibconf
INTEGER, INTENT(IN)      :: nb
REAL, INTENT(IN OUT)     :: b(nb)
INTEGER, INTENT(IN)      :: ndbi
REAL, INTENT(IN OUT)     :: dbi(ndbi)
INTEGER, INTENT(IN OUT)  :: iseed
REAL, INTENT(OUT)        :: a(lda,n)
REAL, INTENT(IN OUT)     :: chisq(n)
REAL, INTENT(OUT)        :: fab
REAL, INTENT(OUT)        :: u(np1)
REAL, INTENT(OUT)        :: w(n)
INTEGER, INTENT(OUT)     :: ifault

REAL    :: ul, wl, wil, wwil, t2, hv, hw
INTEGER :: i, j, l, np1l
INTEGER :: lp1

REAL, PARAMETER  :: zero = 0.0, one = 1.0

!      STATEMENT FUNCTIONS FOR DOUBLE PRECISION
!      DOUBLE PRECISION SQRT, SIGN, ABS
!      SQRT(UL) = DSQRT(UL)
!      SIGN(UL, WL) = DSIGN(UL, WL)
!      ABS(UL) = DABS(UL)

!      SECTION 1 INITIALIZATION

ifault = 0
GO TO 4

!      ERROR RETURNS

2 ifault = 2
RETURN

!      INITIALIZE DENSITY, H(STORED IN A), AND IBSUB

4 fab = one
DO  i=1,n
  a(i,1:n) = zero
  a(i,i) = one
END DO

!      CONSTRUCT  A  ONE COLUMN AT A TIME

DO  l=1,n
  np1l = np1 - l
  
!      SECTION 2
!      CONSTRUCT U(L...N) UNIFORMLY DISTRIBUTED ON THE SURFACE OF THE
!      (N+1-L)-SPHERE OF RADIUS UL = SQRT(CHISQ(L)).
  
  CALL normal(u(l), np1l+1, iseed, chisq(l))
  ul = SQRT(chisq(l))
  IF (ul == zero) GO TO 2
  IF (l == n) GO TO 301
  
!      SECTION 3
!      CALCULATE  W = B*U  AND STORE ITS LENGTH IN WL.
!      SIGN OF WL PREVENTS LOSS OF SIGNIFICANCE IN WWIL.
  
!      B IS IDENTITY
  
  wl = SIGN(ul,-u(l))
  DO  j=l,n
    w(j) = u(j)
  END DO
  
!      SECTION 4
!      PROJECT  W  ONTO ORTHOGOANL COMPLEMENT OF FIRST L-1
!      COLUMNS OF A, NORMALIZE, AND STORE IN L COLUMN OF A.
!      A(L) = H(L) * W * WIL
!      CALCULATE PROJECTION MATRIX H(L) FOR ORTHOGOANL
!      COMPLEMENT OF FIRST L COLUMNS OF A AND STORE IN LAST
!      N-L COLUMNS OF A.
  
  lp1 = l+1
  wil = one/wl
  wwil = one/(wl - w(l))
  DO  i=1,n
    t2 = zero
    DO  j=l,n
      t2 = t2 + a(i,j) * w(j)
    END DO
    hv = t2 * wil
    hw = (hv - a(i,l)) * wwil
    a(i,l) = hv
    DO  j = lp1, n
      a(i,j) = a(i,j) - hw*w(j)
    END DO
  END DO
  
!      END OF L COLUMN OF MATRIX A
  
!      SECTION 5
!      WHEN L = N ONLY SIGN NEEDS TO BE CHOSEN
!      REMEMBER MATRIX H(N) IS STORED IN N COLUMN OF A.
  
  301 IF (u(l) > zero) CYCLE
  DO  i=1,n
    a(i,l) = -a(i,l)
  END DO
  
!      FAB = ABS(B(LAST) * DBI(N)*FAB)
  
END DO
fab = ABS(fab)

RETURN
END SUBROUTINE rnortm



SUBROUTINE zerovc (f, r, r0, r1, n, v, y, ifault)

!  CALCULATES THE SOLUTION R TO F(R*V) = 0.

IMPLICIT NONE
REAL, INTENT(OUT)     :: r
REAL, INTENT(IN OUT)  :: r0
REAL, INTENT(IN OUT)  :: r1
INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN)      :: v(n)
REAL, INTENT(OUT)     :: y(n)
INTEGER, INTENT(OUT)  :: ifault

INTERFACE
  FUNCTION f(n, v) RESULT(fn_val)
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f
END INTERFACE

INTEGER  :: j
REAL     :: f0, f1, fr
REAL, PARAMETER     :: half = 0.5, eta = 1.0E-05
INTEGER, PARAMETER  :: niter = 500

y(1:n) = r0*v(1:n)
f0 = f(n,y)
y(1:n) = r1*v(1:n)
f1 = f(n,y)
IF (SIGN(half,f1) == SIGN(half,f0)) THEN
  ifault = 3
  RETURN
END IF
DO  j = 1,niter
  r = r1 - f1*(r1-r0)/(f1-f0)
  y(1:n) = r*v(1:n)
  fr = f(n, y)
  IF (abs(r-r1) < eta) GO TO 50
  IF (SIGN(half,fr) /= SIGN(half,f1)) THEN
    r0=r1
    f0=f1
  ELSE
    f0 = f0 * half
  END IF
  r1 = r
  f1 = fr
END DO
ifault = 6

50 RETURN
END SUBROUTINE zerovc


FUNCTION zerfp1(a) RESULT(fn_val)
IMPLICIT NONE
REAL, INTENT(IN)  :: a
REAL              :: fn_val

REAL, PARAMETER  :: sqrt2 = 1.41421356237309

fn_val = 2.0*alnorm(sqrt2*a, .false.)
RETURN
END FUNCTION zerfp1



SUBROUTINE normal(u, np1, iseed, chisq)

!  GENERATES AN (NP1-1)-VECTOR OF INDEPENDENT NORMAL VARIATES IN U.
!  ON OUTPUT, CHISQ CONTAINS THE SQUARED LENGTH OF THE VECTOR.

IMPLICIT NONE

INTEGER, INTENT(IN)      :: np1
REAL, INTENT(OUT)        :: u(np1)
INTEGER, INTENT(IN OUT)  :: iseed
REAL, INTENT(OUT)        :: chisq

INTEGER  :: i, ifault

! COMMON /rand/ ix,iy,iz
REAL, PARAMETER  :: zero = 0.0

chisq = zero
DO  i = 1, np1-1
  CALL ppnd(random(), ifault, u(i))
  chisq = chisq + u(i)*u(i)
END DO
RETURN
END SUBROUTINE normal



SUBROUTINE ppnd(p, ifault, fn_val)

!  ALGORITHM AS 111  APPL. STATIST. (1977), VOL.26, NO.1

!  PRODUCES NORMAL DEVIATE CORRESPONDING TO LOWER TAIL AREA OF P
!  REAL VERSION FOR EPS = 2 **(-31)
!  THE HASH SUMS ARE THE SUMS OF THE MODULI OF THE COEFFICIENTS
!  THEY HAVE NO INHERENT MEANINGS BUT ARE INCLUDED FOR USE IN
!  CHECKING TRANSCRIPTIONS
!  STANDARD FUNCTIONS ABS, ALOG AND SQRT ARE USED

IMPLICIT NONE
REAL, INTENT(IN)      :: p
INTEGER, INTENT(OUT)  :: ifault
REAL, INTENT(OUT)     :: fn_val

REAL :: q, r

REAL, PARAMETER  :: zero = 0.0E0, half = 0.5E0, one = 1.0E0
REAL, PARAMETER  :: split = 0.42E0
REAL, PARAMETER  :: a0 = 2.50662823884E0
REAL, PARAMETER  :: a1 = -18.61500062529E0
REAL, PARAMETER  :: a2 = 41.39119773534E0
REAL, PARAMETER  :: a3 = -25.44106049637E0
REAL, PARAMETER  :: b1 = -8.47351093090E0
REAL, PARAMETER  :: b2 = 23.08336743743E0
REAL, PARAMETER  :: b3 = -21.06224101826E0
REAL, PARAMETER  :: b4 = 3.13082909833E0
REAL, PARAMETER  :: c0 = -2.78718931138E0
REAL, PARAMETER  :: c1 = -2.29796479134E0
REAL, PARAMETER  :: c2 = 4.85014127135E0
REAL, PARAMETER  :: c3 = 2.32121276858E0
REAL, PARAMETER  :: d1 = 3.54388924762E0
REAL, PARAMETER  :: d2 = 1.63706781897E0

ifault = 0
q = p - half
IF (ABS(q) > split) GO TO 1
r = q*q
fn_val = q * (((a3*r + a2)*r + a1) * r + a0) /  &
    ((((b4*r + b3)*r + b2) * r + b1) * r + one)
RETURN

1 r = p
IF (q > zero)r = one - p
IF (r <= zero) GO TO 2
r = SQRT(-LOG(r))
fn_val = (((c3 * r + c2) * r + c1) * r + c0)/ ((d2*r + d1) * r + one)
IF (q < zero) fn_val = -fn_val
RETURN

2 ifault = 1
fn_val = zero
RETURN
END SUBROUTINE ppnd



FUNCTION random() RESULT(fn_val)

!  ALGORITHM AS 183  APPL. STATIST. (1982) VOL. 31, NO. 2

!  RETURNS A PSEUDO-RANDOM NUMBER RECTANGULARLY DISTRIBUTED BETWEEN 0 AND 1.

!  IX, IY, AND IZ SHOULD BE SET TO INTEGER VALUES BETWEEN
!  1 AND 30000 BEFORE FIRST ENTRY

!  INTEGER ARITHMETIC UP TO 30323 IS REQUIRED

REAL  :: fn_val

! INTEGER :: ix, iy, iz
! COMMON /rand/ ix,iy,iz

ix = 171 * MOD(ix, 177) - 2 * (ix / 177)
iy = 172 * MOD(iy, 176) - 35 * (iy / 176)
iz = 170 * MOD(iz, 178) - 63 * (iz / 178)

IF (ix < 0) ix = ix + 30269
IF (iy < 0) iy = iy + 30307
IF (iz < 0) iz = iz + 30323

fn_val = MOD(REAL(ix)/30269.0 + REAL(iy)/30307.0 + REAL(iz)/30323.0, 1.0)
RETURN
END FUNCTION random



FUNCTION alnorm(x, upper) RESULT(fn_val)

!  ALGORITHM AS 66  APPL. STATIST. (1973) VOL. 22, NO. 3

!  EVALUATES THE TAIL AREA OF THE STANDARDIZED NORMAL CURVE
!  FROM X TO INFINITY IF UPPER IS .TRUE. OR
!  FROM MINUS INFINITY TO X IF UPPER IS .FALSE.


REAL, INTENT(IN)     :: x
LOGICAL, INTENT(IN)  :: upper
REAL                 :: fn_val

REAL     :: z, y
LOGICAL  :: up

!  LTONE AND UTZERO MUST BE SET TO SUIT THE PARTICULAR COMPUTER
!  (SEE INTRODUCTORY TEXT)

REAL, PARAMETER  :: ltone = 7.0, utzero = 18.66
REAL, PARAMETER  :: zero = 0.0, half = 0.5, one = 1.0, con = 1.28

up = upper
z = x
IF (z >= zero) GO TO 10
up = .NOT. up
z = -z
10    IF (z <= ltone .OR. up .AND. z <= utzero) GO TO 20
fn_val = zero
GO TO 40
20 y = half *z * z
IF (z > con) GO TO 30

fn_val = half - z*(0.398942280444 - 0.399903438504 * y /  &
    (y + 5.75885480458 - 29.8213557808 / (y + 2.62433121679 + 48.6959930692 /  &
    (y + 5.92885724438))))
GO TO 40

30 fn_val = 0.398942280385 * EXP(-y) / (z - 3.8052E-08 + 1.00000615302 /  &
    (z + 3.98064794E-4 + 1.98615381364 / (z - 0.151679116635 + 5.29330324926 /  &
    (z + 4.8385912808 - 15.1508972451 /  &
    (z + 0.742380924027 + 30.789933034 / (z + 3.99019417011))))))

40 IF (.NOT. up) fn_val = one - fn_val
RETURN
END FUNCTION alnorm

END MODULE Mult_Normal



PROGRAM Driver285
!  DRIVING PROGRAM
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-12-15  Time: 12:20:56

USE Mult_Normal
IMPLICIT NONE

REAL     :: covar(1275), work(50,52), t(1275)
INTEGER  :: iseed, n, ifault, j, iter, init, nn, nullty
REAL     :: sdev, prob, rhigh
! COMMON /rand/ ix,iy,iz

REAL, PARAMETER  :: eta = 1.0E-5

INTERFACE
  FUNCTION f1(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f1

  FUNCTION f2(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f2

  FUNCTION f3(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f3

  FUNCTION f4(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f4

  FUNCTION f5(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f5

  FUNCTION f6(n,v) RESULT(fn_val)
    IMPLICIT NONE
    INTEGER, INTENT(IN)  :: n
    REAL, INTENT(IN)     :: v(:)
    REAL                 :: fn_val
  END FUNCTION f6
END INTERFACE

! READ IN COVARIANCE AND INITIALIZE RANDOM NORMAL GENERATOR.

OPEN(1, file='FORT.1', status='old')
READ (1, *) n, sdev, rhigh, init, ix, iy, iz
READ (1, *) (covar(j),j=1,n*(n+1)/2)

!  FIND SMALLEST NONZERO EIGENVALUE OF COVARIANCE MATRIX

nn = n*(n+1)/2
CALL chol(covar, n, nn, t, nullty, ifault)
CALL tdiag(n, eta, covar, work(1,1), work(1,2), work(1,3), ifault)
CALL lrvt(n, eta, work(1,1), work(1,2), work(1,3), ifault)
rhigh = rhigh/SQRT(work(nullty+1,1))

! Run through the supplied functions F1 ... F6

DO n = 2, 5
  CALL mulnor(covar, n, f1, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
  WRITE(*, *) 'Function F1, with n = ', n
  WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
  WRITE(*, *) 'ITER = ',iter
END DO

n = 2
CALL mulnor(covar, n, f2, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
WRITE(*, *)
WRITE(*, *) 'Function F2'
WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
WRITE(*, *) 'ITER = ',iter

n = 4
CALL mulnor(covar, n, f3, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
WRITE(*, *)
WRITE(*, *) 'Function F3'
WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
WRITE(*, *) 'ITER = ',iter

n = 3
CALL mulnor(covar, n, f4, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
WRITE(*, *)
WRITE(*, *) 'Function F4'
WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
WRITE(*, *) 'ITER = ',iter

n = 2
CALL mulnor(covar, n, f5, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
WRITE(*, *)
WRITE(*, *) 'Function F5'
WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
WRITE(*, *) 'ITER = ',iter

n = 2
CALL mulnor(covar, n, f6, sdev, rhigh, init, iseed, t, work, prob, iter, ifault)
WRITE(*, *)
WRITE(*, *) 'Function F6'
WRITE(*, *) 'ESTIMATED PROBABILITY IS', prob, 'WITH IFAULT = ', ifault
WRITE(*, *) 'ITER = ',iter

STOP


CONTAINS


SUBROUTINE tdiag(n, tol, a, d, e, z, ifault)

IMPLICIT NONE
INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN)      :: tol
REAL, INTENT(IN)      :: a(*)
REAL, INTENT(OUT)     :: d(n)
REAL, INTENT(OUT)     :: e(n)
REAL, INTENT(OUT)     :: z(n,n)
INTEGER, INTENT(OUT)  :: ifault

REAL     :: f, g, h, hh
INTEGER  :: i, j, k, l, i1, j1

REAL, PARAMETER  :: zero = 0.0, one = 1.0

ifault = 1
IF (n <= 1) RETURN
ifault = 0
k = 0
DO  i=1,n
  DO  j=1,i
    k = k+1
    z(i,j) = a(k)
  END DO
END DO
i=n
DO  i1=2,n
  l = i-2
  f = z(i,i-1)
  g=zero
  IF (l < 1) GO TO 25
  DO  k=1,l
    g=g + z(i,k)**2
  END DO
  25 h=g + f*f
  
  IF (g > tol) GO TO 30
  e(i)=f
  d(i)=zero
  GO TO 65

  30 l=l+1
  g=SQRT(h)
  IF (f >= zero) g=-g
  e(i) = g
  h = h-f*g
  z(i,i-1)=f-g
  f=zero
  DO  j=1,l
    z(j,i) = z(i,j)/h
    g = zero
    
    DO  k=1,j
      g = g + z(j,k) *z(i,k)
    END DO
    IF (j >= l)GO TO 47
    j1  = j+1
    DO  k=j1,l
      g=g + z(k,j)*z(i,k)
    END DO
    
    47 e(j)=g/h
    f = f+g*z(j,i)
  END DO
  
  hh = f/(h+h)
  
  DO  j=1,l
    f=z(i,j)
    g=e(j)-hh*f
    e(j)=g
    DO  k=1,j
      z(j,k) = z(j,k) - f*e(k) - g*z(i,k)
    END DO
  END DO
  d(i)=h
  65    i=i-1
END DO
d(1) = zero
e(1) = zero

DO  i=1,n
  l=i-1
  IF (d(i) == zero .OR. l == 0) GO TO 100
  DO  j=1,l
    g=zero
    DO  k=1,l
      g=g + z(i,k)*z(k,i)
    END DO
    DO  k=1,l
      z(k,j) = z(k,j) - g*z(k,i)
    END DO
  END DO
  100   d(i) = z(i,i)
  z(i,i) = one
  IF (l == 0) CYCLE
  DO  j=1,l
    z(i,j) = zero
    z(j,i) = zero
  END DO
END DO
RETURN
END SUBROUTINE tdiag



SUBROUTINE lrvt(n, precis, d, e, z, ifault)

IMPLICIT NONE
INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN)      :: precis
REAL, INTENT(IN OUT)  :: d(n)
REAL, INTENT(OUT)     :: e(n)
REAL, INTENT(IN OUT)  :: z(n,n)
INTEGER, INTENT(OUT)  :: ifault

REAL     :: b, c, f, g, h, p, pr, r, s
INTEGER  :: i, j, k, l, n1, jj, m1, m, i1

INTEGER, PARAMETER  :: mits = 30
REAL, PARAMETER     :: zero = 0.0, one = 1.0, two = 2.0

ifault = 2
IF (n <= 1) RETURN
ifault = 1
n1 = n-1
DO  i=2,n
  e(i-1) = e(i)
END DO
e(n) = zero
b = zero
f = zero
DO  l=1,n
  jj=0
  h = precis * (ABS(d(l)) + ABS(e(l)))
  IF (b < h) b = h
  
  DO  m1=l,n
    m=m1
    IF (ABS(e(m)) <= b) EXIT
  END DO
  IF (m == l) GO TO 85
  40 IF (jj == mits) RETURN
  jj = jj+1
  
  p = (d(l+1) - d(l))/(two*e(l))
  r = SQRT(p*p + one)
  pr = p + r
  IF (p < zero) pr = p - r
  h = d(l) - e(l)/pr
  DO  i=l,n
    d(i) = d(i)-h
  END DO
  f=f+h
  
  p=d(m)
  c=one
  s=zero
  m1=m-1
  i=m
  DO  i1=l,m1
    j=i
    i=i-1
    g=c*e(i)
    h=c*p
    IF(ABS(p) >= ABS(e(i))) GO TO 60
    c = p/e(i)
    r = SQRT(c*c+one)
    e(j) = s*e(i)*r
    s = one/r
    c = c/r
    GO TO 70

    60 c=e(i)/p
    r=SQRT(c*c+one)
    e(j) = s*p*r
    s=c/r
    c=one/r
    70 p=c*d(i)-s*g
    d(j)=h + s*(c*g + s*d(i))
    
    DO  k=1,n
      h=z(k,j)
      z(k,j) = s*z(k,i) + c*h
      z(k,i) = c*z(k,i) - s*h
    END DO
  END DO
  e(l)=s*p
  d(l)=c*p
  IF (ABS(e(l)) > b) GO TO 40
  85 d(l) = d(l) + f
END DO

DO  i=1,n1
  k=i
  p=d(i)
  i1=i+1
  DO  j=i1,n
    IF (d(j) >= p) CYCLE
    k=j
    p=d(j)
  END DO
  IF (k == i) CYCLE
  d(k) = d(i)
  d(i) = p
  DO  j=1,n
    p=z(j,i)
    z(j,i)=z(j,k)
    z(j,k)=p
  END DO
END DO
ifault=0
RETURN
END SUBROUTINE lrvt

END PROGRAM Driver285




FUNCTION f1(n,v) RESULT(fn_val)
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-12-15  Time: 13:39:21

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

REAL, PARAMETER  :: cubsz = 3.0
INTEGER  :: i

!  BOUNDARY OF REGION IS AN N-DIMENSIONAL CUBE.
!  CUBSZ IS THE DISTANCE FROM THE ORIGIN TO THE SIDE OF THE CUBE

fn_val = 0.0
DO  i = 1,n
  fn_val = MAX(fn_val, ABS(v(i)))
END DO
fn_val = fn_val - cubsz
RETURN
END FUNCTION f1



FUNCTION f2(n,v) RESULT(fn_val)

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

REAL     :: b(10)
INTEGER  :: i, j, k

!  THE BOUNDARY IS AN ELLIPSOID OF THE FORM X'AX + BX + C = 0
!  THE UPPER TRIANGLE OF A IS STORED COLUMN-WISE IN PACKED FORM.
!  IN THIS EXAMPLE, A = INVERSE (1  .3), B = 0, C = -1.
!                               (.3  1)

REAL, PARAMETER  :: a(3) = (/ 1.0989011, -.32967033, 1.0989011 /)
REAL, PARAMETER  :: c = -1.0

b(1:2) = 0.0
fn_val = c
k = 0
DO  i = 1,n
  k = k + 1
  fn_val = fn_val + (v(i)**2)*a(k)
  DO  j = i+1,n
    k = k + 1
    fn_val = fn_val + 2.0*v(i)*v(j)*a(k)
  END DO
END DO
fn_val = fn_val + DOT_PRODUCT(b(1:n), v(1:n))
RETURN
END FUNCTION f2



FUNCTION f3(n,v) RESULT(fn_val)

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

INTEGER  :: i
REAL     :: temp

!  RECTANGLE IS GIVEN BY A < V < B FOR EACH COORDINATE.
!  THIS SUBROUTINE MAY ALSO BE USED TO CALCULATE THE CDF BY
!  MAKING THE APPROPRIATE SIDES OF THE RECTANGLE SUFFICIENTLY LARGE.
!  THIS EXAMPLE IS THE FOUR-DIMENSIONAL RECTANGLE DEFINED
!  BY  -INFINITY (ALMOST) < V(I) < 2 FOR EACH I.

REAL, PARAMETER  :: a(4) = (/ -20.0, -20.0, -20.0, -20.0 /)
REAL, PARAMETER  :: b(4) = (/ 2.0, 2.0, 2.0, 2.0 /)

fn_val = -1.0
DO  i=1,n
  IF (v(i) > 0.0) THEN
    temp = v(i) - b(i)
  ELSE
    temp = a(i) - v(i)
  END IF
  fn_val  = MAX(fn_val, temp)
END DO
RETURN
END FUNCTION f3



FUNCTION f4(n,v) RESULT(fn_val)

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

REAL, PARAMETER  :: radsq = 4.0
INTEGER  :: i

!  BOUNDARY IS SPHERE OF RADIUS SQRT(RADSQ).

fn_val = 0.0
DO  i = 1,n
  fn_val = fn_val + v(i)*v(i)
END DO
fn_val = fn_val - radsq
RETURN
END FUNCTION f4



FUNCTION f5(n,v) RESULT(fn_val)

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

INTEGER  :: i, j
REAL     :: temp

!  A CONVEX POLYTOPE WITH P FACES IS DEFINED BY THE INEQUALITIES
!  A(1,I)*V(1) + A(2,I)*V(2) +...+ A(N,I)*V(N) <= 1
!  FOR I=1,...,P.

!  THE FIGURE IN THIS DEMONSTRATION PROGRAM IS THE
!  REGULAR PENTAGON WITH DISTANCE FROM CENTER TO POINT = 1.

INTEGER, PARAMETER  :: p = 5
REAL, PARAMETER     :: a(2,5) = RESHAPE( (/ 0.726542528, -1.0,  &
    1.175570504, 0.381966011, 0.0, 1.236067978, -0.726542528, -1.0,  &
    -1.175570504, 0.381966011 /), (/ 2, 5 /) )

fn_val = -1.0
DO  j=1,n
  fn_val = fn_val + a(j,1)*v(j)
END DO
DO  i=2,p
  temp = -1.0
  DO  j=1,n
    temp = temp + a(j,i)*v(j)
  END DO
  fn_val = MAX(fn_val, temp)
END DO
RETURN
END FUNCTION f5


FUNCTION f6(n,v) RESULT(fn_val)

IMPLICIT NONE
INTEGER, INTENT(IN)  :: n
REAL, INTENT(IN)     :: v(:)
REAL                 :: fn_val

INTEGER  :: i, j, INDEX
REAL     :: x(2), tangnt

!  THE FIGURE IN THIS DEMONSTRATION PROGRAM IS THE
!  REGULAR PENTAGRAM WITH DISTANCE FROM CENTER TO POINT = 2.618.

INTEGER, PARAMETER  :: p = 5
REAL, PARAMETER     :: tang(5) = (/  &
    -1.376381921, -0.3249196962329, 0.3249196962329, 1.376381921, 0.0 /)
REAL, PARAMETER     :: a(2,5) = RESHAPE( (/ -0.726542528, -1.0,  &
    1.175570504, 0.381966011, 0.726542528, -1.0, 0.0, 1.236067978,  &
    1.175570504, 0.381966011 /), (/ 2, 5 /) )

x(1) = ABS(v(1))
x(2) = v(2)
INDEX = p
IF (x(1) == 0.0) THEN
  INDEX = 3
ELSE
  tangnt = x(2)/x(1)
  DO  i=1,p
    IF (tangnt <= tang(i)) THEN
      INDEX = i
      EXIT
    END IF
  END DO
END IF

fn_val = -1.0
DO  j=1,2
  fn_val = fn_val + a(j,INDEX)*x(j)
END DO

RETURN
END FUNCTION f6