subroutine array_size ( rank, dim, size )
!
!*******************************************************************************
!
!! ARRAY_SIZE computes the "size" of an object of given shape.
!
!
!  Discussion:
!
!    If A is an unpacked array, then it has both maximum dimensions MAXDIM and
!    used dimensions DIM.  Calling
!
!      array_size ( rank, maxdim, maxsize )
!
!    will return the maximum number of entries that A could hold, while
!
!      array_size ( rank, dim, size )
!
!    will return the actual number of entries that A is being used to hold.
!
!    If A is a "packed array", then the maximum and used dimensions are
!    the same.
!
!  Modified:
!
!    21 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANK, the "rank" of the object.  The rank is the number of
!    indices required to access any particular entry.  A scalar has rank 0,
!    a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer DIM(RANK), the "maximum" or the "used" dimensions.
!    The I-th index of the object must range between 1 and DIM(I).
!
!    Output, integer SIZE, depending on whether the input DIM was the
!    "maximum" or "used" dimensions of the object, SIZE will be the 
!    "maximum" or "used" size of the object.
!
  implicit none
!
  integer rank
!
  integer dim(rank)
  integer i
  integer size
!
!  Watch out for nasty data.
!
  size = 0
  do i = 1, rank
    if ( dim(i) <= 0 ) then
      return
    end if
  end do

  size = product ( dim(1:rank) )

  return
end
subroutine index_next ( n, hi, iarray, more )
!
!*******************************************************************************
!
!! INDEX_NEXT generates all indices within given upper limits.
!
!
!  Example:
!
!    N = 3,
!    HI(1) = 4, HI(2) = 2, HI(3) = 3
!
!    1   2   3
!    ---------
!    1   1   1
!    2   1   1
!    3   1   1
!    4   1   1
!    1   2   1
!    2   2   1
!    3   2   1
!    4   2   1
!    1   1   2
!    2   1   2
!    3   1   2
!    4   1   2
!    1   2   2
!    2   2   2
!    3   2   2
!    4   2   2
!    1   1   3
!    2   1   3
!    3   1   3
!    4   1   3
!    1   2   3
!    2   2   3
!    3   2   3
!    4   2   3
!
!  Modified:
!
!    21 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in IARRAY.  The rank of
!    the object being indexed.
!
!    Input, integer HI(N), the upper limits for the array indices.
!    The lower limit is implicitly 1, and each HI(I) should be at least 1.
!
!    Input/output, integer IARRAY(N).
!    On startup calls, with MORE = .FALSE., the input value of IARRAY
!    doesn't matter, because the routine initializes it.
!    On calls with MORE = .TRUE., the input value of IARRAY must be
!    the output value of IARRAY from the previous call.  (In other words,
!    just leave it alone!).
!    On output, IARRAY contains the successor set of indices to the input
!    value.
!
!    Input/output, logical MORE.  Set this variable .FALSE. before
!    the first call.  Normally, MORE will be returned .TRUE. but
!    once all the vectors have been generated, MORE will be
!    reset .FALSE. and you should stop calling the program.
!
  implicit none
!
  integer n
!
  integer hi(n)
  integer i
  integer iarray(n)
  integer inc
  logical more
!
  if ( .not. more ) then

    iarray(1:n) = 1

    do i = 1, n
      if ( hi(i) < 1 ) then
        more = .false.
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'INDEX_NEXT - Fatal error!'
        write ( *, '(a,i6)' ) '  Entry ', i
        write ( *, '(a,i6)' ) '  of HI is ', hi(i)
        write ( *, '(a)' ) '  but all entries must be at least 1.'
        stop
      end if
    end do

  else

    inc = 0

    do

      inc = inc + 1

      if ( iarray(inc) < hi(inc) ) then
        exit
      end if

      iarray(inc) = 1

    end do

    iarray(inc) = iarray(inc) + 1

  end if
!
!  See if there are more entries to compute.
!
  more = .false.

  do i = 1, n
    if ( iarray(i) < hi(i) ) then
      more = .true.
    end if
  end do

  return
end
subroutine pack_get ( rank, maxsize, a, dim, indx, value )
!
!*******************************************************************************
!
!! PACK_GET gets the value of an element of a packed array.
!
!
!  Modified:
!
!    22 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANK, the "rank" of the data in the packed array A.  
!    The rank is the number of indices used to access any particular entry.  
!    A scalar has rank 0, a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZE, the amount of storage set aside for the packed
!    array.
!
!    Input, real A(MAXSIZE), the packed array.
!
!    Input, integer DIM(RANK), the dimensions of the input packed array A.
!
!    Input, integer INDX(RANK), the index of the entry.
!
!    Output, real VALUE, the value of the indicated entry of A.
!
  implicit none
!
  integer maxsize
  integer rank
!
  real a(maxsize)
  integer dim(rank)
  integer i
  integer indexa
  integer indx(rank)
  real value
!
  if ( rank < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PACK_GET - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if
!
!  Check the dimensions.
!
  do i = 1, rank
    if ( indx(i) < 1 .or. indx(i) > dim(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_GET - Fatal error!'
      write ( *, '(a,i6)' ) '  INDX(I) is illegal for I = ', i
      stop
    end if
  end do
!
!  Compute the corresponding single index.
!
  indexa = 0
  do i = rank, 1, -1
    indexa = indexa * dim(i) + ( indx(i) - 1 ) 
  end do
  indexa = indexa + 1

  value = a(indexa)

  return
end
subroutine pack_inc ( rank, maxsize, a, dim, indx, value )
!
!*******************************************************************************
!
!! PACK_INC increments the value of an element of a packed array.
!
!
!  Modified:
!
!    11 January 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANK, the "rank" of the data in the packed array A.  
!    The rank is the number of indices used to access any particular entry.  
!    A scalar has rank 0, a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZE, the amount of storage set aside for the packed
!    array.
!
!    Input/output, real A(MAXSIZE), the packed array.  On output, the indicated
!    entry has been incremented by the specified value.
!
!    Input, integer DIM(RANK), the dimensions of the input packed array A.
!
!    Input, integer INDX(RANK), the index of the entry.
!
!    Input, real VALUE, the value to be added to the indicated entry of A.
!
  implicit none
!
  integer maxsize
  integer rank
!
  real a(maxsize)
  integer dim(rank)
  integer i
  integer indexa
  integer indx(rank)
  real value
!
  if ( rank < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PACK_INC - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if
!
!  Check the dimensions.
!
  do i = 1, rank
    if ( indx(i) < 1 .or. indx(i) > dim(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_INC - Fatal error!'
      write ( *, '(a,i6)' ) '  INDX(I) is illegal for I = ', i
      stop
    end if
  end do
!
!  Compute the corresponding single index.
!
  indexa = 0
  do i = rank, 1, -1
    indexa = indexa * dim(i) + ( indx(i) - 1 ) 
  end do
  indexa = indexa + 1

  a(indexa) = a(indexa) + value

  return
end
subroutine pack_select ( ranka, maxsizea, a, dima, dimlo, dimhi, rankb, &
  maxsizeb, b, dimb )
!
!*******************************************************************************
!
!! PACK_SELECT selects a subarray from a packed array.
!
!
!  Discussion:
!
!    If we write the packed array, symbolically, as
!      A( I1:J1, I2:J2, I3:J3, ..., IN:JN )
!    then the subarray can be written as:
!      B( K1:L1, K2:L2, K3:L3, ..., KN:LN ). 
!
!    Here, 1 = I1 <= K1 <= L1 <= J1, with the same restriction for
!    the other indices.  The most common cases are:
!
!      I1 = K1 and J1 = L1 
!    and
!      I1 <= K1 = L1 <= J1.
!
!    Thus, we can easily select a row or column, a sub-block, a single
!    entry, and so on.
!
!  Modified:
!
!    11 January 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANKA, the "rank" of the data in the packed array A.  
!    The rank is the number of indices used to access any particular entry.  
!    A scalar has rank 0, a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZEA, the amount of storage set aside for the packed
!    array.
!
!    Input, real A(MAXSIZEA), the packed array.
!
!    Input, integer DIMA(RANKA), the dimensions of the input packed array A.
!
!    Input, integer DIMLO(RANKA), DIMHI(RANKA), the lower and upper limits
!    on the indices of the data to be selected.  It should be the case
!    that 1 <= DIMLO(I) <= DIMHI(I) <= DIMA(I) for each I.
!
!    Output, integer RANKB, the rank of the selected array.  RANKB can
!    have any value between 0 and RANKA, and is the count of the number
!    of indices I for which DIMHI(I) - DIMLO(I) is greater than 1.
!
!    Input, integer MAXSIZEB, the amount of storage set aside for B.
!
!    Output, real B(MAXSIZEB), the packed array that contains a selection
!    of the values in A.
!
!    Output, integer DIMB(RANKA).  On output, entries 1 through RANKB of DIMB
!    contain the dimensions of the packed array B.  (DIMB must have formal
!    rank RANKA since we don't know beforehand how much space to allocate.)
!
  implicit none
!
  integer maxsizea
  integer maxsizeb
  integer ranka
!
  real a(maxsizea)
  real b(maxsizeb)
  integer, allocatable, dimension ( : ) :: digit
  integer dima(ranka)
  integer dimb(ranka)
  integer dimhi(ranka)
  integer dimlo(ranka)
  integer i
  integer iget
  integer iput
  logical inside
  logical more
  integer rankb
!
  if ( ranka < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PACK_SELECT - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if
!
!  Check the dimensions.
!
  do i = 1, ranka
    if ( dimlo(i) < 1 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMLO(I) < 1 for I = ', i
      stop
    else if ( dimlo(i) > dimhi(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMLO(I) > DIMHI(I) for I = ', i
      stop
    else if ( dimhi(i) > dima(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMHI(I) < DIMA(I) for I = ', i
      stop
    end if
  end do
!
!  Simply go through the indices of A in order.  Every time an A
!  index lies within the B range, copy the value into B.
!
  allocate ( digit(1:ranka) )

  more = .false.
  iget = 0
  iput = 0

  do

    call index_next ( ranka, dima, digit, more )

    iget = iget + 1

    inside = .true.
    do i = 1, ranka
      if ( digit(i) < dimlo(i) .or. digit(i) > dimhi(i) ) then
        inside = .false.
      end if
    end do

    if ( inside ) then
      iput = iput + 1
      b(iput) = a(iget)
    end if

    if ( .not. more ) then
      exit
    end if

  end do
!
!  Now compute the rank and dimension of B.
!
  rankb = 0

  do i = 1, ranka
    if ( dimlo(i) < dimhi(i) ) then
      rankb = rankb + 1
      dimb(rankb) = dimhi(i) + 1 - dimlo(i)
    end if
  end do

  deallocate ( digit )

  return
end
subroutine pack_set ( rank, maxsize, a, dim, indx, value )
!
!*******************************************************************************
!
!! PACK_SET sets the value of an element of a packed array.
!
!
!  Modified:
!
!    11 January 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANK, the "rank" of the data in the packed array A.  
!    The rank is the number of indices used to access any particular entry.  
!    A scalar has rank 0, a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZE, the amount of storage set aside for the packed
!    array.
!
!    Input/output, real A(MAXSIZE), the packed array.  On output, the indicated
!    entry has been set to the specified value.
!
!    Input, integer DIM(RANK), the dimensions of the input packed array A.
!
!    Input, integer INDX(RANK), the index of the entry.
!
!    Input, real VALUE, the value to be assigned the indicated entry of A.
!
  implicit none
!
  integer maxsize
  integer rank
!
  real a(maxsize)
  integer dim(rank)
  integer i
  integer indexa
  integer indx(rank)
  real value
!
  if ( rank < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PACK_SET - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if
!
!  Check the dimensions.
!
  do i = 1, rank
    if ( indx(i) < 1 .or. indx(i) > dim(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'PACK_SET - Fatal error!'
      write ( *, '(a,i6)' ) '  INDX(I) is illegal for I = ', i
      stop
    end if
  end do
!
!  Compute the corresponding single index.
!
  indexa = 0
  do i = rank, 1, -1
    indexa = indexa * dim(i) + ( indx(i) - 1 ) 
  end do
  indexa = indexa + 1

  a(indexa) = value

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end
subroutine unpack_select ( ranka, maxsizea, a, maxdima, dima, dimlo, dimhi, &
  rankb, maxsizeb, b, dimb )
!
!*******************************************************************************
!
!! UNPACK_SELECT selects a subarray from an unpacked array.
!
!
!  Discussion:
!
!    If we write the unpacked array, symbolically, as
!      A( I1:J1, I2:J2, I3:J3, ..., IN:JN )
!    then the subarray can be written as:
!      B( K1:L1, K2:L2, K3:L3, ..., KN:LN ). 
!
!    Here, I1 <= K1 <= L1 <= J1, with the same restriction for
!    the other indices.  The most common cases are:
!
!      I1 = K1 and J1 = L1 
!    and
!      I1 <= K1 = L1 <= J1.
!
!  Modified:
!
!    21 July 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANKA, the "rank" of the data in the unpacked array A.  
!    The rank is the number of indices used to access any particular entry.  
!    A scalar has rank 0, a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZEA, the amount of storage set aside for the unpacked
!    array.
!
!    Input, real A(MAXSIZEA), the unpacked array.
!
!    Input, integer MAXDIMA(RANKA), the maximum dimensions of the unpacked 
!    array A.  
!
!    Input, integer DIMA(RANKA), the used dimensions of the input unpacked 
!    array A.
!
!    Input, integer DIMLO(RANKA), DIMHI(RANKA), the lower and upper limits
!    on the indices of the data to be selected.  It should be the case
!    that 1 <= DIMLO(I) <= DIMHI(I) <= DIMA(I) for each I.
!
!    Output, integer RANKB, the rank of the selected array.  RANKB can
!    have any value between 0 and RANKA, and is the count of the number
!    of indices I for which DIMHI(I) - DIMLO(I) is greater than 1.
!
!    Input, integer MAXSIZEB, the amount of storage set aside for B.
!
!    Output, real B(MAXSIZEB), the packed array that contains a selection
!    of the values in A.
!
!    Output, integer DIMB(RANKA).  On output, entries 1 through RANKB of DIMB
!    contain the dimensions of the packed array B.  (DIMB must have formal
!    rank RANKA since we don't know beforehand how much space to allocate.)
!
  implicit none
!
  integer maxsizea
  integer maxsizeb
  integer ranka
!
  real a(maxsizea)
  real b(maxsizeb)
  integer, allocatable, dimension ( : ) :: digit
  integer dima(ranka)
  integer dimb(ranka)
  integer dimhi(ranka)
  integer dimlo(ranka)
  integer i
  integer iget
  integer iput
  logical inside
  integer maxdima(ranka)
  logical more
  integer rankb
!
  if ( ranka < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'UNPACK_SELECT - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if
!
!  Check the dimensions.
!
  do i = 1, ranka
    if ( dimlo(i) < 1 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'UNPACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMLO(I) < 1 for I = ', i
      stop
    else if ( dimlo(i) > dimhi(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'UNPACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMLO(I) > DIMHI(I) for I = ', i
      stop
    else if ( dimhi(i) > dima(i) ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'UNPACK_SELECT - Fatal error!'
      write ( *, '(a,i6)' ) '  DIMHI(I) < DIMA(I) for I = ', i
      stop
    end if
  end do
!
!  Simply go through the indices of A in order.  Every time an A
!  index lies within the B range, copy the value into B.
!

  allocate ( digit(1:ranka) )

  more = .false.
  iget = 0
  iput = 0

  do

    call index_next ( ranka, maxdima, digit, more )

    iget = iget + 1

    inside = .true.
    do i = 1, ranka
      if ( digit(i) < dimlo(i) .or. digit(i) > dimhi(i) ) then
        inside = .false.
      end if
    end do

    if ( inside ) then
      iput = iput + 1
      b(iput) = a(iget)
    end if

    if ( .not. more ) then
      exit
    end if

  end do
!
!  Now compute the rank and dimension of B.
!
  rankb = 0

  do i = 1, ranka
    if ( dimlo(i) < dimhi(i) ) then
      rankb = rankb + 1
      dimb(rankb) = dimhi(i) + 1 - dimlo(i)
    end if
  end do

  deallocate ( digit )

  return
end
subroutine unpack_to_pack ( rank, maxsize, a, maxdim, dim )
!
!*******************************************************************************
!
!! UNPACK_TO_PACK converts an unpacked array to a packed array.
!
!
!  Discussion:
!
!    The difference between the unpacked and packed arrays is that the
!    unpacked array has maximum dimensions that are larger than the used 
!    dimensions.  This means that there are various unused "placeholder" 
!    entries in the stored version of the unpacked array which are not
!    being used to store data.  The packed array is a copy of the data in 
!    the unpacked array, with the unused entries "squeezed out".
!
!  Example:
!
!    Here is how the routine would handle a 2D unpacked array which contains
!    3 by 3 data inside a 5 by 4 space:
!
!    Input:
!
!      A = ( 11, 21, 31, 00, 00, 12, 22, 32, 00, 00, 13, 23, 33, 00, 00,
!            00, 00, 00, 00, 00 )
!      MAXSIZE = 20
!      RANK = 2
!      MAXDIM = ( 5, 4 )
!      DIM = ( 3, 3 )
!
!    Output:
!
!      A = ( 11, 21, 31, 12, 22, 32, 13, 23, 33, ...garbage... )
!
!  Modified:
!
!    11 January 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer RANK, the "rank" of the data.  The rank is the number of
!    indices used to access any particular entry.  A scalar has rank 0,
!    a vector rank 1, a matrix rank 2, and so on.
!
!    Input, integer MAXSIZE, the amount of storage set aside for the array.
!
!    Input/output, real A(MAXSIZE), the set of consecutive storage locations
!    in which the array is stored. 
!    On input, the A data is stored as an unpacked array.
!    On output, the A data is stored as a packed array.
!
!    Input, integer MAXDIM(RANK), the maximum dimensions of the unpacked 
!    array A.  
!
!    Input, integer DIM(RANK), the "used" dimensions of the input unpacked 
!    array A, and the dimensions of the output packed array A.
!
  implicit none
!
  integer maxsize
  integer rank
!
  real a(maxsize)
  integer, allocatable, dimension ( : ) :: digit
  integer dim(rank)
  integer i
  integer idig
  integer iget
  integer iput
  integer maxdim(rank)
!
  if ( rank < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'UNPACK_TO_PACK - Fatal error!'
    write ( *, '(a)' ) '  Object rank less than 1.'
    stop
  end if

  allocate ( digit(1:rank) )

  iget = 0
  iput = 0

  digit(1:rank) = 0

  do

    idig = 1
    digit(idig) = digit(idig) + 1

2   continue

    if ( digit(idig) <= dim(idig) ) then
      iget = iget + 1
      iput = iput + 1
      a(iput) = a(iget)
    else if ( digit(idig) <= maxdim(idig) ) then
      iget = iget + 1
    else if ( idig < rank ) then
      digit(idig) = 1
      idig = idig + 1
      digit(idig) = digit(idig) + 1
      go to 2
    else
      exit
    end if
  
  end do

  deallocate ( digit )

  return
end