function product_factor ( list_filename, quad_filename )

%% MAIN is the main program for PRODUCT_FACTOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    John Burkardt
%
  timestamp;
  fprintf ( 1, '\n');
  fprintf ( 1, 'PRODUCT_FACTOR\n');
  fprintf ( 1, '  MATLAB version\n');
  fprintf ( 1, '\n');
  fprintf ( 1, '  Create a multidimensional product rule\n');
  fprintf ( 1, '  as a product of distinct 1D integration rules.\n');
%
%  Get the list filename.
%
  if ( nargin < 1 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_FACTOR:\n' );
    list_filename = input ( '  Enter the name of the file listing the factors.' );

  end
%
%  Get the product file prefix.
%
  if ( nargin < 2 )

    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_FACTOR:\n' );
    quad_filename = input ( '  Enter the product file prefix to use.' );

  end
%
%  Count the items in the list file.
%
  list_num = file_row_count ( list_filename );
%
%  Determine the spatial dimension and number of points in the product.
%
  dim_num = list_num;
  point_num = product_rule_size ( list_filename, list_num );
%
%  Allocate the product items.
%
  x(1:dim_num,1:point_num) = 0.0;
  w(1:point_num) = 1.0;
  r(1:dim_num,1:2) = 0.0;

  list_unit = fopen ( list_filename, 'rt' );

  if ( list_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
    fprintf ( 1, '  Nonzero value of IOS while opening list file.\n' );
    error ( 'PRODUCT_FACTOR - Fatal error!' );
  end
%
%  Read the factor information and apply it.
%
  for dim = 1 : dim_num

    quad_1d_filename = fgetl ( list_unit );

    quad_x_1d_filename = strcat ( quad_1d_filename, '_x.txt' );
    quad_w_1d_filename = strcat ( quad_1d_filename, '_w.txt' );
    quad_r_1d_filename = strcat ( quad_1d_filename, '_r.txt' );
%
%  Read the X file.
%
    [ dim_num_1d, point_num_1d ] = dtable_header_read ( quad_x_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature abscissa file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_FACTOR - Fatal error!' );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Number of points in 1D rule = %d\n', point_num_1d );

    x_1d = dtable_data_read ( quad_x_1d_filename, dim_num_1d, point_num_1d );
%
%  Read the W file.
%
    [ dim_num_1d, point_num_1d2 ] = dtable_header_read ( quad_w_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature weight file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_FACTOR - Fatal error!' );
    end

    if ( point_num_1d2 ~= point_num_1d )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature weight file should have exactly\n' );
      fprintf ( 1, '  the same number of lines as the abscissa file.\n' );
      error ( 'PRODUCT_FACTOR - Fatal error!' );
    end

    w_1d = dtable_data_read ( quad_w_1d_filename, dim_num_1d, point_num_1d );
%
%  Read the R file.
%
    [ dim_num_1d, point_num_1d2 ] = dtable_header_read ( quad_r_1d_filename );

    if ( dim_num_1d ~= 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature region file should have exactly\n' );
      fprintf ( 1, '  one value on each line.\n' );
      error ( 'PRODUCT_FACTOR - Fatal error!' );
    end

    if ( point_num_1d2 ~= 2 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'PRODUCT_FACTOR - Fatal error!\n' );
      fprintf ( 1, '  The 1D quadrature region file should have exactly\n' );
      fprintf ( 1, '  two lines.\n' );
      error ( 'PRODUCT_FACTOR - Fatal error!' );
    end

    r_1d = dtable_data_read ( quad_r_1d_filename, 1, 2 );
%
%  Update the X, W, and R of the product rule.
%
    x = r8vec_direct_product ( dim, point_num_1d, x_1d, dim_num, point_num, x );

    w = r8vec_direct_product2 ( dim, point_num_1d, w_1d, dim_num, point_num, w );

    r(dim,1) = r_1d(1);
    r(dim,2) = r_1d(2);

  end

  fclose ( list_unit );
%
%  Write the product rule.
%
  quad_filename = 'product';

  quad_x_filename = strcat ( quad_filename, '_x.txt' );
  quad_w_filename = strcat ( quad_filename, '_w.txt' );
  quad_r_filename = strcat ( quad_filename, '_r.txt' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule X file = "%s".\n', ...
    quad_x_filename );

  header = 0;
  dtable_write ( quad_x_filename, dim_num, point_num, x, header );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule W file = "%s".\n', ...
    quad_w_filename );

  header = 0;
  dtable_write ( quad_w_filename, 1, point_num, w, header );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Creating product quadrature rule R file = "%s".\n', ...
    quad_r_filename );

  header = 0;
  dtable_write ( quad_r_filename, dim_num, 2, r, header );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRODUCT_FACTOR:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp;