function sparse_grid_gl_dataset ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_GL_DATASET is the main program.
%
%  Discussion:
%
%    This program computes a sparse grid quadrature rule based on a 1D
%    Gauss-Legendre rule and writes it to a file.. 
%
%    The user specifies:
%    * the spatial dimension of the quadrature region,
%    * the level that defines the Smolyak grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    08 October 2007
%
%  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.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE_GRID_GL_DATASET\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Compute the abscissas and weights of a quadrature rule\n' );
  fprintf ( 1, '  associated with a sparse grid derived from a Smolyak\n' );
  fprintf ( 1, '  construction based on 1D Gauss-Legendre rules.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Inputs to the program include:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    DIM_NUM, the spatial dimension.\n' );
  fprintf ( 1, '    (typically in the range of 2 to 10)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    LEVEL_MAX, the "level" of the sparse grid.\n' );
  fprintf ( 1, '    (typically in the range of 0, 1, 2, 3, ...\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Output from the program includes:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    A printed table of the abscissas and weights.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    * A set of 3 files that define the quadrature rule.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    "gl_d?_level?_r.txt", a file of the ranges.\n' );
  fprintf ( 1, '    "gl_d?_level?_w.txt", a file of the weights;\n' );
  fprintf ( 1, '    "gl_d?_level?_x.txt", a file of the abscissas;\n' );
%
%  Get the spatial dimension.
%
  if ( nargin < 1)
    fprintf ( 1, '\n' );
    dim_num = input ( '  Enter the value of DIM_NUM (1 or greater): ' );
  elseif ( ischar ( dim_num ) )
    dim_num = str2num ( dim_num );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension requested is = %d\n', dim_num );
%
%  Get the level.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    level_max = input ( '  Enter the value of LEVEL_MAX (0 or greater): ' );
  elseif ( ischar ( level_max ) )
    level_max = str2num ( level_max );
  end

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

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEVEL_MIN = %d\n', level_min );
  fprintf ( 1, '  LEVEL_MAX = %d\n', level_max );
%
%  How many distinct points will there be?
%
  point_num = sparse_grid_gl_size ( dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of distinct abscissas in the\n' );
  fprintf ( 1, '  quadrature rule is determined from the spatial\n' );
  fprintf ( 1, '  dimension DIM_NUM and the level LEVEL_MAX.\n' );
  fprintf ( 1, '  For the given input, this value will be = %d\n', point_num );

  r = zeros ( dim_num, 2 );
%
%  Compute the weights and points.
%
  r(1:dim_num,1) = -1.0;
  r(1:dim_num,2) = +1.0;

  [ w, x ] = sparse_grid_gl ( dim_num, level_max, point_num );

  r8mat_transpose_print_some ( dim_num, point_num, x, 1, 1, ...
    dim_num, 10, '  First 10 grid points:' );

  r8vec_print_some ( point_num, w, 1, 10, '  First 10 weights:' );

  weight_sum = sum ( w(1:point_num) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Weights sum to   %24.16f\n', weight_sum );
  fprintf ( 1, '  Correct value is %24.16f\n', 2.0^dim_num );
%
%  Write the rule to files.
%
  r_filename = sprintf ( 'gl_d%d_level%d_r.txt', dim_num, level_max );
  w_filename = sprintf ( 'gl_d%d_level%d_w.txt', dim_num, level_max );
  x_filename = sprintf ( 'gl_d%d_level%d_x.txt', dim_num, level_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating R file = "%s".\n', r_filename );

  r8mat_write ( r_filename, dim_num, 2, r );

  fprintf ( 1, '  Creating W file = "%s".\n', w_filename );

  r8mat_write ( w_filename, 1, point_num, w );

  fprintf ( 1, '  Creating X file = "%s".\n', x_filename );

  r8mat_write ( x_filename, dim_num, point_num, x );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPARSE)GRID_GL_DATASET:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function [ a, more, h, t ] = comp_next ( n, k, a, more, h, t )

%*****************************************************************************80
%
%% COMP_NEXT computes the compositions of the integer N into K parts.
%
%  Discussion:
%
%    A composition of the integer N into K parts is an ordered sequence
%    of K nonnegative integers which sum to N.  The compositions (1,2,1)
%    and (1,1,2) are considered to be distinct.
%
%    The routine computes one composition on each call until there are no more.
%    For instance, one composition of 6 into 3 parts is
%    3+2+1, another would be 6+0+0.
%
%    On the first call to this routine, set MORE = FALSE.  The routine
%    will compute the first element in the sequence of compositions, and
%    return it, as well as setting MORE = TRUE.  If more compositions
%    are desired, call again, and again.  Each time, the routine will
%    return with a new composition.
%
%    However, when the LAST composition in the sequence is computed 
%    and returned, the routine will reset MORE to FALSE, signaling that
%    the end of the sequence has been reached.
%
%    This routine originally used a SAVE statement to maintain the
%    variables H and T.  I have decided that it is safer
%    to pass these variables as arguments, even though the user should
%    never alter them.  This allows this routine to safely shuffle
%    between several ongoing calculations.
%
%    There are 28 compositions of 6 into three parts.  This routine will
%    produce those compositions in the following order:
%
%     I         A
%     -     ---------
%     1     6   0   0
%     2     5   1   0
%     3     4   2   0
%     4     3   3   0
%     5     2   4   0
%     6     1   5   0
%     7     0   6   0
%     8     5   0   1
%     9     4   1   1
%    10     3   2   1
%    11     2   3   1
%    12     1   4   1
%    13     0   5   1
%    14     4   0   2
%    15     3   1   2
%    16     2   2   2
%    17     1   3   2
%    18     0   4   2
%    19     3   0   3
%    20     2   1   3
%    21     1   2   3
%    22     0   3   3
%    23     2   0   4
%    24     1   1   4
%    25     0   2   4
%    26     1   0   5
%    27     0   1   5
%    28     0   0   6
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    02 July 2008
%
%  Author:
%
%    Original FORTRAN77 version by Albert Nijenhuis and Herbert Wilf
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Albert Nijenhuis, Herbert Wilf,
%    Combinatorial Algorithms for Computers and Calculators,
%    Second Edition,
%    Academic Press, 1978,
%    ISBN: 0-12-519260-6,
%    LC: QA164.N54.
%
%  Parameters:
%
%    Input, integer N, the integer whose compositions are desired.
%
%    Input, integer K, the number of parts in the composition.
%
%    Input, integer A(K), the previous composition.  On the first call,
%    with MORE = FALSE, set A = [].  Thereafter, A should be the 
%    value of A output from the previous call.
%
%    Input, logical MORE.  The input value of MORE on the first
%    call should be FALSE, which tells the program to initialize.
%    On subsequent calls, MORE should be TRUE, or simply the
%    output value of MORE from the previous call.
%
%    Input, integer H, T, two internal parameters needed for the
%    computation.  The user may need to initialize these before the
%    very first call, but these initial values are not important.
%    The user should not alter these parameters once the computation
%    begins.
%
%    Output, integer A(K), the next composition.
%
%    Output, logical MORE, will be TRUE unless the composition 
%    that is being returned is the final one in the sequence.
%
%    Output, integer H, T, the updated values of the two internal 
%    parameters.
%
  if ( ~more )

    t = n;
    h = 0;
    a(1) = n;
    a(2:k) = 0;

  else
      
    if ( 1 < t )
      h = 0;
    end

    h = h + 1;
    t = a(h);
    a(h) = 0;
    a(1) = t - 1;
    a(h+1) = a(h+1) + 1;

  end

  more = ( a(k) ~= n );

  return
end
function grid_point = gl_abscissa ( dim_num, point_num, grid_index, grid_base )

%*****************************************************************************80
%
%% GL_ABSCISSA sets abscissas for multidimensional Gauss-Legendre quadrature.
%
%  Discussion:
%
%    The "nesting" as it occurs for Gauss-Legendre sparse grids simply
%    involves the use of a specified set of permissible orders for the
%    rule.  
%
%    The X array lists the (complete) Gauss-Legendre abscissas for rules 
%    of order 1, 3, 7, 15, 31, 63 and 127, in order. 
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), for each
%    point and dimension, the index of the Gauss-Legendre abscissa.
%
%    Input, integer GRID_BASE(DIM_NUM), the order of the 
%    Gauss-Legendre rule being used in each dimension.
%
%    Output, real GRID_POINT(DIM_NUM), the grid points of
%    Gauss-Legendre abscissas.
%
  skip = [ 0, 1, 4, 11, 26, 57, 120, 247 ];
  x = [ ...
       0.0E+00, ...
     - 0.774596669241483377035853079956E+00, ...
       0.0E+00, ...
       0.774596669241483377035853079956E+00, ...
     - 0.949107912342758524526189684048E+00, ...
     - 0.741531185599394439863864773281E+00, ...
     - 0.405845151377397166906606412077E+00, ...
       0.0E+00, ...
       0.405845151377397166906606412077E+00, ...
       0.741531185599394439863864773281E+00, ...
       0.949107912342758524526189684048E+00, ...
      - 0.987992518020485428489565718587E+00, ...
      - 0.937273392400705904307758947710E+00, ...
      - 0.848206583410427216200648320774E+00, ...
      - 0.724417731360170047416186054614E+00, ...
      - 0.570972172608538847537226737254E+00, ...
      - 0.394151347077563369897207370981E+00, ...
      - 0.201194093997434522300628303395E+00, ...
        0.0E+00, ...
        0.201194093997434522300628303395E+00, ...
       0.394151347077563369897207370981E+00, ...
       0.570972172608538847537226737254E+00, ...
       0.724417731360170047416186054614E+00, ...
       0.848206583410427216200648320774E+00, ...
       0.937273392400705904307758947710E+00, ...
       0.987992518020485428489565718587E+00, ...
      -0.99708748181947707454263838179654, ...   
      -0.98468590966515248400211329970113, ...   
      -0.96250392509294966178905249675943, ...   
      -0.93075699789664816495694576311725, ... 
      -0.88976002994827104337419200908023, ... 
      -0.83992032014626734008690453594388, ... 
      -0.78173314841662494040636002019484, ... 
      -0.71577678458685328390597086536649, ...
      -0.64270672292426034618441820323250, ...  
      -0.56324916140714926272094492359516, ...   
      -0.47819378204490248044059403935649, ...  
      -0.38838590160823294306135146128752, ...    
      -0.29471806998170161661790389767170, ... 
      -0.19812119933557062877241299603283, ... 
      -0.99555312152341520325174790118941E-01, ...
       0.00000000000000000000000000000000, ... 
       0.99555312152341520325174790118941E-01, ...
       0.19812119933557062877241299603283, ...   
       0.29471806998170161661790389767170, ...   
       0.38838590160823294306135146128752, ...  
       0.47819378204490248044059403935649, ...   
       0.56324916140714926272094492359516, ...
       0.64270672292426034618441820323250, ... 
       0.71577678458685328390597086536649, ...
       0.78173314841662494040636002019484, ...   
       0.83992032014626734008690453594388, ...    
       0.88976002994827104337419200908023, ...    
       0.93075699789664816495694576311725, ...    
       0.96250392509294966178905249675943, ...    
       0.98468590966515248400211329970113, ...
       0.99708748181947707454263838179654, ...
      -0.99928298402912378050701628988630E+00, ...    
      -0.99622401277797010860209018267357E+00, ...    
      -0.99072854689218946681089469460884E+00, ...    
      -0.98280881059372723486251140727639E+00, ...    
      -0.97248403469757002280196067864927E+00, ...    
      -0.95977944975894192707035416626398E+00, ...    
      -0.94472613404100980296637531962798E+00, ...    
      -0.92736092062184320544703138132518E+00, ...    
      -0.90772630277853155803695313291596E+00, ...    
      -0.88587032850785342629029845731337E+00, ...    
      -0.86184648236412371953961183943106E+00, ...    
      -0.83571355431950284347180776961571E+00, ...    
      -0.80753549577345676005146598636324E+00, ...    
      -0.77738126299037233556333018991104E+00, ...    
      -0.74532464831784741782932166103759E+00, ...    
      -0.71144409958484580785143153770401E+00, ...    
      -0.67582252811498609013110331596954E+00, ...    
      -0.63854710582136538500030695387338E+00, ...    
      -0.59970905187762523573900892686880E+00, ...    
      -0.55940340948628501326769780007005E+00, ...    
      -0.51772881329003324812447758452632E+00, ...    
      -0.47478724799480439992221230985149E+00, ...    
      -0.43068379879511160066208893391863E+00, ...    
      -0.38552639421224789247761502227440E+00, ...    
      -0.33942554197458440246883443159432E+00, ...    
      -0.29249405858625144003615715555067E+00, ...    
      -0.24484679324595336274840459392483E+00, ...    
      -0.19660034679150668455762745706572E+00, ...    
      -0.14787278635787196856983909655297E+00, ...    
      -0.98783356446945279529703669453922E-01, ...
      -0.49452187116159627234233818051808E-01, ...
       0.00000000000000000000000000000000E+00, ...    
       0.49452187116159627234233818051808E-01, ...
       0.98783356446945279529703669453922E-01, ...
       0.14787278635787196856983909655297E+00, ...    
       0.19660034679150668455762745706572E+00, ...    
       0.24484679324595336274840459392483E+00, ...    
       0.29249405858625144003615715555067E+00, ...    
       0.33942554197458440246883443159432E+00, ...    
       0.38552639421224789247761502227440E+00, ...    
       0.43068379879511160066208893391863E+00, ...    
       0.47478724799480439992221230985149E+00, ...    
       0.51772881329003324812447758452632E+00, ...    
       0.55940340948628501326769780007005E+00, ...    
       0.59970905187762523573900892686880E+00, ...    
       0.63854710582136538500030695387338E+00, ...    
       0.67582252811498609013110331596954E+00, ...    
       0.71144409958484580785143153770401E+00, ...    
       0.74532464831784741782932166103759E+00, ...    
       0.77738126299037233556333018991104E+00, ...    
       0.80753549577345676005146598636324E+00, ...    
       0.83571355431950284347180776961571E+00, ...    
       0.86184648236412371953961183943106E+00, ...    
       0.88587032850785342629029845731337E+00, ...    
       0.90772630277853155803695313291596E+00, ...    
       0.92736092062184320544703138132518E+00, ...    
       0.94472613404100980296637531962798E+00, ...    
       0.95977944975894192707035416626398E+00, ...    
       0.97248403469757002280196067864927E+00, ...    
       0.98280881059372723486251140727639E+00, ...    
       0.99072854689218946681089469460884E+00, ...    
       0.99622401277797010860209018267357E+00, ... 
       0.99928298402912378050701628988630E+00, ...
      -0.99982213041530614629963254927125E+00, ...    
      -0.99906293435531189513828920479421E+00, ...    
      -0.99769756618980462107441703193392E+00, ...    
      -0.99572655135202722663543337085008E+00, ...    
      -0.99315104925451714736113079489080E+00, ...    
      -0.98997261459148415760778669967548E+00, ...    
      -0.98619317401693166671043833175407E+00, ...    
      -0.98181502080381411003346312451200E+00, ...    
      -0.97684081234307032681744391886221E+00, ...    
      -0.97127356816152919228894689830512E+00, ...    
      -0.96511666794529212109082507703391E+00, ...    
      -0.95837384942523877114910286998060E+00, ...    
      -0.95104920607788031054790764659636E+00, ...    
      -0.94314718462481482734544963026201E+00, ...    
      -0.93467258232473796857363487794906E+00, ...    
      -0.92563054405623384912746466814259E+00, ...    
      -0.91602655919146580931308861741716E+00, ...    
      -0.90586645826182138280246131760282E+00, ...    
      -0.89515640941708370896904382642451E+00, ...    
      -0.88390291468002656994525794802849E+00, ...    
      -0.87211280599856071141963753428864E+00, ...    
      -0.85979324109774080981203134414483E+00, ...    
      -0.84695169913409759845333931085437E+00, ...    
      -0.83359597615489951437955716480123E+00, ...    
      -0.81973418036507867415511910167470E+00, ...    
      -0.80537472720468021466656079404644E+00, ...    
      -0.79052633423981379994544995252740E+00, ...    
      -0.77519801587020238244496276354566E+00, ...    
      -0.75939907785653667155666366659810E+00, ...    
      -0.74313911167095451292056688997595E+00, ...   
      -0.72642798867407268553569290153270E+00, ...    
      -0.70927585412210456099944463906757E+00, ...    
      -0.69169312100770067015644143286666E+00, ...    
      -0.67369046373825048534668253831602E+00, ...    
      -0.65527881165548263027676505156852E+00, ...    
      -0.63646934240029724134760815684175E+00, ...    
      -0.61727347512685828385763916340822E+00, ...    
      -0.59770286357006522938441201887478E+00, ...    
      -0.57776938897061258000325165713764E+00, ...    
      -0.55748515286193223292186190687872E+00, ...    
      -0.53686246972339756745816636353452E+00, ...    
      -0.51591385950424935727727729906662E+00, ...    
      -0.49465204002278211739494017368636E+00, ...    
      -0.47308991924540524164509989939699E+00, ...    
      -0.45124058745026622733189858020729E+00, ...    
      -0.42911730928019337626254405355418E+00, ...    
      -0.40673351568978256340867288124339E+00, ...    
      -0.38410279579151693577907781452239E+00, ...    
      -0.36123888860586970607092484346723E+00, ...    
      -0.33815567472039850137600027657095E+00, ...    
      -0.31486716786289498148601475374890E+00, ...    
      -0.29138750639370562079451875284568E+00, ...    
      -0.26773094472238862088834352027938E+00, ...    
      -0.24391184465391785797071324453138E+00, ...    
      -0.21994466666968754245452337866940E+00, ...    
      -0.19584396114861085150428162519610E+00, ...    
      -0.17162435953364216500834492248954E+00, ...    
      -0.14730056544908566938932929319807E+00, ...    
      -0.12288734577408297172603365288567E+00, ...    
      -0.98399521677698970751091751509101E-01, ...
      -0.73851959621048545273440409360569E-01, ...
      -0.49259562331926630315379321821927E-01, ...
      -0.24637259757420944614897071846088E-01, ...
       0.00000000000000000000000000000000E+00, ...    
       0.24637259757420944614897071846088E-01, ...
       0.49259562331926630315379321821927E-01, ...
       0.73851959621048545273440409360569E-01, ...
       0.98399521677698970751091751509101E-01, ...
       0.12288734577408297172603365288567E+00, ...    
       0.14730056544908566938932929319807E+00, ...    
       0.17162435953364216500834492248954E+00, ...    
       0.19584396114861085150428162519610E+00, ...    
       0.21994466666968754245452337866940E+00, ...    
       0.24391184465391785797071324453138E+00, ...    
       0.26773094472238862088834352027938E+00, ...    
       0.29138750639370562079451875284568E+00, ...    
       0.31486716786289498148601475374890E+00, ...    
       0.33815567472039850137600027657095E+00, ...    
       0.36123888860586970607092484346723E+00, ...    
       0.38410279579151693577907781452239E+00, ...    
       0.40673351568978256340867288124339E+00, ...    
       0.42911730928019337626254405355418E+00, ...    
       0.45124058745026622733189858020729E+00, ...    
       0.47308991924540524164509989939699E+00, ...    
       0.49465204002278211739494017368636E+00, ...    
       0.51591385950424935727727729906662E+00, ...    
       0.53686246972339756745816636353452E+00, ...    
       0.55748515286193223292186190687872E+00, ...    
       0.57776938897061258000325165713764E+00, ...    
       0.59770286357006522938441201887478E+00, ...    
       0.61727347512685828385763916340822E+00, ...    
       0.63646934240029724134760815684175E+00, ...    
       0.65527881165548263027676505156852E+00, ...    
       0.67369046373825048534668253831602E+00, ...    
       0.69169312100770067015644143286666E+00, ...   
       0.70927585412210456099944463906757E+00, ...    
       0.72642798867407268553569290153270E+00, ...    
       0.74313911167095451292056688997595E+00, ...    
       0.75939907785653667155666366659810E+00, ...    
       0.77519801587020238244496276354566E+00, ...    
       0.79052633423981379994544995252740E+00, ...    
       0.80537472720468021466656079404644E+00, ...    
       0.81973418036507867415511910167470E+00, ...    
       0.83359597615489951437955716480123E+00, ...    
       0.84695169913409759845333931085437E+00, ...    
       0.85979324109774080981203134414483E+00, ...    
       0.87211280599856071141963753428864E+00, ...    
       0.88390291468002656994525794802849E+00, ...    
       0.89515640941708370896904382642451E+00, ...    
       0.90586645826182138280246131760282E+00, ...    
       0.91602655919146580931308861741716E+00, ...    
       0.92563054405623384912746466814259E+00, ...    
       0.93467258232473796857363487794906E+00, ...    
       0.94314718462481482734544963026201E+00, ...    
       0.95104920607788031054790764659636E+00, ...    
       0.95837384942523877114910286998060E+00, ...    
       0.96511666794529212109082507703391E+00, ...    
       0.97127356816152919228894689830512E+00, ...    
       0.97684081234307032681744391886221E+00, ...    
       0.98181502080381411003346312451200E+00, ...    
       0.98619317401693166671043833175407E+00, ...    
       0.98997261459148415760778669967548E+00, ...    
       0.99315104925451714736113079489080E+00, ...    
       0.99572655135202722663543337085008E+00, ...    
       0.99769756618980462107441703193392E+00, ...    
       0.99906293435531189513828920479421E+00, ...    
       0.99982213041530614629963254927125E+00  ...    
  ];

  if ( any ( grid_base(1:dim_num) < 0 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GL_ABSCISSA - Fatal error!\n' );
    fprintf ( 1, '  Some base values are less than 0.\n' );
    error ( 'GL_ABSCISSA - Fatal error!' );
  end

  if ( any ( 63 < grid_base(1:dim_num) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'GL_ABSCISSA - Fatal error!\n' );
    fprintf ( 1, '  Some base values are greater than 63.\n' );
    error ( 'GL_ABSCISSA - Fatal error!' );
  end

  for point = 1 : point_num
    for dim = 1 : dim_num

      level = i4_log_2 ( grid_base(dim) + 1 );

      pointer = skip(level+1) + ( grid_index(dim,point) + grid_base(dim) + 1 );

      grid_point(dim,point) = x(pointer);

    end
  end

  return
end
function weight = gl_weights ( order )

%*****************************************************************************80
%
%% GL_WEIGHTS returns weights for certain Gauss-Legendre quadrature rules.
%
%  Discussion:
%
%    The allowed orders are 1, 3, 7, 15, 31, 63 or 127.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    02 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Arthur Stroud, Don Secrest,
%    Gaussian Quadrature Formulas,
%    Prentice Hall, 1966,
%    LC: QA299.4G3S7.
%
%  Parameters:
%
%    Input, integer ORDER, the order of the rule.
%    ORDER must be 1, 3, 7, 15, 31, 63, or 127.
%
%    Output, real WEIGHT(ORDER), the weights.
%    The weights are positive, symmetric and should sum to 2.
%
  if ( order == 1 )

    weight(1) = 2.0E+00;

  elseif ( order == 3 )

    weight(1) = 0.555555555555555555555555555555E+00;
    weight(2) = 0.888888888888888888888888888888E+00;
    weight(3) = 0.555555555555555555555555555555E+00;

  elseif ( order == 7 )

    weight(1) = 0.129484966168869693270611432679E+00;
    weight(2) = 0.279705391489276667901467771424E+00;
    weight(3) = 0.381830050505118944950369775489E+00;
    weight(4) = 0.417959183673469387755102040816E+00;
    weight(5) = 0.381830050505118944950369775489E+00;
    weight(6) = 0.279705391489276667901467771424E+00;
    weight(7) = 0.129484966168869693270611432679E+00;

  elseif ( order == 15 )

    weight(1) =  0.307532419961172683546283935772E-01;
    weight(2) =  0.703660474881081247092674164507E-01;
    weight(3) =  0.107159220467171935011869546686E+00;
    weight(4) =  0.139570677926154314447804794511E+00;
    weight(5) =  0.166269205816993933553200860481E+00;
    weight(6) =  0.186161000015562211026800561866E+00;
    weight(7) =  0.198431485327111576456118326444E+00;
    weight(8) =  0.202578241925561272880620199968E+00;
    weight(9) =  0.198431485327111576456118326444E+00;
    weight(10) = 0.186161000015562211026800561866E+00;
    weight(11) = 0.166269205816993933553200860481E+00;
    weight(12) = 0.139570677926154314447804794511E+00;
    weight(13) = 0.107159220467171935011869546686E+00;
    weight(14) = 0.703660474881081247092674164507E-01;
    weight(15) = 0.307532419961172683546283935772E-01;

  elseif ( order == 31 )
 
    weight( 1) =   0.74708315792487746093913218970494E-02;
    weight( 2) =   0.17318620790310582463552990782414E-01;
    weight( 3) =   0.27009019184979421800608642617676E-01;
    weight( 4) =   0.36432273912385464024392008749009E-01;
    weight( 5) =   0.45493707527201102902315857856518E-01;
    weight( 6) =   0.54103082424916853711666259085477E-01;
    weight( 7) =   0.62174786561028426910343543686657E-01;
    weight( 8) =   0.69628583235410366167756126255124E-01;
    weight( 9) =   0.76390386598776616426357674901331E-01;
    weight(10) =   0.82392991761589263903823367431962E-01;
    weight(11) =   0.87576740608477876126198069695333E-01;
    weight(12) =   0.91890113893641478215362871607150E-01;
    weight(13) =   0.95290242912319512807204197487597E-01;
    weight(14) =   0.97743335386328725093474010978997E-01;
    weight(15) =   0.99225011226672307874875514428615E-01;
    weight(16) =   0.99720544793426451427533833734349E-01;
    weight(17) =   0.99225011226672307874875514428615E-01;
    weight(18) =   0.97743335386328725093474010978997E-01;
    weight(19) =   0.95290242912319512807204197487597E-01;
    weight(20) =   0.91890113893641478215362871607150E-01;
    weight(21) =   0.87576740608477876126198069695333E-01;
    weight(22) =   0.82392991761589263903823367431962E-01;
    weight(23) =   0.76390386598776616426357674901331E-01;
    weight(24) =   0.69628583235410366167756126255124E-01;
    weight(25) =   0.62174786561028426910343543686657E-01;
    weight(26) =   0.54103082424916853711666259085477E-01;
    weight(27) =   0.45493707527201102902315857856518E-01;
    weight(28) =   0.36432273912385464024392008749009E-01;
    weight(29) =   0.27009019184979421800608642617676E-01;
    weight(30) =   0.17318620790310582463552990782414E-01;
    weight(31) =   0.74708315792487746093913218970494E-02;

  elseif ( order == 63 )

    weight( 1) =   0.18398745955770837880499331680577E-02;
    weight( 2) =   0.42785083468637618661951422543371E-02;
    weight( 3) =   0.67102917659601362519069109850892E-02;
    weight( 4) =   0.91259686763266563540586445877022E-02;
    weight( 5) =   0.11519376076880041750750606118707E-01;
    weight( 6) =   0.13884612616115610824866086365937E-01;
    weight( 7) =   0.16215878410338338882283672974995E-01;
    weight( 8) =   0.18507464460161270409260545805144E-01;
    weight( 9) =   0.20753761258039090775341953421471E-01;
    weight(10) =   0.22949271004889933148942319561770E-01;
    weight(11) =   0.25088620553344986618630138068443E-01;
    weight(12) =   0.27166574359097933225189839439413E-01;
    weight(13) =   0.29178047208280526945551502154029E-01;
    weight(14) =   0.31118116622219817508215988557189E-01;
    weight(15) =   0.32982034883779341765683179672459E-01;
    weight(16) =   0.34765240645355877697180504642788E-01;
    weight(17) =   0.36463370085457289630452409787542E-01;
    weight(18) =   0.38072267584349556763638324927889E-01;
    weight(19) =   0.39587995891544093984807928149202E-01;
    weight(20) =   0.41006845759666398635110037009072E-01;
    weight(21) =   0.42325345020815822982505485403028E-01;
    weight(22) =   0.43540267083027590798964315704401E-01;
    weight(23) =   0.44648638825941395370332669516813E-01;
    weight(24) =   0.45647747876292608685885992608542E-01;
    weight(25) =   0.46535149245383696510395418746953E-01;
    weight(26) =   0.47308671312268919080604988338844E-01;
    weight(27) =   0.47966421137995131411052756195132E-01;
    weight(28) =   0.48506789097883847864090099145802E-01;
    weight(29) =   0.48928452820511989944709361549215E-01;
    weight(30) =   0.49230380423747560785043116988145E-01;
    weight(31) =   0.49411833039918178967039646116705E-01;
    weight(32) =   0.49472366623931020888669360420926E-01;
    weight(33) =   0.49411833039918178967039646116705E-01;
    weight(34) =   0.49230380423747560785043116988145E-01;
    weight(35) =   0.48928452820511989944709361549215E-01;
    weight(36) =   0.48506789097883847864090099145802E-01;
    weight(37) =   0.47966421137995131411052756195132E-01;
    weight(38) =   0.47308671312268919080604988338844E-01;
    weight(39) =   0.46535149245383696510395418746953E-01;
    weight(40) =   0.45647747876292608685885992608542E-01;
    weight(41) =   0.44648638825941395370332669516813E-01;
    weight(42) =   0.43540267083027590798964315704401E-01;
    weight(43) =   0.42325345020815822982505485403028E-01;
    weight(44) =   0.41006845759666398635110037009072E-01;
    weight(45) =   0.39587995891544093984807928149202E-01;
    weight(46) =   0.38072267584349556763638324927889E-01;
    weight(47) =   0.36463370085457289630452409787542E-01;
    weight(48) =   0.34765240645355877697180504642788E-01;
    weight(49) =   0.32982034883779341765683179672459E-01;
    weight(50) =   0.31118116622219817508215988557189E-01;
    weight(51) =   0.29178047208280526945551502154029E-01;
    weight(52) =   0.27166574359097933225189839439413E-01;
    weight(53) =   0.25088620553344986618630138068443E-01;
    weight(54) =   0.22949271004889933148942319561770E-01;
    weight(55) =   0.20753761258039090775341953421471E-01;
    weight(56) =   0.18507464460161270409260545805144E-01;
    weight(57) =   0.16215878410338338882283672974995E-01;
    weight(58) =   0.13884612616115610824866086365937E-01;
    weight(59) =   0.11519376076880041750750606118707E-01;
    weight(60) =   0.91259686763266563540586445877022E-02;
    weight(61) =   0.67102917659601362519069109850892E-02;
    weight(62) =   0.42785083468637618661951422543371E-02;
    weight(63) =   0.18398745955770837880499331680577E-02;
 
  elseif ( order == 127 )
  
    weight(  1) =   0.45645726109586654495731936146574E-03;
    weight(  2) =   0.10622766869538486959954760554099E-02;
    weight(  3) =   0.16683488125171936761028811985672E-02;
    weight(  4) =   0.22734860707492547802810838362671E-02;
    weight(  5) =   0.28772587656289004082883197417581E-02;
    weight(  6) =   0.34792893810051465908910894094105E-02;
    weight(  7) =   0.40792095178254605327114733456293E-02;
    weight(  8) =   0.46766539777779034772638165662478E-02;
    weight(  9) =   0.52712596565634400891303815906251E-02;
    weight( 10) =   0.58626653903523901033648343751367E-02;
    weight( 11) =   0.64505120486899171845442463868748E-02;
    weight( 12) =   0.70344427036681608755685893032552E-02;
    weight( 13) =   0.76141028256526859356393930849227E-02;
    weight( 14) =   0.81891404887415730817235884718726E-02;
    weight( 15) =   0.87592065795403145773316804234385E-02;
    weight( 16) =   0.93239550065309714787536985834029E-02;
    weight( 17) =   0.98830429087554914716648010899606E-02;
    weight( 18) =   0.10436130863141005225673171997668E-01;
    weight( 19) =   0.10982883090068975788799657376065E-01;
    weight( 20) =   0.11522967656921087154811609734510E-01;
    weight( 21) =   0.12056056679400848183529562144697E-01;
    weight( 22) =   0.12581826520465013101514365424172E-01;
    weight( 23) =   0.13099957986718627426172681912499E-01;
    weight( 24) =   0.13610136522139249906034237533759E-01;
    weight( 25) =   0.14112052399003395774044161633613E-01;
    weight( 26) =   0.14605400905893418351737288078952E-01;
    weight( 27) =   0.15089882532666922992635733981431E-01;
    weight( 28) =   0.15565203152273955098532590262975E-01;
    weight( 29) =   0.16031074199309941802254151842763E-01;
    weight( 30) =   0.16487212845194879399346060358146E-01;
    weight( 31) =   0.16933342169871654545878815295200E-01;
    weight( 32) =   0.17369191329918731922164721250350E-01;
    weight( 33) =   0.17794495722974774231027912900351E-01;
    weight( 34) =   0.18208997148375106468721469154479E-01;
    weight( 35) =   0.18612443963902310429440419898958E-01;
    weight( 36) =   0.19004591238555646611148901044533E-01;
    weight( 37) =   0.19385200901246454628112623489471E-01;
    weight( 38) =   0.19754041885329183081815217323169E-01;
    weight( 39) =   0.20110890268880247225644623956287E-01;
    weight( 40) =   0.20455529410639508279497065713301E-01;
    weight( 41) =   0.20787750081531811812652137291250E-01;
    weight( 42) =   0.21107350591688713643523847921658E-01;
    weight( 43) =   0.21414136912893259295449693233545E-01;
    weight( 44) =   0.21707922796373466052301324695331E-01;
    weight( 45) =   0.21988529885872983756478409758807E-01;
    weight( 46) =   0.22255787825930280235631416460158E-01;
    weight( 47) =   0.22509534365300608085694429903050E-01;
    weight( 48) =   0.22749615455457959852242553240982E-01;
    weight( 49) =   0.22975885344117206754377437838947E-01;
    weight( 50) =   0.23188206663719640249922582981729E-01;
    weight( 51) =   0.23386450514828194170722043496950E-01;
    weight( 52) =   0.23570496544381716050033676844306E-01;
    weight( 53) =   0.23740233018760777777714726703424E-01;
    weight( 54) =   0.23895556891620665983864481754172E-01;
    weight( 55) =   0.24036373866450369675132086026456E-01;
    weight( 56) =   0.24162598453819584716522917710986E-01;
    weight( 57) =   0.24274154023278979833195063936748E-01;
    weight( 58) =   0.24370972849882214952813561907241E-01;
    weight( 59) =   0.24452996155301467956140198471529E-01;
    weight( 60) =   0.24520174143511508275183033290175E-01;
    weight( 61) =   0.24572466031020653286354137335186E-01;
    weight( 62) =   0.24609840071630254092545634003360E-01;
    weight( 63) =   0.24632273575707679066033370218017E-01;
    weight( 64) =   0.24639752923961094419579417477503E-01;
    weight( 65) =   0.24632273575707679066033370218017E-01;
    weight( 66) =   0.24609840071630254092545634003360E-01;
    weight( 67) =   0.24572466031020653286354137335186E-01;
    weight( 68) =   0.24520174143511508275183033290175E-01;
    weight( 69) =   0.24452996155301467956140198471529E-01;
    weight( 70) =   0.24370972849882214952813561907241E-01;
    weight( 71) =   0.24274154023278979833195063936748E-01;
    weight( 72) =   0.24162598453819584716522917710986E-01;
    weight( 73) =   0.24036373866450369675132086026456E-01;
    weight( 74) =   0.23895556891620665983864481754172E-01;
    weight( 75) =   0.23740233018760777777714726703424E-01;
    weight( 76) =   0.23570496544381716050033676844306E-01;
    weight( 77) =   0.23386450514828194170722043496950E-01;
    weight( 78) =   0.23188206663719640249922582981729E-01;
    weight( 79) =   0.22975885344117206754377437838947E-01;
    weight( 80) =   0.22749615455457959852242553240982E-01;
    weight( 81) =   0.22509534365300608085694429903050E-01;
    weight( 82) =   0.22255787825930280235631416460158E-01;
    weight( 83) =   0.21988529885872983756478409758807E-01;
    weight( 84) =   0.21707922796373466052301324695331E-01;
    weight( 85) =   0.21414136912893259295449693233545E-01;
    weight( 86) =   0.21107350591688713643523847921658E-01;
    weight( 87) =   0.20787750081531811812652137291250E-01;
    weight( 88) =   0.20455529410639508279497065713301E-01;
    weight( 89) =   0.20110890268880247225644623956287E-01;
    weight( 90) =   0.19754041885329183081815217323169E-01;
    weight( 91) =   0.19385200901246454628112623489471E-01;
    weight( 92) =   0.19004591238555646611148901044533E-01;
    weight( 93) =   0.18612443963902310429440419898958E-01;
    weight( 94) =   0.18208997148375106468721469154479E-01;
    weight( 95) =   0.17794495722974774231027912900351E-01;
    weight( 96) =   0.17369191329918731922164721250350E-01;
    weight( 97) =   0.16933342169871654545878815295200E-01;
    weight( 98) =   0.16487212845194879399346060358146E-01;
    weight( 99) =   0.16031074199309941802254151842763E-01;
    weight(100) =   0.15565203152273955098532590262975E-01;
    weight(101) =   0.15089882532666922992635733981431E-01;
    weight(102) =   0.14605400905893418351737288078952E-01;
    weight(103) =   0.14112052399003395774044161633613E-01;
    weight(104) =   0.13610136522139249906034237533759E-01;
    weight(105) =   0.13099957986718627426172681912499E-01;
    weight(106) =   0.12581826520465013101514365424172E-01;
    weight(107) =   0.12056056679400848183529562144697E-01;
    weight(108) =   0.11522967656921087154811609734510E-01;
    weight(109) =   0.10982883090068975788799657376065E-01;
    weight(110) =   0.10436130863141005225673171997668E-01;
    weight(111) =   0.98830429087554914716648010899606E-02;
    weight(112) =   0.93239550065309714787536985834029E-02;
    weight(113) =   0.87592065795403145773316804234385E-02;
    weight(114) =   0.81891404887415730817235884718726E-02;
    weight(115) =   0.76141028256526859356393930849227E-02;
    weight(116) =   0.70344427036681608755685893032552E-02;
    weight(117) =   0.64505120486899171845442463868748E-02;
    weight(118) =   0.58626653903523901033648343751367E-02;
    weight(119) =   0.52712596565634400891303815906251E-02;
    weight(120) =   0.46766539777779034772638165662478E-02;
    weight(121) =   0.40792095178254605327114733456293E-02;
    weight(122) =   0.34792893810051465908910894094105E-02;
    weight(123) =   0.28772587656289004082883197417581E-02;
    weight(124) =   0.22734860707492547802810838362671E-02;
    weight(125) =   0.16683488125171936761028811985672E-02;
    weight(126) =   0.10622766869538486959954760554099E-02;
    weight(127) =   0.45645726109586654495731936146574E-03;
 
  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'GL_WEIGHTS - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of ORDER = %d\n', order );
    fprintf ( 1, '  Legal values are 1, 3, 7, 15, 31, 63 or 127.\n' );
    error ( 'GL_WEIGHTS - Fatal error!' );

  end

  return
end
function value = i4_choose ( n, k )

%*****************************************************************************80
%
%% I4_CHOOSE computes the binomial coefficient C(N,K).
%
%  Discussion:
%
%    The value is calculated in such a way as to avoid overflow and
%    roundoff.  The calculation is done in integer arithmetic.
%
%    The formula used is:
%
%      C(N,K) = N! / ( K! * (N-K)! )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    ML Wolfson, HV Wright,
%    Algorithm 160:
%    Combinatorial of M Things Taken N at a Time,
%    Communications of the ACM,
%    Volume 6, Number 4, April 1963, page 161.
%
%  Parameters:
%
%    Input, integer N, K, are the values of N and K.
%
%    Output, integer VALUE, the number of combinations of N
%    things taken K at a time.
%
  mn = min ( k, n - k );

  if ( mn < 0 )

    value = 0;

  elseif ( mn == 0 )

    value = 1;

  else

    mx = max ( k, n - k );
    value = mx + 1;

    for i = 2 : mn
      value = ( value * ( mx + i ) ) / i;
    end

  end

  return
end
function value = i4_log_2 ( i )

%*****************************************************************************80
%
%% I4_LOG_2 returns the integer part of the logarithm base 2 of |I|.
%
%  Discussion:
%
%    For positive I4_LOG_2(I), it should be true that
%      2**I4_LOG_2(X) <= |I| < 2**(I4_LOG_2(I)+1).
%    The special case of I4_LOG_2(0) returns -HUGE().
%
%  Example:
%
%     I  Value
%
%     0  -1
%     1,  0
%     2,  1
%     3,  1
%     4,  2
%     5,  2
%     6,  2
%     7,  2
%     8,  3
%     9,  3
%    10,  3
%   127,  6
%   128,  7
%   129,  7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    09 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number whose logarithm base 2 is desired.
%
%    Output, integer VALUE, the integer part of the logarithm base 2 of
%    the absolute value of I.
%
  i = floor ( i );

  if ( i == 0 )

    value = -inf;

  else

    value = 0;

    i = abs ( i );

    while ( 2 <= i )
      i = i / 2;
      value = value + 1;
    end

  end

  return
end
function value = i4_modp ( i, j )

%*****************************************************************************80
%
%% I4_MODP returns the nonnegative remainder of I4 division.
%
%  Discussion:
%
%    If
%      NREM = I4_MODP ( I, J )
%      NMULT = ( I - NREM ) / J
%    then
%      I = J * NMULT + NREM
%    where NREM is always nonnegative.
%
%    The MOD function computes a result with the same sign as the
%    quantity being divided.  Thus, suppose you had an angle A,
%    and you wanted to ensure that it was between 0 and 360.
%    Then mod(A,360) would do, if A was positive, but if A
%    was negative, your result would be between -360 and 0.
%
%    On the other hand, I4_MODP(A,360) is between 0 and 360, always.
%
%  Example:
%
%        I     J     MOD  I4_MODP    Factorization
%
%      107    50       7       7    107 =  2 *  50 + 7
%      107   -50       7       7    107 = -2 * -50 + 7
%     -107    50      -7      43   -107 = -3 *  50 + 43
%     -107   -50      -7      43   -107 =  3 * -50 + 43
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    26 August 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number to be divided.
%
%    Input, integer J, the number that divides I.
%
%    Output, integer VALUE, the nonnegative remainder when I is
%    divided by J.
%
  if ( j == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_MODP - Fatal error!\n' );
    fprintf ( 1, '  Illegal divisor J = %d\n', j );
    error ( 'I4_MODP - Fatal error!' );
  end

  value = mod ( i, j );

  if ( value < 0 )
    value = value + abs ( j );
  end

  return
end
function grid_level = index_level_gl ( level, level_max, dim_num, ...
  point_num, grid_index, grid_base )

%*****************************************************************************80
%
%% INDEX_LEVEL_GL: determine first level at which given index is generated.
%
%  Discussion:
%
%    We are constructing a sparse grid of Gauss-Legendre points.  The grid
%    is built up of product grids, with a characteristic LEVEL.  
%
%    We are concerned with identifying points in this product grid which
%    have actually been generated previously, on a lower value of LEVEL.
%
%    This routine determines the lowest value of LEVEL at which each of
%    the input points would be generated.
%
%    In 1D, given LEVEL, the number of points is ORDER = 2**(LEVEL+1) + 1,
%    (except that LEVEL = 0 implies ORDER = 1%), the BASE is (ORDER-1)/2, 
%    and the point INDEX values range from -BASE to +BASE.
%
%    The values of INDEX and BASE allow us to determine the abstract
%    properties of the point.  In particular, if INDEX is 0, the corresponding
%    Gauss-Legendre abscissa is 0, the special "nested" value we need
%    to take care of.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    26 September 2007
%
%  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 LEVEL, the level at which these points were 
%    generated.  LEVEL_MIN <= LEVEL <= LEVEL_MAX.
%
%    Input, integer LEVEL_MAX, the maximum level.
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer POINT_NUM, the number of points to be tested.
%
%    Input, integer GRID_INDEX(DIM_NUM,POINT_NUM), the indices of the 
%    points to be tested.
%
%    Input, integer GRID_BASE(DIM_NUM), the "base", which is essentially
%    the denominator of the index.
%
%    Output, integer GRID_LEVEL(POINT_NUM), the value of LEVEL at 
%    which the point would first be generated.  This will be the same as
%    the input value of LEVEL, unless the point has an INDEX of 0 and
%    a corresponding BASE that is NOT zero.
%
  level_min = max ( 0, level_max + 1 - dim_num );
%
%  If a point has a DIM-th component whose INDEX is 0, then the 
%  value of LEVEL at which this point would first be generated is
%  less than LEVEL, unless the DIM-th component of GRID_BASE is 0.
%
  for point = 1 : point_num

    grid_level(point) = max ( level, level_min );

    for dim = 1 : dim_num
      if ( grid_index(dim,point) == 0 )
        grid_level(point) = max ( grid_level(point) - grid_base(dim), level_min );
      end
    end

  end

  return
end
function order = level_to_order_open ( dim_num, level )

%*****************************************************************************80
%
%% LEVEL_TO_ORDER converts a level to an order for open rules.
%
%  Discussion:
%
%    Sparse grids can naturally be nested.  A natural scheme is to use
%    a series of one-dimensional rules arranged in a series of "levels"
%    whose order roughly doubles with each step.
%
%    The arrangement described here works naturally for the Fejer Type 1,
%    Fejer Type 2, Newton Cotes Open, Newton Cotes Half Open,
%    and Gauss-Patterson rules.  It also can be used, partially, to describe
%    the growth of Gauss-Legendre rules.
%
%    The idea is that we start with LEVEL = 0, ORDER = 1 indicating the single 
%    point at the center, and for all values afterwards, we use the relationship
%
%      ORDER = 2**(LEVEL+1) - 1.
%
%    The following table shows how the growth will occur:
%
%    Level    Order
%
%    0          1
%    1          3 =  4 - 1
%    2          7 =  8 - 1
%    3         15 = 16 - 1
%    4         31 = 32 - 1
%    5         63 = 64 - 1
%
%    For the Fejer Type 1, Fejer Type 2, Newton Cotes Open, 
%    Newton Cotes Open Half, and Gauss-Patterson rules, the point growth is
%    nested.  If we have ORDER points on a particular LEVEL, the next level 
%    includes all these old points, plus ORDER+1 new points, formed in the 
%    gaps between successive pairs of old points plus an extra point at each 
%    end.
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  4  +  3
%    3         15   =  8  +  7
%    4         31   = 16  + 15
%    5         63   = 32  + 31
%
%    If we use a series of Gauss-Legendre rules, then there is almost no 
%    nesting, except that the central point is shared.  If we insist on 
%    producing a comparable series of such points, then the "nesting" behavior
%    is as follows:
%
%    Level    Order = New + Old
%
%    0          1   =  1  +  0
%    1          3   =  2  +  1
%    2          7   =  6  +  1
%    3         15   = 14  +  1
%    4         31   = 30  +  1
%    5         63   = 62  +  1
%
%    Moreover, if we consider ALL the points used in such a set of "nested" 
%    Gauss-Legendre rules, then we must sum the "NEW" column, and we see that
%    we get roughly twice as many points as for the truly nested rules.
%
%    In this routine, we assume that a vector of levels is given,
%    and the corresponding orders are desired.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    18 April 2007
%
%  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(DIM_NUM), the nesting level.
%
%    Output, integer ORDER(DIM_NUM), the order (number of points) of the rule.
%
  for dim = 1 : dim_num

    if ( level(dim) < 0 )
      order(dim) = -1;
    elseif ( level(dim) == 0 )
      order(dim) = 1;
    else
      order(dim) = 2^( level(dim) + 1 ) - 1;
    end

  end

  return
end
function indx = multigrid_index_z ( dim_num, order_1d, order_nd )

%*****************************************************************************80
%
%% MULTIGRID_INDEX_Z returns an indexed multidimensional grid.
%
%  Discussion:
%
%    For dimension DIM, the number of points is ORDER_1D(DIM).
%
%    We assume that ORDER_1D(DIM) is an odd number,
%      ORDER_1D(DIM) = N = 2 * M + 1
%    so that the points have coordinates
%      -M/M, -(M-1)/M, ..., -1/M, 0/M, 1/M, 2/M, 3/M, ..., (M-1)/M, M/M.
%    and we index them as
%      -M,   -(M-1),        -1,   0,   1,   2,   3,   ...,  M-1,    M.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    26 September 2007
%
%  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 of the points.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the
%    rule in each dimension.
%
%    Input, integer ORDER_ND, the product of the entries of ORDER_1D.
%
%    Output, integer INDX(DIM_NUM,ORDER_ND), the indices of the points in
%    the grid.  The second dimension of this array is equal to the
%    product of the entries of ORDER_1D.
%
  a = [];
  more = 0;
  p = 0;

  while ( 1 )

    [ a, more ] = vec_colex_next2 ( dim_num, order_1d, a, more );

    if ( ~more )
      break
    end

    p = p + 1;
%
%  The values of A(DIM) are between 0 and ORDER_1D(DIM)-1 = N - 1 = 2 * M.
%  Subtracting M sets the range to -M to +M, as we wish.
%
    indx(1:dim_num,p) = a(1:dim_num) - round ( ( order_1d(1:dim_num) - 1 ) / 2 );

  end

  return
end
function w_nd = product_weight_gl ( dim_num, order_1d, order_nd )

%*****************************************************************************80
%
%% PRODUCT_WEIGHT_GL: weights for a product Gauss-Legendre rule.
%
%  Discussion:
%
%    This routine computes the weights for a quadrature rule which is
%    a product of 1D Gauss-Legendre rules of varying order.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    26 September 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer ORDER_1D(DIM_NUM), the order of the 1D rules.
%
%    Input, integer ORDER_ND, the order of the product rule.
%
%    Output, real W_ND(ORDER_ND), the product rule weights.
%
  w_nd(1:order_nd) = 1.0;

  for dim = 1 : dim_num

    w_1d = gl_weights ( order_1d(dim) );

    w_nd = r8vec_direct_product2 ( dim, order_1d(dim), w_1d, dim_num, ...
      order_nd, w_nd );
 
  end

  return
end
function r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

%*****************************************************************************80
%
%% R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns.
%
%    Input, real A(M,N), an M by N matrix to be printed.
%
%    Input, integer ILO, JLO, the first row and column to print.
%
%    Input, integer IHI, JHI, the last row and column to print.
%
%    Input, string TITLE, an optional title.
%
  incx = 5;

  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  for i2lo = max ( ilo, 1 ) : incx : min ( ihi, m )

    i2hi = i2lo + incx - 1;
    i2hi = min ( i2hi, m );
    i2hi = min ( i2hi, ihi );

    inc = i2hi + 1 - i2lo;
    
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Row: ' );
    for i = i2lo : i2hi
      fprintf ( 1, '%7d       ', i );
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Col\n' );

    j2lo = max ( jlo, 1 );
    j2hi = min ( jhi, n );

    for j = j2lo : j2hi

      fprintf ( 1, '%5d ', j );
      for i2 = 1 : inc
        i = i2lo - 1 + i2;
        fprintf ( 1, '%12f', a(i,j) );
      end
      fprintf ( 1, '\n' );

    end

  end

  return
end
function r8mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% R8MAT_WRITE writes an R8MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, real TABLE(M,N), the points.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 ) 
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'R8MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
%  For greater precision, try:
%
%     fprintf ( output_unit, '  %24,16f', table(i,j) );
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %14f', table(i,j) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end
function w = r8vec_direct_product2 ( factor_index, factor_order, ...
  factor_value, factor_num, point_num, w )

%*****************************************************************************80
%
%% R8VEC_DIRECT_PRODUCT2 creates a direct product of R8VEC's.
%
%  Discussion:
%
%    To explain what is going on here, suppose we had to construct
%    a multidimensional quadrature rule as the product of K rules
%    for 1D quadrature.
%
%    The product rule will be represented as a list of points and weights.
%
%    The J-th item in the product rule will be associated with
%      item J1 of 1D rule 1,
%      item J2 of 1D rule 2, 
%      ..., 
%      item JK of 1D rule K.
%
%    In particular, 
%      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
%    and
%      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
%
%    So we can construct the quadrature rule if we can properly
%    distribute the information in the 1D quadrature rules.
%
%    This routine carries out that task for the weights W.
%
%    Another way to do this would be to compute, one by one, the
%    set of all possible indices (J1,J2,...,JK), and then index
%    the appropriate information.  An advantage of the method shown
%    here is that you can process the K-th set of information and
%    then discard it.
%
%  Example:
%
%    Rule 1: 
%      Order = 4
%      W(1:4) = ( 2, 3, 5, 7 )
%
%    Rule 2:
%      Order = 3
%      W(1:3) = ( 11, 13, 17 )
%
%    Rule 3:
%      Order = 2
%      W(1:2) = ( 19, 23 )
%
%    Product Rule:
%      Order = 24
%      W(1:24) =
%        ( 2 * 11 * 19 )
%        ( 3 * 11 * 19 )
%        ( 4 * 11 * 19 )
%        ( 7 * 11 * 19 )
%        ( 2 * 13 * 19 )
%        ( 3 * 13 * 19 )
%        ( 5 * 13 * 19 )
%        ( 7 * 13 * 19 )
%        ( 2 * 17 * 19 )
%        ( 3 * 17 * 19 )
%        ( 5 * 17 * 19 )
%        ( 7 * 17 * 19 )
%        ( 2 * 11 * 23 )
%        ( 3 * 11 * 23 )
%        ( 5 * 11 * 23 )
%        ( 7 * 11 * 23 )
%        ( 2 * 13 * 23 )
%        ( 3 * 13 * 23 )
%        ( 5 * 13 * 23 )
%        ( 7 * 13 * 23 )
%        ( 2 * 17 * 23 )
%        ( 3 * 17 * 23 )
%        ( 5 * 17 * 23 )
%        ( 7 * 17 * 23 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR_INDEX, the index of the factor being processed.
%    The first factor processed must be factor 1.
%
%    Input, integer FACTOR_ORDER, the order of the factor.
%
%    Input, real FACTOR_VALUE(FACTOR_ORDER), the factor values for
%    factor FACTOR_INDEX.
%
%    Input, integer FACTOR_NUM, the number of factors.
%
%    Input, integer POINT_NUM, the number of elements in the direct product.
%
%    Input/output, real W(POINT_NUM), the elements of the
%    direct product, updated by the latest factor.
%
%  Local Parameters:
%
%    Local, integer START, the first location of a block of values to set.
%
%    Local, integer CONTIG, the number of consecutive values to set.
%
%    Local, integer SKIP, the distance from the current value of START
%    to the next location of a block of values to set.
%
%    Local, integer REP, the number of blocks of values to set.
%
  persistent contig;
  persistent rep;
  persistent skip;

  if ( factor_index == 1 )
    contig = 1;
    skip = 1;
    rep = point_num;
    w(1:point_num) = 1.0;
  end

  rep = rep / factor_order;
  skip = skip * factor_order;

  for j = 1 : factor_order

    start = 1 + ( j - 1 ) * contig;

    for k = 1 : rep
      w(start:start+contig-1) = w(start:start+contig-1) * factor_value(j);
      start = start + skip;
    end

  end

  contig = contig * factor_order;

  return
end
function r8vec_print_some ( n, a, i_lo, i_hi, title )

%*****************************************************************************80
%
%% R8VEC_PRINT_SOME prints "some" of an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    16 October 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the dimension of the vector.
%
%    Input, real A(N), the vector to be printed.
%
%    Input, integer MAX_PRINT, the maximum number of lines to print.
%
%    Input, string TITLE, an optional title.
%
  if ( 0 < s_len_trim ( title ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, '%s\n', title );
  end

  fprintf ( 1, '\n' );

  for i = max ( 1, i_lo ) : min ( n, i_hi )
    fprintf ( 1, '  %8d  %12f\n', i, a(i) );
  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LENGTH, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function [ grid_weight, grid_point ] = sparse_grid_gl ( dim_num, level_max, ...
  point_num )
  
%*****************************************************************************80
%
%% SPARSE_GRID_GL computes a sparse grid of Gauss-Legendre points. 
% 
%  Discussion: 
% 
%    The quadrature rule is associated with a sparse grid derived from 
%    a Smolyak construction using a 1D Gauss-Legendre quadrature rule.  
% 
%    The user specifies: 
%    * the spatial dimension of the quadrature region, 
%    * the level that defines the Smolyak grid. 
% 
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    04 July 2008
% 
%  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 
%    sparse grid. 
% 
%    Input, integer POINT_NUM, the number of points in the grid, 
%    as determined by SPARSE_GRID_GL_SIZE. 
% 
%    Output, real GRID_WEIGHT(POINT_NUM), the weights. 
% 
%    Output, real GRID_POINT(DIM_NUM,POINT_NUM), the points. 
% 
  grid_weight(1:point_num) = 0.0;
%
%  The outer loop generates LEVELs from LEVEL_MIN to LEVEL_MAX.
%
  point_num2 = 0;

  level_min = max ( 0, level_max + 1 - dim_num );
  
  for level = level_min : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = 0;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%  The relationship is the same as for other OPEN rules.
%  The GL rule differs from the other OPEN rules only in the nesting behavior.
%
      order_1d = level_to_order_open ( dim_num, level_1d );

      grid_base2(1:dim_num) = round ( ( order_1d(1:dim_num) - 1 ) / 2 );
%
%  The product of the 1D orders gives us the number of points in this subgrid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  Compute the weights for this product grid.
%
      grid_weight2 = product_weight_gl ( dim_num, order_1d, order_nd );
%
%  Now determine the coefficient of the weight.
%
      coeff = (-1)^( level_max - level ) ...
        * i4_choose ( dim_num - 1, level_max - level );
%
%  The inner (hidden) loop generates all points corresponding to given grid.
%  The grid indices will be between -M to +M, where 2*M + 1 = ORDER_1D(DIM).
%
      grid_index2 = multigrid_index_z ( dim_num, order_1d, order_nd );
%
%  Determine the first level of appearance of each of the points.
%  This allows us to flag certain points as being repeats of points
%  generated on a grid of lower level.  
%
%  This is SLIGHTLY tricky.
%
      grid_level = index_level_gl ( level, level_max, dim_num, order_nd, ...
        grid_index2, grid_base2 );
%
%  Only keep those points which first appear on this level.
%
      for point = 1 : order_nd
%
%  Either a "new" point (increase count, create point, create weight)
%
        if ( grid_level(point) == level )

          point_num2 = point_num2 + 1;

          grid_point(1:dim_num,point_num2) = gl_abscissa ( dim_num, 1, ...
            grid_index2(1:dim_num,point), grid_base2(1:dim_num) );

          grid_weight(point_num2) = coeff * grid_weight2(point);
%
%  or an already existing point (create point temporarily, find match,
%  add weight to matched point's weight).
%
        else
 
          grid_point_temp(1:dim_num) = gl_abscissa ( dim_num, 1, ...
            grid_index2(1:dim_num,point), grid_base2(1:dim_num) );
            
          point3 = -1;
          
          for point2 = 1 : point_num2
            if ( grid_point(1:dim_num,point2) == grid_point_temp(1:dim_num)' )
              point3 = point2;
              break
            end
          end
          
          if ( point3 == -1 )
            fprintf ( 1, '\n' );
            fprintf ( 1, 'SPARSE_GRID_GL - Fatal error!\n' );
            fprintf ( 1, '  Could not match point.\n' );
            error ( 'SPARSE_GRID_GL - Fatal error!' );
          end
          
          grid_weight(point3) = grid_weight(point3) + coeff * grid_weight2(point);
                   
        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end
function [ grid_index, grid_base ] = sparse_grid_gl_index ( dim_num, ...
  level_max, point_num )

%*****************************************************************************80
%
%% SPARSE_GRID_GL_INDEX indexes the points forming a sparse grid of GL points.
%
%  Discussion:
%
%    The sparse grid is assumed to be formed from 1D Gauss-Legendre rules
%    of ODD order, which have the property that only the central abscissa,
%    X = 0.0, is "nested".
%
%    The necessary dimensions of GRID_INDEX can be determined by 
%    calling SPARSE_GRID_GL_SIZE first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    04 July 2008
%
%  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 POINT_NUM, the total number of points in 
%    the grids.
%
%    Output, integer GRID_INDEX(DIM_NUM,POINT_NUM), a list of 
%    point indices, representing a subset of the product grid of level 
%    LEVEL_MAX, representing (exactly once) each point that will show up in a
%    sparse grid of level LEVEL_MAX.
%
%    Output, integer GRID_BASE(DIM_NUM,POINT_NUM), a list of 
%    the orders of the Gauss-Legendre rules associated with each point 
%    and dimension.
%

%
%  The outer loop generates LEVELs from LEVEL_MIN to LEVEL_MAX.
%
  point_num2 = 0;

  level_min = max ( 0, level_max + 1 - dim_num );
  
  for level = level_min : level_max
%
%  The middle loop generates the next partition LEVEL_1D(1:DIM_NUM)
%  that adds up to LEVEL.
%
    level_1d = [];
    more = false;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%  The relationship is the same as for other OPEN rules.
%  The GL rule differs from the other OPEN rules only in the nesting behavior.
%
      order_1d = level_to_order_open ( dim_num, level_1d );

      grid_base2(1:dim_num) = round ( ( order_1d(1:dim_num) - 1 ) / 2 );
%
%  The product of the 1D orders gives us the number of points in this subgrid.
%
      order_nd = prod ( order_1d(1:dim_num) );
%
%  The inner (hidden) loop generates all points corresponding to given grid.
%  The grid indices will be between -M to +M, where 2*M + 1 = ORDER_1D(DIM).
%
      grid_index2 = multigrid_index_z ( dim_num, order_1d, order_nd );
%
%  Determine the first level of appearance of each of the points.
%  This allows us to flag certain points as being repeats of points
%  generated on a grid of lower level.  
%
%  This is SLIGHTLY tricky.
%
      grid_level = index_level_gl ( level, level_max, dim_num, order_nd, ...
        grid_index2, grid_base2 );
%
%  Only keep those points which first appear on this level.
%
      for point = 1 : order_nd

        if ( grid_level(point) == level )

          point_num2 = point_num2 + 1;

          grid_index(1:dim_num,point_num2) = grid_index2(1:dim_num,point);
          grid_base(1:dim_num,point_num2) = grid_base2(1:dim_num);

        end

      end

      if ( ~more )
        break
      end

    end

  end

  return
end
function point_num = sparse_grid_gl_size ( dim_num, level_max )

%*****************************************************************************80
%
%% SPARSE_GRID_GL_SIZE sizes a sparse grid of Gauss-Legendre points.
%
%  Discussion:
%
%    The grid is defined as the sum of the product rules whose LEVEL
%    satisfies:
%
%      LEVEL_MIN <= LEVEL <= LEVEL_MAX.
%
%    where LEVEL_MAX is user specified, and 
%
%      LEVEL_MIN = max ( 0, LEVEL_MAX + 1 - DIM_NUM ).
%
%    The grids are only very weakly nested, since Gauss-Legendre rules
%    only have the origin in common.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    04 July 2008
%
%  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.
%
%    Output, integer POINT_NUM, the number of points in the grid.
%

%
%  Special case.
%
  if ( level_max == 0 )
    point_num = 1;
    return;
  end
%
%  The outer loop generates LEVELs from LEVEL_MIN to LEVEL_MAX.
%
  point_num = 0;

  level_min = max ( 0, level_max + 1 - dim_num );
  
  for level = level_min : level_max
%
%  The middle loop generates the next partition that adds up to LEVEL.
%
    level_1d = [];
    more = false;
    h = 0;
    t = 0;

    while ( 1 )

      [ level_1d, more, h, t ] = comp_next ( level, dim_num, level_1d, more, h, t );
%
%  Transform each 1D level to a corresponding 1D order.
%
      order_1d = level_to_order_open ( dim_num, level_1d );

      for dim = 1 : dim_num
%
%  If we can reduce the level in this dimension by 1 and
%  still not go below LEVEL_MIN.
%
        if ( level_min < level & 1 < order_1d(dim) )
          order_1d(dim) = order_1d(dim) - 1;
        end
		
      end

      point_num = point_num + prod ( order_1d(1:dim_num) );

      if ( ~more )
        break
      end
	  
    end
  end

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
function [ a, more ] = vec_colex_next2 ( dim_num, base, a, more )

%*****************************************************************************80
%
%% VEC_COLEX_NEXT2 generates vectors in colex order.
%
%  Discussion:
%
%    The vectors are produced in colexical order, starting with
%    (0,0,...,0),
%    (1,0,...,0),
%    ...
%    (BASE(1)-1,BASE(2)-1,...,BASE(DIM_NUM)-1).
%
%  Example:
%
%    DIM_NUM = 2, 
%    BASE = [ 3, 3]
%
%    0   0
%    1   0
%    2   0
%    0   1
%    1   1
%    2   1
%    0   2
%    1   2
%    2   2
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dennis Stanton, Dennis White,
%    Constructive Combinatorics,
%    Springer, 1986,
%    ISBN: 0387963472,
%    LC: QA164.S79.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension.
%
%    Input, integer BASE(DIM_NUM), the base to be used in each dimension.
%
%    Input, integer A(DIM_NUM), except on the first call, this should
%    be the output value of A on the last call.
%
%    Input, logical MORE, should be FALSE on the first call, and
%    thereafter should be the output value of MORE from the previous call.  
%
%    Output, integer A(DIM_NUM), the next vector.
%
%    Output, logical MORE, is TRUE if another vector was computed.
%    If MORE is FALSE on return, then ignore the output value A, and
%    stop calling the routine.
%
  if ( ~more )

    a(1:dim_num) = 0;
    more = 1;

  else
      
    for i = 1 : dim_num

      a(i) = a(i) + 1;

      if ( a(i) < base(i) )
        return
      end

      a(i) = 0;

    end

    more = 0;

  end

  return
end