MODULE MV_Normal

!      ALGORITHM 725, COLLECTED ALGORITHMS FROM ACM.
 
! Code converted using TO_F90 by Alan Miller
! Date: 2003-06-05  Time: 09:51:33
 
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 19, NO. 4, DECEMBER, 1993, P. 546.

IMPLICIT NONE
INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(14, 60)


CONTAINS


SUBROUTINE dmv(m, k, h, r, prob, eps, ier, ERR, cut)

INTEGER, INTENT(IN)      :: m
INTEGER, INTENT(IN OUT)  :: k
REAL (dp), INTENT(IN)    :: h(20)
REAL (dp), INTENT(IN)    :: r(20,20)
REAL (dp), INTENT(OUT)   :: prob
REAL (dp), INTENT(IN)    :: eps
INTEGER, INTENT(OUT)     :: ier
REAL (dp), INTENT(OUT)   :: ERR
REAL (dp), INTENT(IN)    :: cut


!   DMV CALCULATES THE MULTIVARIATE NORMAL INTEGRAL.

!   THE INTEGRAL LOWER LIMITS ARE -INFINITY FOR ALL VARIABLES.
!   THE UPPER LIMITS ARE GIVEN BY THE VECTOR H (INPUT).
!   THE CORRELATION MATRIX IS R (INPUT).

!   PARAMETERS:

!   M IS THE NUMBER OF VARIABLES.  M <= 20 (INPUT).

!   K IS A FLAG CONTROLLING THE INTEGRATION METHOD (INPUT).
!     2<=K<=10 ONE GAUSIAN QUADRATURE WITH K POINTS FOR EACH
!              VARIABLE IS PERFORMED.  EPS IS DISREGARDED.
!    12<=K<=20 PROGRESSIVE GAUSIAN QUADRATURES ARE PERFORMED FOR
!              K=2,3,...,KMAX, WHERE KMAX IS 10 LESS THAN THE INPUT
!              VALUE OF K, UNTIL THE DIFFERENCE BETWEEN THE CALCULATIONS
!              FOR TWO SUCCESSIVE K'S DOES NOT EXCEED EPS.
!              THE LAST VALUE OF K IS RETURNED IN IER.
!              IF KMAX IS REACHED BEFORE A DIFFERENCE OF EPS IS ACHIEVED,
!              THEN IER=0 AND THIS DIFFERENCE IS RETURNED IN ERR.
!     EVERY OTHER K IS EQUIVALENT TO K=20.

!   H IS THE VECTOR OF UPPER LIMITS FOR THE INTEGRATION (INPUT).

!   R IS THE CORRELATION MATRIX DIMENSIONED R(20,20) (INPUT).
!     IT MUST BE A SYMMETRIC, POSITIVE SEMI-DEFINITE MATRIX.
!     THE CODE DOES NOT CHECK FOR THAT.
!     IF R IS SINGULAR, NUMERICAL PROBLEMS MAY OCCUR.

!   PROB IS THE CALCULATED PROBABILITY (OUTPUT).

!   EPS IS THE REQUIRED ERROR LIMIT (INPUT).

!   IER IS AN ERROR INDICATOR. ON NORMAL TERMINATION 0 < IER <= 10 IS THE
!       VALUE OF K CORRESPONDING TO THE RETURNED VALUE PROB, (OUTPUT).
!       IER= 0 IS RETURNED WHEN THE ERROR IS GREATER THAN EPS.
!       IER=98 MEANS AN ILLEGAL M VALUE (M SHOULD NOT EXCEED 20).
!       IER=99 MEANS THAT R IS SINGULAR OR ITS DETERMINANT IS NEGATIVE.

!   ERR IS THE DIFFERENCE IN CALCULATING BY K-1 AND K. CALCULATED
!       ONLY WHEN THE INPUT K > 10 (OUTPUT).

!   CUT - IF X<-CUT THEN EXP(X) IS IGNORED (ASSUMED ZERO).
!         SUGGESTED VALUE CUT=15 (INPUT).

!   TWO SUBROUTINES (DPOFA, DPODI) ARE TAKEN FROM THE NAAS:LINPACK
!   LIBRARY. THEY CAN BE REPLACED BY A SUBROUTINE THAT INVERTS A
!   SYMMETRIC MATRIX AND CALCULATES ITS DETERMINANT.

INTEGER    :: ker, kmax
REAL (dp)  :: pold

IF(k < 2 .OR. k > 10) GO TO 10

!     DMV1 CALCULATES THE PROBABILITY FOR A GIVEN K

CALL dmv1(m, k, h, r, cut, prob, ier)
IF(ier == 0) ier=k
RETURN

10   kmax=10
IF(k > 11 .AND. k < 21) kmax = k - 10
pold=-1.
DO  ker=2,kmax
  CALL dmv1(m, ker, h, r, cut, prob, ier)
  IF(ier > 0) RETURN
  ier=ker
  ERR=prob-pold
  IF(ABS(ERR) <= eps) RETURN
  pold=prob
END DO
ier=0
RETURN
END SUBROUTINE dmv


SUBROUTINE dmv1(m, k, h, r, cut, res, ier)

!   DMV1 CALCULATES THE PROBABILITY FOR A GIVEN K BY APPLYING THEOREM 1.
!   IT CALLS SUBROUTINE SOLVE TO CALCULATE INDIVIDUAL VALUES BY (1) and (2).
!   DMV1 CAN BE CALLED DIRECTLY.

!   K IS THE NUMBER OF POINTS IN THE GAUSS QUADRATURE (INPUT).
!       IF K < 2 OR K > 10 THEN K IS SET TO 10.

!   IER IS AN ERROR INDICATOR.
!       IER=0 IS NORMAL TERMINATION.
!       IER=98 INDICATES THAT M IS OUT OF RANGE.
!       IER=99 MEANS THAT R IS SINGULAR OR ITS DETERMINANT IS NEGATIVE.

!   THE OTHER VARIABLES ARE DESCRIBED ABOVE IN DMV.

INTEGER, INTENT(IN)      :: m
INTEGER, INTENT(IN OUT)  :: k
REAL (dp), INTENT(IN)    :: h(20)
REAL (dp), INTENT(IN)    :: r(20,20)
REAL (dp), INTENT(IN)    :: cut
REAL (dp), INTENT(OUT)   :: res
INTEGER, INTENT(OUT)     :: ier

REAL (dp)  :: h1(20), r1(20,20), s
INTEGER    :: i, ipos, j, l, l1, l2, list(20), list1(20)

ier=0
IF(m < 0 .OR. m > 20) ier=98
IF(ier > 0) RETURN
res=0.

!   WE EXPAND UP TO 2 TO THE POWER OF M TERMS DESCRIBED IN THEOREM 1.
!   THE VECTOR LIST REPRESENTS THE "+" "-" AND INFINITY IN THE TERM.
!   LIST(I) HAS THE VALUE OF 0 IF H(I) <= 0, 1 IF H(I) > 0, AND
!   THE VALUE OF 2 IF INFINITY IS USED AS THE UPPER LIMIT AND THE
!   VARIABLE IS DROPPED FROM THE INTEGRATION.
!   LIST1 IS THE LIST OF VARIABLES THAT PARTICIPATE IN THE
!   INTEGRATION (THAT DO NOT HAVE INFINITY AS AN UPPER LIMIT).

DO  i=1,m
  list(i)=0
  IF(h(i) > 0.) list(i)=1
END DO
20   ipos=0
l=0
DO  i=1,m
  IF(list(i) == 2) CYCLE
  l=l+1
  list1(l)=i
END DO
IF(l == 0) GO TO 50
DO  i=1,l
  l1=list1(i)
  h1(i)=h(l1)
  IF(list(l1) == 1) h1(i)=-h1(i)
  IF(list(l1) == 1) ipos=1-ipos
  DO  j=1,l
    l2=list1(j)
    r1(i,j)=r(l1,l2)
    IF(list(l1) == 1) r1(i,j)=-r1(i,j)
    IF(list(l2) == 1) r1(i,j)=-r1(i,j)
  END DO
END DO

50 CALL solve(l, k, h1, r1, cut, s, ier)
IF(ier > 0) RETURN
IF(ipos == 1) s = -s
res = res + s
DO  i=1,m
  IF(list(i) == 0) CYCLE
  list(i) = list(i) + 1
  IF(list(i) == 2) GO TO 20
  list(i) = 1
END DO
RETURN
END SUBROUTINE dmv1


SUBROUTINE solve(m, k, h, r, cut, solution, ier)

!     SOLVE CALCULATES THE MULTIVARIATE PROBABILITY OF A GIVEN VECTOR H
!     AND GIVEN K.  SOLVE CAN BE CALLED DIRECTLY.

!     ALL THE ARGUMENTS ARE AS DESCRIBED ABOVE IN DMV1 OR DMV.

!     THE INTEGRATION PARAMETERS TAKEN FROM [8] HAVE 15 DIGITS ACCURACY.

INTEGER, INTENT(IN)      :: m
INTEGER, INTENT(IN OUT)  :: k
REAL (dp), INTENT(IN)    :: h(20)
REAL (dp), INTENT(IN)    :: r(20,20)
REAL (dp), INTENT(IN)    :: cut
REAL (dp), INTENT(OUT)   :: solution
INTEGER, INTENT(OUT)     :: ier

REAL (dp)  :: d(2), det, h1(20), h2(20), pi, prod, r1(20,20), sum, y(20)
INTEGER    :: list(20)
INTEGER    :: i, j, l

!     THE VECTORS COEF AND X CONTAIN THE GAUSS QUADRATURE COEFFICIENTS
!     AND POINTS RESPECTIVELY AND ARE OBTAINED FROM [6].

REAL(dp), PARAMETER  :: coef(10,10) = RESHAPE( (/ 0.0D0, 0.0D0, 0.0D0, 0.0D0, &
    0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .6405291796843790D0,  &
    .245697745768379D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,  &
    0.0D0, .446029770466658D0, .396468266998335D0, 4.37288879877644D-2,  &
    0.d0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .325302999756919D0,  &
    .421107101852062D0, .133442500357520D0, .637432348625728D-2, 0.d0, 0.0D0, &
    0.0D0, 0.0D0, 0.0D0, 0.0D0, .248406152028443D0, .392331066652399D0,   &
    .211418193076057D0, .332466603513439D-1, .824853344515628D-3, 0.0D0,  &
    0.0D0, 0.0D0, 0.0D0, 0.0D0, .196849675488598D0, .349154201525395D0, &
    .257259520584421D0, .0760131375840057D0, .685191862513596D-2,  &
    9.84716452019267D-5, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .160609965149261D0,  &
    .306319808158099D0, .275527141784905D0, .120630193130784D0,   &
    .0218922863438067D0, .00123644672831056D0, 1.10841575911059D-5, 0.0D0,  &
    0.0D0, 0.0D0, .13410918845336D0, .26833075447264D0, .275953397988422D0,  &
    .15744828261879D0, .0448141099174625D0, .00536793575602526D0,   &
    2.02063649132407D-4, 1.19259692659532D-6, 0.d0, 0.0D0, .114088970242118D0, &
    .235940791223685D0, .266425473630253D0, .183251679101663D0,   &
    .0713440493066916D0, .0139814184155604D0, .00116385272078519D0,  &
    .305670214897831D-4, 1.23790511337496D-7, 0.d0, 0.0985520975191087D0,  &
    .208678066608185D0, .252051688403761D0, .198684340038387D0,  &
    .097198422760062D0, .0270244164355446D0, .00380464962249537D0,   &
    2.28886243044656D-4, 4.34534479844469D-6, 1.24773714817825D-8 /),  &
    (/ 10, 10 /) )
REAL(dp), PARAMETER  :: x(10,10) = RESHAPE( (/ 0.0D0, 0.0D0, 0.0D0, 0.0D0, &
    0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .300193931060839D0,  &
    1.25242104533372D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,  &
    0.0D0, .190554149798192D0, .848251867544577D0, 1.79977657841573D0,  &
    0.d0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .133776446996068D0,  &
    .624324690187190D0, 1.34253782564499D0, 2.26266447701036D0, 0.d0,  &
    0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .100242151968216D0,  &
    .482813966046201D0, 1.06094982152572D0, 1.77972941852026D0,  &
    2.66976035608766D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0, .0786006594130979D0,&
    .386739410270631D0, .866429471682044D0, 1.46569804966352D0,   &
    2.172707796939D0, 3.03682016932287D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0,   &
    .0637164846067008D0, .318192018888619D0, .724198989258373D0,   &
    1.23803559921509D0, 1.83852822027095D0, 2.53148815132768D0,   &
    3.37345643012458D0, 0.0D0, 0.0D0, 0.0D0, .0529786439318514D0,  &
    .267398372167767D0, .616302884182402D0, 1.06424631211623D0,  &
    1.58885586227006D0, 2.18392115309586D0, 2.86313388370808D0,  &
    3.6860071627244D0, 0.d0, 0.0D0, .0449390308011934D0, .228605305560535D0, &
    .532195844331646D0, .927280745338081D0, 1.39292385519588D0,  &
    1.91884309919743D0, 2.50624783400574D0, 3.17269213348124D0,  &
    3.97889886978978D0, 0.d0, .0387385243257289D0, .198233304013083D0,  &
    .465201111814767D0, .816861885592273D0, 1.23454132402818D0,  &
    1.70679814968913D0, 2.22994008892494D0, 2.80910374689875D0,  &
    3.46387241949586D0, 4.25536180636608D0 /), (/ 10, 10 /) )

pi = 4.0D0*ATAN(1.d0)
ier=0
IF(m < 0 .OR. m > 20) ier=98
IF(ier > 0) RETURN
IF(k < 2 .OR. k > 10) k=10
IF(m > 0) GO TO 10
solution=1.
ier=0
RETURN

10 r1(1:m,1:m) = r(1:m,1:m)
CALL dpofa(r1, 20, m, ier)
IF(ier > 0) ier=99
IF(ier == 99) RETURN
CALL dpodi(r1, 20, m, d, 11)
IF(d(1) <= 0. .OR. d(2) <= -30) ier=99
IF(ier == 99) RETURN

!   NOTE THAT ANOTHER SUBROUTINE MAY GIVE THE DETERMINANT DIRECTLY
!   IN DET AND THEREFORE THE NEXT STATEMENT MAY BE UNNECESSARY.
!   R1 CONTAINS THE INVERSE OF R.  LOOP 30 CATERS TO THE POSSIBILITY
!   THAT THE INVERSE IS STORED ONLY IN THE LOWER TRIANGLE.

det=d(1)*10.d0**d(2)
DO  i=1,m
  DO  j=i,m
    r1(j,i)=r1(i,j)
  END DO
END DO
prod=1.d0
DO  i=1,m
  h1(i)=SQRT(2./r1(i,i))
  h2(i)=h(i)/h1(i)
  prod=prod*pi*r1(i,i)
  list(i)=1
END DO
prod=1.d0 / SQRT(det*prod)
DO  i=1,m
  DO  j=1,m
    r1(i,j)=r1(i,j)*h1(i)*h1(j)*0.5D0
  END DO
END DO

solution=0.
60   sum=0.
DO  i=1,m
  l=list(i)
  sum=sum + x(l,k)*x(l,k)
  y(i) = x(l,k) - h2(i)
END DO
DO  i=1,m
  DO  j=1,m
    sum = sum - y(i)*y(j)*r1(i,j)
  END DO
END DO
IF(sum < -cut) GO TO 100
sum=EXP(sum)
DO  i=1,m
  sum=sum*coef(list(i),k)
END DO
solution = solution + sum
100  DO  i=1,m
  list(i)=list(i)+1
  IF(list(i) <= k) GO TO 60
  list(i)=1
END DO
solution = solution*prod
RETURN
END SUBROUTINE solve


SUBROUTINE dpofa(a, lda, n, info)

INTEGER, INTENT(IN)        :: lda
REAL (dp), INTENT(IN OUT)  :: a(lda,:)
INTEGER, INTENT(IN)        :: n
INTEGER, INTENT(OUT)       :: info

!   dpofa factors a double precision symmetric positive definite matrix.

!   dpofa is usually called by dpoco, but it can be called
!   directly with a saving in time if  rcond  is not needed.
!   (time for dpoco) = (1 + 18/n)*(time for dpofa) .

!   on entry

!      a       double precision(lda, n)
!              the symmetric matrix to be factored.  Only the
!              diagonal and upper triangle are used.

!      lda     integer
!              the leading dimension of the array  a .

!      n       integer
!              the order of the matrix  a .

!   on return

!      a       an upper triangular matrix  r  so that  a = trans(r)*r
!              where  trans(r)  is the transpose.
!              the strict lower triangle is unaltered.
!              if  info .ne. 0 , the factorization is not complete.

!      info    integer
!              = 0  for normal return.
!              = k  signals an error condition.  The leading minor
!                   of order  k  is not positive definite.

!   linpack.  this version dated 08/14/78 .
!   cleve moler, university of new mexico, argonne national lab.

!   subroutines and functions

!   blas ddot
!   fortran dsqrt

!   internal variables

REAL (dp)  :: s, t
INTEGER    :: j, jm1, k

!     begin block with ...exits to 40

DO  j = 1, n
  info = j
  s = 0.0D0
  jm1 = j - 1
  IF (jm1 < 1) GO TO 20
  DO  k = 1, jm1
    t = a(k,j) - ddot(k-1, a(1:,k), 1, a(1:,j), 1)
    t = t/a(k,k)
    a(k,j) = t
    s = s + t*t
  END DO
  20 s = a(j,j) - s
!     ......exit
  IF (s <= 0.0D0) GO TO 40
  a(j,j) = SQRT(s)
END DO
info = 0
40 RETURN
END SUBROUTINE dpofa


SUBROUTINE dpodi(a, lda, n, det, job)

INTEGER, INTENT(IN)        :: lda
REAL (dp), INTENT(IN OUT)  :: a(lda,:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(OUT)     :: det(2)
INTEGER, INTENT(IN)        :: job

!   dpodi computes the determinant and inverse of a certain
!   double precision symmetric positive definite matrix (see below)
!   using the factors computed by dpoco, dpofa or dqrdc.

!   on entry

!      a       double precision(lda, n)
!              the output  a  from dpoco or dpofa
!              or the output  x  from dqrdc.

!      lda     integer
!              the leading dimension of the array  a .

!      n       integer
!              the order of the matrix  a .

!      job     integer
!              = 11   both determinant and inverse.
!              = 01   inverse only.
!              = 10   determinant only.

!   on return

!      a       if dpoco or dpofa was used to factor  a  then
!              dpodi produces the upper half of inverse(a) .
!              if dqrdc was used to decompose  x  then
!              dpodi produces the upper half of inverse(trans(x)*x)
!              where trans(x) is the transpose.
!              elements of  a  below the diagonal are unchanged.
!              if the units digit of job is zero,  a  is unchanged.

!      det     double precision(2)
!              determinant of  a  or of  trans(x)*x  if requested.
!              otherwise not referenced.
!              determinant = det(1) * 10.0**det(2)
!              with  1.0 .le. det(1) .lt. 10.0
!              or  det(1) .eq. 0.0 .

!   error condition

!      a division by zero will occur if the input factor contains
!      a zero on the diagonal and the inverse is requested.
!      it will not occur if the subroutines are called correctly
!      and if dpoco or dpofa has set info .eq. 0 .

!   linpack.  this version dated 08/14/78 .
!   cleve moler, university of new mexico, argonne national lab.

!   subroutines and functions

!   blas daxpy,dscal
!   fortran mod

!   internal variables

REAL (dp)  :: s, t
INTEGER    :: i, j, jm1, k, kp1

!     compute determinant

IF (job/10 == 0) GO TO 70
det(1) = 1.0D0
det(2) = 0.0D0
s = 10.0D0
DO  i = 1, n
  det(1) = a(i,i)**2*det(1)
!        ...exit
  IF (det(1) == 0.0D0) EXIT
  10  IF (det(1) >= 1.0D0) GO TO 30
  det(1) = s*det(1)
  det(2) = det(2) - 1.0D0
  GO TO 10
  30  IF (det(1) < s) CYCLE
  det(1) = det(1)/s
  det(2) = det(2) + 1.0D0
  GO TO 30
END DO

!     compute inverse(r)

70 IF (MOD(job,10) == 0) GO TO 140
DO  k = 1, n
  a(k,k) = 1.0D0/a(k,k)
  t = -a(k,k)
  CALL dscal(k-1, t, a(1:,k), 1)
  kp1 = k + 1
  IF (n >= kp1) THEN
    DO  j = kp1, n
      t = a(k,j)
      a(k,j) = 0.0D0
      CALL daxpy(k, t, a(1:,k), 1, a(1:,j), 1)
    END DO
  END IF
END DO

!        form  inverse(r) * trans(inverse(r))

DO  j = 1, n
  jm1 = j - 1
  IF (jm1 >= 1) THEN
    DO  k = 1, jm1
      t = a(k,j)
      CALL daxpy(k, t, a(1:,j), 1, a(1:,k), 1)
    END DO
  END IF
  t = a(j,j)
  CALL dscal(j, t, a(1:,j), 1)
END DO

140 RETURN
END SUBROUTINE dpodi


SUBROUTINE daxpy(n, da, dx, incx, dy, incy)

!     constant times a vector plus a vector.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: da
REAL (dp), INTENT(IN)   :: dx(:)
INTEGER, INTENT(IN)     :: incx
REAL (dp), INTENT(OUT)  :: dy(:)
INTEGER, INTENT(IN)     :: incy

INTEGER  :: i, ix, iy, m, mp1

IF(n <= 0) RETURN
IF (da == 0.0D0) RETURN
IF(incx == 1 .AND. incy == 1) GO TO 20

!        code for unequal increments or equal increments not equal to 1

ix = 1
iy = 1
IF(incx < 0) ix = (-n+1)*incx + 1
IF(incy < 0) iy = (-n+1)*incy + 1
DO  i = 1,n
  dy(iy) = dy(iy) + da*dx(ix)
  ix = ix + incx
  iy = iy + incy
END DO
RETURN

!        code for both increments equal to 1

!        clean-up loop

20 m = MOD(n,4)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  dy(i) = dy(i) + da*dx(i)
END DO
IF( n < 4 ) RETURN
40 mp1 = m + 1
DO  i = mp1,n,4
  dy(i) = dy(i) + da*dx(i)
  dy(i + 1) = dy(i + 1) + da*dx(i + 1)
  dy(i + 2) = dy(i + 2) + da*dx(i + 2)
  dy(i + 3) = dy(i + 3) + da*dx(i + 3)
END DO
RETURN
END SUBROUTINE daxpy


SUBROUTINE dscal(n, da, dx, incx)

!     scales a vector by a constant.
!     uses unrolled loops for increment equal to one.
!     jack dongarra, linpack, 3/11/78.

INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(IN)   :: da
REAL (dp), INTENT(OUT)  :: dx(:)
INTEGER, INTENT(IN)     :: incx

INTEGER :: i, m, mp1, nincx

IF(n <= 0) RETURN
IF(incx == 1) GO TO 20

!        code for increment not equal to 1

nincx = n*incx
DO  i = 1,nincx,incx
  dx(i) = da*dx(i)
END DO
RETURN

!        code for increment equal to 1


!        clean-up loop

20 m = MOD(n,5)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  dx(i) = da*dx(i)
END DO
IF( n < 5 ) RETURN
40 mp1 = m + 1
DO  i = mp1,n,5
  dx(i) = da*dx(i)
  dx(i + 1) = da*dx(i + 1)
  dx(i + 2) = da*dx(i + 2)
  dx(i + 3) = da*dx(i + 3)
  dx(i + 4) = da*dx(i + 4)
END DO
RETURN
END SUBROUTINE  dscal


FUNCTION ddot(n, dx, incx, dy, incy) RESULT(fn_val)

!     forms the dot product of two vectors.
!     uses unrolled loops for increments equal to one.
!     jack dongarra, linpack, 3/11/78.

INTEGER, INTENT(IN)    :: n
REAL (dp), INTENT(IN)  :: dx(:)
INTEGER, INTENT(IN)    :: incx
REAL (dp), INTENT(IN)  :: dy(:)
INTEGER, INTENT(IN)    :: incy
REAL (dp)              :: fn_val

REAL (dp)  :: dtemp
INTEGER    :: i, ix, iy, m, mp1

fn_val = 0.0D0
dtemp = 0.0D0
IF(n <= 0) RETURN
IF(incx == 1 .AND. incy == 1) GO TO 20

!        code for unequal increments or equal increments
!          not equal to 1

ix = 1
iy = 1
IF(incx < 0)ix = (-n+1)*incx + 1
IF(incy < 0)iy = (-n+1)*incy + 1
DO  i = 1,n
  dtemp = dtemp + dx(ix)*dy(iy)
  ix = ix + incx
  iy = iy + incy
END DO
fn_val = dtemp
RETURN

!        code for both increments equal to 1

!        clean-up loop

20 m = MOD(n,5)
IF( m == 0 ) GO TO 40
DO  i = 1,m
  dtemp = dtemp + dx(i)*dy(i)
END DO
IF( n < 5 ) GO TO 60
40 mp1 = m + 1
DO  i = mp1,n,5
  dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +  &
      dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
END DO
60 fn_val = dtemp
RETURN
END FUNCTION ddot

END MODULE MV_Normal



PROGRAM test725

!   THIS TEST PROGRAM CHECKS DMV ON A SET OF PRE-SPECIFIED PROBLEMS.
!   FOR EACH PROBLEM, THE RESULT CALCULATED BY THE PROGRAM IS PRINTED
!   NEXT TO THE CORRECT RESULT SO THAT ONE CAN CHECK HIS MACHINE.
!   NO INPUT IS REQUIRED.

USE MV_Normal
IMPLICIT NONE

REAL (dp)  :: diff, err, h(20), prob, r(20,20)
INTEGER    :: ier, ih, ir, j, k, l, m, number
REAL(dp), PARAMETER  :: result(12) = (/ .72437, .74520, .95618, .95855,   &
                                        .63439, .67778, .93656, .94253,   &
                                        .56299, .62670, .91819, .92845 /)

!       All output is to the terminal.  To modify this, change all
!       following "WRITE(*" to "WRITE(LUN" and precede this comment
!       by code to initialize unit LUN for output.

WRITE(*,101)
101  FORMAT('   M    CALCULATED     TRUE       DIFF.    IER')
NUMBER=0
k=0
DO  m=2,4
  DO  ih=1,2
    DO  ir=1,2
      DO  l=1,m
        h(l)=ih
        DO  j=1,m
          r(l,j)=0.25*ir
        END DO
        r(l,l)=1.
      END DO
      NUMBER=NUMBER+1
      CALL dmv(m, k, h, r, prob, 1.d-5, ier, ERR, 15.d0)
      diff=prob - result(NUMBER)
      WRITE(*, 34) m, prob, result(NUMBER), diff, ier
      34  FORMAT(i4, 2F12.5, ' ', g12.2, i5)
    END DO
  END DO
END DO
STOP
END PROGRAM test725