subroutine sparse_grid_mixed_index ( dim_num, level_max, rule, point_num, &
  point_total_num, sparse_unique_index, sparse_order, sparse_index )

!*****************************************************************************80
!
!! sparse_grid_mixed_index() indexes a sparse grid made from mixed 1D rules.
!
!  Discussion:
!
!    For each "unique" point in the sparse grid, we return its INDEX and ORDER.
!
!    That is, for the I-th unique point P, we determine the product grid which
!    first generated this point, and we return in SPARSE_ORDER the orders of
!    the 1D rules in that grid, and in SPARSE_INDEX the component indexes in
!    those rules that generated this specific point.
!
!    For instance, say P was first generated by a rule which was a 3D product
!    of a 9th order CC rule and a 15th order GL rule, and that to generate P,
!    we used the 7-th point of the CC rule and the 3rd point of the GL rule.
!    Then the SPARSE_ORDER information would be (9,15) and the SPARSE_INDEX
!    information would be (7,3).  This, combined with the information in RULE,
!    is enough to regenerate the value of P.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 December 2009
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, the maximum value of LEVEL.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, integer POINT_NUM, the number of unique points
!    in the grid.
!
!    Input, integer POINT_TOTAL_NUM, the total number of points
!    in the grid.
!
!    Input, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM),
!    associates each point in the grid with its unique representative.
!
!    Output, integer SPARSE_ORDER(DIM_NUM,POINT_NUM), lists,
!    for each point, the order of the 1D rules used in the grid that
!    generated it.
!
!    Output, integer SPARSE_INDEX(DIM_NUM,POINT_NUM), lists, for
!    each point, its index in each of the 1D rules in the grid that generated
!    it.  The indices are 1-based.
!
  implicit none

  integer dim_num
  integer point_num
  integer point_total_num

  integer h
  integer level
  integer level_1d(dim_num)
  integer level_max
  integer level_min
  logical more_grids
  logical more_points
  integer order_1d(dim_num)
  integer point_count
  integer point_index(dim_num)
  integer point_unique
  integer rule(dim_num)
  integer sparse_index(dim_num,point_num)
  integer sparse_order(dim_num,point_num)
  integer sparse_unique_index(point_total_num)
  integer t
!
!  Special cases.
!
  if ( level_max < 0 ) then
    return
  end if

  if ( level_max == 0 ) then
    sparse_order(1:dim_num,1) = 1
    sparse_index(1:dim_num,1) = 1
    return
  end if
!
!  Initialize to -1 to help catch errors.
!
  sparse_order(1:dim_num,1:point_num) = -1
  sparse_index(1:dim_num,1:point_num) = -1

  point_count = 0
!
!  The outer loop generates values of LEVEL.
!
  level_min = max ( 0, level_max + 1 - dim_num )

  do level = level_min, level_max
!
!  The middle loop generates a GRID,
!  based on the next partition that adds up to LEVEL.
!
    more_grids = .false.
    h = 0
    t = 0

    do

      call comp_next ( level, dim_num, level_1d, more_grids, h, t )

      call level_to_order_default ( dim_num, level_1d, rule, order_1d )
!
!  The inner loop generates a POINT of the GRID of the LEVEL.
!
      more_points = .false.

      do

        call vec_colex_next3 ( dim_num, order_1d, point_index, more_points )

        if ( .not. more_points ) then
          exit
        end if

        point_count = point_count + 1
        point_unique = sparse_unique_index(point_count)
        sparse_order(1:dim_num,point_unique) = order_1d(1:dim_num)
        sparse_index(1:dim_num,point_unique) = point_index(1:dim_num)

      end do

      if ( .not. more_grids ) then
        exit
      end if

    end do

  end do

  return
end
subroutine sparse_grid_mixed_point ( dim_num, level_max, rule, alpha, beta, &
  point_num, sparse_order, sparse_index, sparse_point )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_POINT computes the points of a sparse grid rule.
!
!  Discussion:
!
!    The sparse grid is the logical sum of low degree product rules.
!
!    Each product rule is the product of 1D factor rules.
!
!    The user specifies:
!    * the spatial dimension of the quadrature region,
!    * the level that defines the Smolyak grid.
!    * the quadrature rules.
!    * the number of points.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 March 2011
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, controls the size of the final
!    sparse grid.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, real ( kind = rk ) ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
!    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss
!    Jacobi rules.
!
!    Input, integer POINT_NUM, the number of unique points
!    in the grid.
!
!    Input, integer SPARSE_ORDER(DIM_NUM,POINT_NUM), lists, for
!    each point, the order of the 1D rules used in the grid that generated it.
!
!    Input, integer SPARSE_INDEX(DIM_NUM,POINT_NUM), lists, for
!    each point, its index in each of the 1D rules in the grid that generated
!    it.  The indices are 1-based.
!
!    Output, real ( kind = rk ) SPARSE_POINT(DIM_NUM,POINT_NUM), the points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer dim_num
  integer point_num

  real ( kind = rk ) alpha(dim_num)
  real ( kind = rk ) beta(dim_num)
  integer dim
  integer level
  integer level_max
  integer order
  integer point
  real ( kind = rk ), allocatable, dimension ( : ) :: points
  real ( kind = rk ) r8_huge
  integer rule(dim_num)
  integer sparse_index(dim_num,point_num)
  integer sparse_order(dim_num,point_num)
  real ( kind = rk ) sparse_point(dim_num,point_num)
!
!  Compute the point coordinates.
!
  sparse_point(1:dim_num,1:point_num) = - r8_huge ( )

  do dim = 1, dim_num

    do level = 0, level_max

      call level_to_order_default ( 1, level, rule(dim), order )

      allocate ( points(1:order) )

      if ( rule(dim) == 1 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 2 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 3 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 4 ) then
        call legendre_compute_points ( order, points )
      else if ( rule(dim) == 5 ) then
        call hermite_compute_points ( order, points )
      else if ( rule(dim) == 6 ) then
        call gen_hermite_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 7 ) then
        call laguerre_compute_points ( order, points )
      else if ( rule(dim) == 8 ) then
        call gen_laguerre_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 9 ) then
        call jacobi_compute_points ( order, alpha(dim), beta(dim), points )
      else if ( rule(dim) == 10 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_POINT - Fatal error!'
        write ( *, '(a)' ) '  Do not know how to assign points for rule 10.'
        stop
      else if ( rule(dim) == 11 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 12 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 13 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 14 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 15 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 16 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 17 ) then
        call ccn_compute_points ( order, points )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_POINT - Fatal error!'
        write ( *, '(a,i8)' ) '  Unexpected value of RULE = ', rule(dim)
        stop
      end if

      do point = 1, point_num
        if ( sparse_order(dim,point) == order ) then
          sparse_point(dim,point) = points ( sparse_index(dim,point) )
        end if
      end do

      deallocate ( points )

    end do
  end do

  return
end
subroutine sparse_grid_mixed_size ( dim_num, level_max, rule, alpha, beta, &
  tol, point_num )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_SIZE sizes a sparse grid, discounting duplicate points.
!
!  Discussion:
!
!    The sparse grid is the logical sum of product grids with total LEVEL
!    between LEVEL_MIN and LEVEL_MAX.
!
!    Depending on the 1D rules involved, there may be many duplicate points
!    in the sparse grid.
!
!    This routine counts the unique points in the sparse grid.  It does this
!    in a straightforward way, by actually generating all the points, and
!    comparing them, with a tolerance for equality.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 March 2011
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, the maximum value of LEVEL.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, real ( kind = rk ) ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
!    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss
!    Jacobi rules.
!
!    Input, real ( kind = rk ) TOL, the tolerance for point equality.
!
!    Output, integer POINT_NUM, the number of unique points
!    in the grid.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer dim_num

  real ( kind = rk ) alpha(dim_num)
  real ( kind = rk ) beta(dim_num)
  integer dim
  integer h
  integer level
  integer level_1d(dim_num)
  integer level_max
  integer level_min
  logical more_grids
  logical more_points
  integer order
  integer order_1d(dim_num)
  integer point
  integer point_index(dim_num)
  integer point_num
  integer point_total_num
  integer point_total_num2
  real ( kind = rk ), allocatable, dimension ( : ) :: points
  real ( kind = rk ) r8_huge
  integer rule(dim_num)
  integer, allocatable, dimension ( :, : ) :: sparse_total_index
  integer, allocatable, dimension ( :, : ) :: sparse_total_order
  real ( kind = rk ), allocatable, dimension ( :, : ) :: sparse_total_point
  integer t
  real ( kind = rk ) tol
!
!  Special cases.
!
  if ( level_max < 0 ) then
    point_num = -1
    return
  end if

  if ( level_max == 0 ) then
    point_num = 1
    return
  end if
!
!  Get total number of points, including duplicates.
!
  call sparse_grid_mixed_size_total ( dim_num, level_max, rule, &
    point_total_num )
!
!  Generate SPARSE_TOTAL_ORDER and SPARSE_TOTAL_INDEX arrays
!  for the TOTAL set of points.
!
  allocate ( sparse_total_order(1:dim_num,1:point_total_num ) )
  allocate ( sparse_total_index(1:dim_num,1:point_total_num ) )

  point_total_num2 = 0
!
!  The outer loop generates values of LEVEL.
!
  level_min = max ( 0, level_max + 1 - dim_num )

  do level = level_min, level_max
!
!  The middle loop generates a GRID,
!  based on the next partition that adds up to LEVEL.
!
    more_grids = .false.
    h = 0
    t = 0

    do

      call comp_next ( level, dim_num, level_1d, more_grids, h, t )

      call level_to_order_default ( dim_num, level_1d, rule, order_1d )
!
!  The inner loop generates a POINT of the GRID of the LEVEL.
!
      more_points = .false.

      do

        call vec_colex_next3 ( dim_num, order_1d, point_index, more_points )

        if ( .not. more_points ) then
          exit
        end if

        point_total_num2 = point_total_num2 + 1
        sparse_total_order(1:dim_num,point_total_num2) = order_1d(1:dim_num)
        sparse_total_index(1:dim_num,point_total_num2) = point_index(1:dim_num)

      end do

      if ( .not. more_grids ) then
        exit
      end if

    end do

  end do
!
!  Now compute the coordinates of the TOTAL set of points.
!
  allocate ( sparse_total_point(1:dim_num,1:point_total_num) )
  sparse_total_point(1:dim_num,1:point_total_num) = r8_huge ( )

  do dim = 1, dim_num

    do level = 0, level_max

      call level_to_order_default ( 1, level, rule(dim), order )

      allocate ( points(1:order) )

      if ( rule(dim) == 1 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 2 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 3 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 4 ) then
        call legendre_compute_points ( order, points )
      else if ( rule(dim) == 5 ) then
        call hermite_compute_points ( order, points )
      else if ( rule(dim) == 6 ) then
        call gen_hermite_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 7 ) then
        call laguerre_compute_points ( order, points )
      else if ( rule(dim) == 8 ) then
        call gen_laguerre_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 9 ) then
        call jacobi_compute_points ( order, alpha(dim), beta(dim), points )
      else if ( rule(dim) == 10 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_SIZE - Fatal error!'
        write ( *, '(a)' ) '  Do not know how to assign points for rule 10.'
        stop
      else if ( rule(dim) == 11 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 12 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 13 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 14 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 15 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 16 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 17 ) then
        call ccn_compute_points ( order, points )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_SIZE - Fatal error!'
        write ( *, '(a,i8)' ) '  Unexpected value of RULE = ', rule(dim)
        stop
      end if

      do point = 1, point_total_num
        if ( sparse_total_order(dim,point) == order ) then
          sparse_total_point(dim,point) = &
            points ( sparse_total_index(dim,point) )
        end if
      end do

      deallocate ( points )

    end do

  end do
!
!  Sort the columns.
!
! call r8col_sort_heap_a ( dim_num, point_total_num, sparse_total_point )
!
!  Count the unique columns.
!
  call r8col_sorted_unique_count ( dim_num, point_total_num, &
    sparse_total_point, tol, point_num )

! call r8col_tol_unique_count ( dim_num, point_total_num, &
!   sparse_total_point, tol, point_num )

  deallocate ( sparse_total_index )
  deallocate ( sparse_total_order )
  deallocate ( sparse_total_point )

  return
end
subroutine sparse_grid_mixed_size_total ( dim_num, level_max, rule, &
  point_total_num )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_SIZE_TOTAL sizes a sparse grid, counting duplicate points.
!
!  Discussion:
!
!    The sparse grid is the logical sum of product grids with total LEVEL
!    between LEVEL_MIN and LEVEL_MAX.
!
!    In some cases, the same point may occur in different product grids
!    used to form the sparse grid.
!
!    This routine counts the total number of points used to construct the sparse
!    grid; if the same point occurs several times, each occurrence is added
!    to the sum.
!
!    This computation is useful in order to be able to allocate enough
!    space for the full set of points, before they are compressed by removing
!    duplicates.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 December 2009
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, the maximum value of LEVEL.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Output, integer POINT_TOTAL_NUM, the total number of points
!    in the grid.
!
  implicit none

  integer dim_num

  integer h
  integer level
  integer level_1d(dim_num)
  integer level_max
  integer level_min
  logical more_grids
  integer order_1d(dim_num)
  integer point_total_num
  integer rule(dim_num)
  integer t
!
!  Special case.
!
  if ( level_max == 0 ) then
    point_total_num = 1
    return
  end if

  point_total_num = 0
!
!  The outer loop generates values of LEVEL.
!
  level_min = max ( 0, level_max + 1 - dim_num )

  do level = level_min, level_max
!
!  The middle loop generates a GRID,
!  based on the next partition that adds up to LEVEL.
!
    more_grids = .false.
    h = 0
    t = 0

    do

      call comp_next ( level, dim_num, level_1d, more_grids, h, t )

      call level_to_order_default ( dim_num, level_1d, rule, order_1d )

      point_total_num = point_total_num + product ( order_1d(1:dim_num) )

      if ( .not. more_grids ) then
        exit
      end if

    end do

  end do

  return
end
subroutine sparse_grid_mixed_unique_index ( dim_num, level_max, rule, &
  alpha, beta, tol, point_num, point_total_num, sparse_unique_index )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_UNIQUE_INDEX maps nonunique points to unique points.
!
!  Discussion:
!
!    The sparse grid usually contains many points that occur in more
!    than one product grid.
!
!    When generating the point locations, it is easy to realize that a point
!    has already been generated.
!
!    But when it's time to compute the weights of the sparse grids, it is
!    necessary to handle situations in which weights corresponding to
!    the same point generated in multiple grids must be collected together.
!
!    This routine generates ALL the points, including their multiplicities,
!    and figures out a mapping from them to the collapsed set of unique points.
!
!    This mapping can then be used during the weight calculation so that
!    a contribution to the weight gets to the right place.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 July 2010
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, the maximum value of LEVEL.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, real ( kind = rk ) ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
!    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss
!    Jacobi rules.
!
!    Input, real ( kind = rk ) TOL, the tolerance for point equality.
!
!    Input, integer POINT_NUM, the number of unique points in
!    the grid.
!
!    Input, integer POINT_TOTAL_NUM, the total number of points
!    in the grid.
!
!    Output, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM), lists,
!    for each (nonunique) point, the corresponding index of the same point in
!    the unique listing.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer dim_num
  integer point_total_num

  real ( kind = rk ) alpha(dim_num)
  real ( kind = rk ) beta(dim_num)
  integer dim
  integer h
  integer level
  integer level_1d(dim_num)
  integer level_max
  integer level_min
  logical more_grids
  logical more_points
  integer order
  integer order_1d(dim_num)
  integer point
  integer point_index(dim_num)
  integer point_num
  integer point_total_num2
  real ( kind = rk ), allocatable, dimension ( : ) :: points
  real ( kind = rk ) r8_huge
  integer rule(dim_num)
  integer, allocatable, dimension ( :, : ) :: sparse_total_index
  integer, allocatable, dimension ( :, : ) :: sparse_total_order
  real ( kind = rk ), allocatable, dimension ( :, : ) :: sparse_total_point
  integer sparse_unique_index(point_total_num)
  integer t
  real ( kind = rk ) tol
  integer, allocatable, dimension ( : ) :: undx
!
!  Special cases.
!
  if ( level_max < 0 ) then
    return
  end if

  if ( level_max == 0 ) then
    sparse_unique_index(1) = 1
    return
  end if
!
!  Generate SPARSE_TOTAL_ORDER and SPARSE_TOTAL_INDEX arrays
!  for the TOTAL set of points.
!
  allocate ( sparse_total_order(1:dim_num,1:point_total_num ) )
  allocate ( sparse_total_index(1:dim_num,1:point_total_num ) )

  point_total_num2 = 0
!
!  The outer loop generates values of LEVEL.
!
  level_min = max ( 0, level_max + 1 - dim_num )

  do level = level_min, level_max
!
!  The middle loop generates a GRID,
!  based on the next partition that adds up to LEVEL.
!
    more_grids = .false.
    h = 0
    t = 0

    do

      call comp_next ( level, dim_num, level_1d, more_grids, h, t )

      call level_to_order_default ( dim_num, level_1d, rule, order_1d )
!
!  The inner loop generates a POINT of the GRID of the LEVEL.
!
      more_points = .false.

      do

        call vec_colex_next3 ( dim_num, order_1d, point_index, more_points )

        if ( .not. more_points ) then
          exit
        end if

        point_total_num2 = point_total_num2 + 1
        sparse_total_order(1:dim_num,point_total_num2) = order_1d(1:dim_num)
        sparse_total_index(1:dim_num,point_total_num2) = point_index(1:dim_num)

      end do

      if ( .not. more_grids ) then
        exit
      end if

    end do

  end do
!
!  Now compute the coordinates of the TOTAL set of points.
!
  allocate ( sparse_total_point(1:dim_num,1:point_total_num) )
  sparse_total_point(1:dim_num,1:point_total_num) = r8_huge ( )

  do dim = 1, dim_num

    do level = 0, level_max

      call level_to_order_default ( 1, level, rule(dim), order )

      allocate ( points(1:order) )

      if ( rule(dim) == 1 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 2 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 3 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 4 ) then
        call legendre_compute_points ( order, points )
      else if ( rule(dim) == 5 ) then
        call hermite_compute_points ( order, points )
      else if ( rule(dim) == 6 ) then
        call gen_hermite_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 7 ) then
        call laguerre_compute_points ( order, points )
      else if ( rule(dim) == 8 ) then
        call gen_laguerre_compute_points ( order, alpha(dim), points )
      else if ( rule(dim) == 9 ) then
        call jacobi_compute_points ( order, alpha(dim), beta(dim), points )
      else if ( rule(dim) == 10 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_UNIQUE_INDEX - Fatal error!'
        write ( *, '(a)' ) '  Do not know how to assign points for rule 10.'
        stop
      else if ( rule(dim) == 11 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 12 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 13 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 14 ) then
        call clenshaw_curtis_compute_points ( order, points )
      else if ( rule(dim) == 15 ) then
        call fejer2_compute_points ( order, points )
      else if ( rule(dim) == 16 ) then
        call patterson_lookup_points ( order, points )
      else if ( rule(dim) == 17 ) then
        call ccn_compute_points ( order, points )
      else
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'SPARSE_GRID_MIXED_UNIQUE_INDEX - Fatal error!'
        write ( *, '(a,i8)' ) '  Unexpected value of RULE = ', rule(dim)
        stop
      end if

      do point = 1, point_total_num
        if ( sparse_total_order(dim,point) == order ) then
          sparse_total_point(dim,point) = &
            points ( sparse_total_index(dim,point) )
        end if
      end do

      deallocate ( points )

    end do

  end do
!
!  Now determine the mapping from nonunique points to unique points.
!  We can't really use the UNDX output right now.
!
  allocate ( undx(1:point_num) )

  call r8col_undex ( dim_num, point_total_num, sparse_total_point, point_num, &
    tol, undx, sparse_unique_index )

! call r8col_tol_undex ( dim_num, point_total_num, sparse_total_point, &
!   point_num, tol, undx, sparse_unique_index )

  deallocate ( sparse_total_index )
  deallocate ( sparse_total_order )
  deallocate ( sparse_total_point )
  deallocate ( undx )

  return
end
subroutine sparse_grid_mixed_weight ( dim_num, level_max, rule, alpha, beta, &
  point_num, point_total_num, sparse_unique_index, sparse_weight )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_WEIGHT computes sparse grid weights for a mix of 1D rules.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 December 2009
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer LEVEL_MAX, the maximum value of LEVEL.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, real ( kind = rk ) ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
!    Generalized Gauss Hermite, Generalized Gauss Laguerre, and
!    Gauss Jacobi rules.
!
!    Input, integer POINT_NUM, the number of unique points in
!    the grid.
!
!    Input, integer POINT_TOTAL_NUM, the total number of points
!    in the grid.
!
!    Input, integer SPARSE_UNIQUE_INDEX(POINT_TOTAL_NUM), lists,
!    for each (nonunique) point, the corresponding index of the same point in
!    the unique listing.
!
!    Output, real ( kind = rk ) SPARSE_WEIGHT(POINT_NUM), the weights
!    associated with the sparse grid points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer dim_num
  integer point_num
  integer point_total_num

  real ( kind = rk ) alpha(dim_num)
  real ( kind = rk ) beta(dim_num)
  real ( kind = rk ) coeff
  real ( kind = rk ), allocatable, dimension ( : ) :: grid_weight
  integer h
  integer level
  integer level_1d(dim_num)
  integer level_max
  integer level_min
  logical more_grids
  integer order
  integer order_1d(dim_num)
  integer order_nd
  integer point_total
  integer point_unique
  real ( kind = rk ) r8_choose
  real ( kind = rk ) r8_mop
  integer rule(dim_num)
  integer sparse_unique_index(point_total_num)
  real ( kind = rk ) sparse_weight(point_num)
  integer t

  sparse_weight(1:point_num) = 0.0D+00

  point_total = 0

  level_min = max ( 0, level_max + 1 - dim_num )

  do level = level_min, level_max
!
!  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
!  that adds up to LEVEL.
!
    more_grids = .false.
    h = 0
    t = 0

    do

      call comp_next ( level, dim_num, level_1d, more_grids, h, t )
!
!  Transform each 1D level to a corresponding 1D order.
!
      call level_to_order_default ( dim_num, level_1d, rule, order_1d )
!
!  The product of the 1D orders gives us the number of points in this grid.
!
      order_nd = product ( order_1d(1:dim_num) )
!
!  Compute the weights for this grid.
!
!  The correct transfer of data from the product grid to the sparse grid
!  depends on the fact that the product rule weights are stored under colex
!  order of the points, and this is the same ordering implicitly used in
!  generating the SPARSE_UNIQUE_INDEX array.
!
      allocate ( grid_weight(1:order_nd) )

      call product_mixed_weight ( dim_num, order_1d, order_nd, rule, &
        alpha, beta, grid_weight )
!
!  Compute Smolyak's binomial coefficient for this grid.
!
      coeff = r8_mop ( level_max - level ) &
        * r8_choose ( dim_num - 1, level_max - level )
!
!  Add these weights to the rule.
!
      do order = 1, order_nd

        point_total = point_total + 1

        point_unique = sparse_unique_index(point_total)

        sparse_weight(point_unique) = sparse_weight(point_unique) &
          + coeff * grid_weight(order)

      end do

      deallocate ( grid_weight )

      if ( .not. more_grids ) then
        exit
      end if

    end do

  end do

  return
end
subroutine sparse_grid_mixed_write ( dim_num, rule, alpha, beta, point_num, &
  sparse_weight, sparse_point, file_name )

!*****************************************************************************80
!
!! SPARSE_GRID_MIXED_WRITE writes a sparse grid rule to five files.
!
!  Discussion:
!
!    The files are:
!    * the "A" file stores the ALPHA values, as a DIM_NUM x 1 list.
!    * the "B" file stores the BETA values, as a DIM_NUM x 1 list.
!    * the "R" file stores the region, as a DIM_NUM x 2 list.
!    * the "W" file stores the weights as a POINT_NUM list;
!    * the "X" file stores the abscissas as a DIM_NUM x POINT_NUM list;
!
!    The entries in the "R" file are the two corners of the DIM_NUM dimensional
!    rectangle that constitutes the integration region.  Coordinates that
!    should be infinite are set to 1.0E+30.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 February 2010
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Fabio Nobile, Raul Tempone, Clayton Webster,
!    A Sparse Grid Stochastic Collocation Method for Partial Differential
!    Equations with Random Input Data,
!    SIAM Journal on Numerical Analysis,
!    Volume 46, Number 5, 2008, pages 2309-2345.
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer RULE(DIM_NUM), the rule in each dimension.
!     1, "CC",  Clenshaw Curtis, Closed Fully Nested.
!     2, "F2",  Fejer Type 2, Open Fully Nested.
!     3, "GP",  Gauss Patterson, Open Fully Nested.
!     4, "GL",  Gauss Legendre, Open Weakly Nested.
!     5, "GH",  Gauss Hermite, Open Weakly Nested.
!     6, "GGH", Generalized Gauss Hermite, Open Weakly Nested.
!     7, "LG",  Gauss Laguerre, Open Non Nested.
!     8, "GLG", Generalized Gauss Laguerre, Open Non Nested.
!     9, "GJ",  Gauss Jacobi, Open Non Nested.
!    10, "GW",  Golub Welsch, (presumed) Open Non Nested.
!    11, "CC_SE", Clenshaw Curtis Slow Exponential, Closed Fully Nested.
!    12, "F2_SE", Fejer Type 2 Slow Exponential, Closed Fully Nested.
!    13, "GP_SE", Gauss Patterson Slow Exponential, Closed Fully Nested.
!    14, "CC_ME", Clenshaw Curtis Moderate Exponential, Closed Fully Nested.
!    15, "F2_ME", Fejer Type 2 Moderate Exponential, Closed Fully Nested.
!    16, "GP_ME", Gauss Patterson Moderate Exponential, Closed Fully Nested.
!    17, "CCN", Clenshaw Curtis Nested, Linear, Closed Fully Nested rule.
!
!    Input, real ( kind = rk ) ALPHA(DIM_NUM), BETA(DIM_NUM), parameters used for
!    Generalized Gauss Hermite, Generalized Gauss Laguerre, and Gauss
!    Jacobi rules.
!
!    Input, integer POINT_NUM, the number of unique points
!    in the grid.
!
!    Input, real ( kind = rk ) SPARSE_WEIGHT(POINT_NUM), the weights.
!
!    Input, real ( kind = rk ) SPARSE_POINT(DIM_NUM,POINT_NUM), the points.
!
!    Input, character ( len = * ) FILE_NAME, the main part of the file name.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer dim_num
  integer point_num

  real ( kind = rk ) alpha(dim_num)
  real ( kind = rk ) beta(dim_num)
  integer dim
  character ( len = *  ) file_name
  character ( len = 80 ) file_name_a
  character ( len = 80 ) file_name_b
  character ( len = 80 ) file_name_r
  character ( len = 80 ) file_name_w
  character ( len = 80 ) file_name_x
  real ( kind = rk ) r8_huge
  integer rule(dim_num)
  real ( kind = rk ) sparse_point(dim_num,point_num)
  real ( kind = rk ) sparse_region(dim_num,2)
  real ( kind = rk ) sparse_weight(point_num)

  do dim = 1, dim_num
    if ( rule(dim) == 1 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 2 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 3 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 4 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 5 ) then
      sparse_region(dim,1) = - r8_huge ( )
      sparse_region(dim,2) = + r8_huge ( )
    else if ( rule(dim) == 6 ) then
      sparse_region(dim,1) = - r8_huge ( )
      sparse_region(dim,2) = + r8_huge ( )
    else if ( rule(dim) == 7 ) then
      sparse_region(dim,1) = 0.0D+00
      sparse_region(dim,2) = r8_huge ( )
    else if ( rule(dim) == 8 ) then
      sparse_region(dim,1) = 0.0D+00
      sparse_region(dim,2) = r8_huge ( )
    else if ( rule(dim) == 9 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 10 ) then
      sparse_region(dim,1) = minval ( sparse_point(dim,1:point_num) )
      sparse_region(dim,2) = maxval ( sparse_point(dim,1:point_num) )
    else if ( rule(dim) == 11 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 12 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 13 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 14 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 15 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 16 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else if ( rule(dim) == 17 ) then
      sparse_region(dim,1) = -1.0D+00
      sparse_region(dim,2) = +1.0D+00
    else
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'SPARSE_GRID_MIXED_WRITE - Fatal error!'
      write ( *, '(a,i8)' ) '  Unexpected value of RULE = ', rule(dim)
      stop
    end if
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'SPARSE_GRID_MIXED_WRITE:'

  file_name_a = trim ( file_name ) // '_a.txt'
  call r8mat_write ( file_name_a, dim_num, 1, alpha )
  write ( *, '(a)' ) '  Wrote the A file = "' // trim ( file_name_a ) // '".'

  file_name_b = trim ( file_name ) // '_b.txt'
  call r8mat_write ( file_name_b, dim_num, 1, beta )
  write ( *, '(a)' ) '  Wrote the B file = "' // trim ( file_name_b ) // '".'

  file_name_r = trim ( file_name ) // '_r.txt'
  call r8mat_write ( file_name_r, dim_num, 2, sparse_region )
  write ( *, '(a)' ) '  Wrote the R file = "' // trim ( file_name_r ) // '".'

  file_name_w = trim ( file_name ) // '_w.txt'
  call r8mat_write ( file_name_w, 1, point_num, sparse_weight )
  write ( *, '(a)' ) '  Wrote the W file = "' // trim ( file_name_w ) // '".'

  file_name_x = trim ( file_name ) // '_x.txt'
  call r8mat_write ( file_name_x, dim_num, point_num, sparse_point )
  write ( *, '(a)' ) '  Wrote the X file = "' // trim ( file_name_x ) // '".'

  return
end