!   Program for computing the control limit of a MEWMC
!   User input -- dimension, smoothing constant, ARL target, and
!   required number of simulations.
!
!   Output -- estimated control limit with standard error
!
!   Subroutine required:
!   The program calls 'randn', a random normal generator.  The call
!   randn(data, length, seed1, seed2) fills 'length' elements of the
!   real(8) vector 'data' with pseudorandom N(0,1) quantities.
!   'seed1' and 'seed2' are two random number generator seeds.
!   Users will have to replaced this call with one to their own
!   random normal generator.
!

    implicit none

    integer, parameter :: grid=2000

    real(8), allocatable :: data(:), eye(:,:), update(:,:), &
                            covinv(:,:), dotter(:), vars(:)

    real(4), allocatable :: sucmax(:)

    real(8) :: h, lambda, onelam, logonelam, crit, lndet, &
               one=1.d0, zero=0.d0, range,   &
               arl(grid), serl(grid), &
               dev(grid), gmax, count, big=1.e20, c, &
               detupd, target, guessh, guessse, step

    real(4) :: timenow, lasttime

    integer :: i, j, isim, nsim, nord, ij=0, kl=0, &
               ios, inx, last, lastx, maxrun


    logical :: hit

    write(*,*) 'Program for control limit of a MEWMC control chart'

! forever:  do
      write(*,*) ' Enter dimensions, lambda, target ARL and # of sims :   '
      read(*,*, iostat=ios) nord, lambda, target, nsim

      if (ios /= 0) then
        write(*,*) 'Invalid entry.  Exiting'
        pause
        stop
        end if

      if (lambda <= zero .or. lambda >=  one) then
        write(*,*) 'Lambda must be between 0 and 1'
        pause
        stop
        end if

      if (nsim < 20) then
        write(*,*) '20 simulations minimum'
        pause
        stop
        end if

      write(*,*) 'Starting run.... '

      maxrun = 200 + 14 * target       ! Covers all but 1 run in a million

      allocate (data(nord), eye(nord,nord), covinv(nord,nord), &
        dotter(nord), update(nord,nord), vars(nord), sucmax(0:maxrun))

      onelam = one - lambda
      logonelam = log(onelam)
      c = lambda / onelam
      eye = 0
      do i = 1, nord
        eye(i,i) = one
        end do

      last = maxrun

      gmax = big
      range = -big
      arl = 0
      serl = 0
      count = 0

      call cpu_time(lasttime)

main:       do isim = 1, nsim
        call cpu_time(timenow)
        if (timenow > lasttime + 60) then
          write(*,'('' Doing iteration'',i7)') isim
          lasttime = timenow
          end if
        vars = one
        covinv = eye
        sucmax = -big
        lndet = zero

        do i = 1, last
          call randn(data, nord, ij, kl)
          vars = onelam * vars + lambda * data**2
          dotter = matmul(covinv, data)
          update = matmul(reshape(dotter, (/ nord, 1 /)), reshape(dotter, (/ 1, nord /)))
          detupd = one + c * dot_product(dotter, data)
          covinv = (covinv - c * update / detupd) / onelam
          if (detupd <= zero) then
            write(*,'('' loop'',2i6,'' lndet and update'',2g15.7)') isim, i, &
              lndet, detupd
            write(*,'('' Precision lost; can not handle these settings'')')
            stop
            end if
          lndet = lndet + log(detupd) + nord * logonelam
          crit = -nord - lndet + sum(vars)
          sucmax(i) = max(sucmax(i-1), crit)
          if (isim <= 10) range = max(range, crit)
          if (crit > gmax) exit
          end do

        select case (isim)

          case (11)
          gmax = range

          case (12:)
          lastx = 0
          count = count + 1
          do i = 1, last
            inx = max(zero, grid * min(one, sucmax(i) / range)) + one
            inx = max(1, min(inx, grid))
            if (sucmax(i) > sucmax(i-1)) then
              dev(lastx+1:inx) = i - arl(lastx+1:inx)
              arl(lastx+1:inx) = arl(lastx+1:inx) + dev(lastx+1:inx) / count
              serl(lastx+1:inx) = serl(lastx+1:inx) + (count-1)/count * &
                dev(lastx+1:inx) ** 2
              lastx = inx
              end if
            end do
          end select
!
!         do some ranging
!
        if (isim > 30 .and. mod(isim,15) == 0) then
          do i = 1, grid
            if (arl(i) > target * (1.1+4./sqrt(count-10))) then
              gmax = min(gmax, (i-1)*range / grid)
              end if
            end do
          end if
        end do main

      step = range / grid
      hit = .false.
      do i = 2, grid
        h = (i-1) * step
        if (h <= gmax) then
          serl(i) = sqrt(serl(i) / (count * (count - one)))

          if (mod(i,10) == 0 .or. arl(i) > target) &
            write( *,'('' p, lambda, h'',i5, f6.3, f10.4,'' ARL'', 2f10.2)') &
            nord, lambda, h, arl(i), serl(i)

          if (arl(i) > target .and. arl(i-1) <= target) then  !  Straddled
            guessh = h - (target - arl(i)) / (arl(i-1) - arl(i)) * step
            guessse = serl(i) * 20 * step / (arl(i) - arl(i-20))
            hit = .true.
            exit
            end if
          end if
        end do

      if (hit) then
        write(*,'('' control limit, s.e.'',2f10.4, '',   s.e. of ARL'',f9.1)') &
          guessh, guessse, serl(i)
        else
        write(*,'('' run failed; target ARL exceeded program dimensions.'', &
          &''  Initial range'',2f9.4)') range, gmax
        end if

      deallocate (data, eye, covinv, dotter, update, vars, sucmax)

!     end do forever

    pause

    end program