subroutine simdo ( qind, qfor, iprod, kdim, jsub, ivec, ifault )

!*****************************************************************************80
!
!! SIMDO generates multi-indices, simulating nested DO-loops.
!
!  Discussion:
!
!    The loops are assumed to be nested to a depth of K.
!
!    The R-th loop is assumed to have upper limit N(R) and increment Inc(R).
!
!    The total number of executions of the innermost loop is 
!
!      N = product ( 1 <= R <= K ) N(R).
!
!    Let these executions be indexed by the single integer ( kind = 4 ) J, which
!    we call the index subscript.
!
!    Each value of J corresponds to a particular set of loop indices,
!    which we call the subscript vector I(J).
!
!    This routine can start with J and find I(J), or determine
!    J from I(J).
!    
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    26 July 2008
!
!  Author:
!
!    M O'Flaherty, G MacKenzie
!
!  Reference:
!
!    M O'Flaherty, G MacKenzie,
!    Algorithm AS 172:
!    Direct Simulation of Nested Fortran DO-LOOPS,
!    Applied Statistics,
!    Volume 31, Number 1, 1982, pages 71-74.
!
!  Parameters:
!
!    Input, logical QIND.
!    TRUE to convert an index subscript J to the subscript vector I(J).
!    FALSE to convert the subscript vector I(J) to the index subscript J.
!
!    Input, logical QFOR,
!    TRUE if conversion is required in standard Fortran subscripting order,
!    FALSE otherwise.
!
!    Input, integer ( kind = 4 ) IPROD(KDIM), contains the partial products.
!    If QFOR is FALSE, then
!      IPROD(S) = product ( 1 <= R <= S ) N(R).
!    If QFOR is TRUE, then
!      IPROD(S) = product ( 1 <= R <= S ) N(KDIM+1-R).
!
!    Input, integer ( kind = 4 ) KDIM, the nesting depth of the loops.
!
!    Input/output, integer ( kind = 4 ) JSUB.
!    If QIND is TRUE, then JSUB is an input quantity, an index subscript J
!    to be converted into the subscript vector I(J).
!    If QIND is FALSE, then JSUB is an output quantity, the index subscript J
!    corresponding to the subscript vector I(J).
!
!    Input/output, integer ( kind = 4 ) IVEC(KDIM).
!    if QIND is TRUE, then IVEC is an output quantity, the subscript vector I(J)
!    corresponding to the index subscript J.
!    If QIND is FALSE, then IVEC is an input quantity, a subscript vector I(J)
!    for which the corresponding index subscript J is to be computed.
!
!    Output, integer ( kind = 4 ) IFAULT, error flag.
!    0, no error was detected.
!    1, if QIND is TRUE, and the input value of JSUB exceeds IPROD(KDIM).
!    2, if QIND is FALSE, and IVEC contains an illegal component.
!
  implicit none

  integer ( kind = 4 ) kdim

  integer ( kind = 4 ) i
  integer ( kind = 4 ) ifault
  integer ( kind = 4 ) ik
  integer ( kind = 4 ) iprod(kdim)
  integer ( kind = 4 ) itempv
  integer ( kind = 4 ) ivec(kdim)
  integer ( kind = 4 ) jsub
  logical qfor
  logical qind

  ifault = 0
!
!  Index subscript to subscript vector conversion.
!
  if ( qind )  then

    if ( iprod(kdim) < jsub ) then
      ifault = 1
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SIMDO - Fatal error!'
      write ( *, '(a)' ) '  JSUB is out of bounds.'
      stop
    end if

    itempv = jsub - 1

    do i = 1, kdim - 1
      ik = kdim - i
      ivec(i) = itempv / iprod(ik)
      itempv = itempv - iprod(ik) * ivec(i)
      ivec(i) = ivec(i) + 1
    end do

    ivec(kdim) = itempv + 1
    if ( qfor ) then
      ivec(1:kdim) = ivec(kdim:1:-1)
    end if
!
!  Subscript vector to index subscript conversion.
!
  else

    if ( .not. qfor ) then
      ivec(1:kdim) = ivec(kdim:1:-1)
    end if

    if ( iprod(1) < ivec(1) ) then
      ifault = 2
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SIMDO - Fatal error!'
      write ( *, '(a)' ) '  An entry of IVEC is out of bounds.'
      stop
    end if

    do i = 2, kdim
      if ( iprod(i) / iprod(i-1) < ivec(i) ) then
        ifault = 2
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SIMDO - Fatal error!'
        write ( *, '(a)' ) '  An entry of IVEC is out of bounds.'
        stop
      end if
    end do

    jsub = ivec(1)

    do i = 2, kdim
      jsub = jsub + ( ivec(i) - 1 ) * iprod(i-1)
    end do
!
!  As a courtesy to the caller, UNREVERSE the IVEC vector
!  if you reversed it!
!
    if ( .not. qfor ) then
      ivec(1:kdim) = ivec(kdim:1:-1)
    end if

  end if

  return
end