function svd_basis ( DUMMY )

%% MAIN is the main program for SVD_BASIS.
%
%  Discussion:
%
%    SVD_BASIS forms a basis from the SVD of a set of data vectors.
%
%    This program uses the singular value decomposition (SVD) to analyze
%    a set of data, and extract a number of dominant modes.
%
%    This program is intended as an intermediate application, in
%    the following situation:
%
%    A) a "high fidelity" or "high resolution" PDE solver is used
%       to determine many (say N = 500) solutions of a discretized
%       PDE at various times, or parameter values.  Each solution
%       may be regarded as an M vector.  Typically, each solution
%       involves an M by M linear system, greatly reduced in
%       complexity because of bandedness or sparsity.
%
%    B) This program is applied to extract L dominant modes from
%       the N solutions.  This is done using the singular value
%       decomposition of the M by N matrix, each of whose columns
%       is one of the original solution vectors.
%
%    C) a "reduced order model" program may then attempt to solve
%       a discretized version of the PDE, using the L dominant
%       modes as basis vectors.  Typically, this means that a dense
%       L by L linear system will be involved.
%
%    Thus, the program might read in 500 files, and write out
%    5 or 10 files of the corresponding size and "shape", representing
%    the dominant solution modes.
%
%    An option has been added to compute the average of the vectors
%    and subtract it before SVD processing.
%
%  Modified:
%
%    15 May 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gal Berkooz, Philip Holmes, John Lumley,
%    The proper orthogonal decomposition in the analysis of turbulent flows,
%    Annual Review of Fluid Mechanics,
%    Volume 25, 1993, pages 539-575.
%
%    John Burkardt, Max Gunzburger, Hyung-Chun Lee,
%    Centroidal Voronoi Tessellation-Based Reduced-Order
%    Modelling of Complex Systems,
%    SIAM Journal on Scientific Computing,
%    Volume 28, Number 2, 2006, pages 459-484.
%
%    Lawrence Sirovich,
%    Turbulence and the dynamics of coherent structures, Parts I-III,
%    Quarterly of Applied Mathematics,
%    Volume XLV, Number 3, 1987, pages 561-590.
%
  
  clean = 1;

  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_BASIS\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Given a PDE for which:\n' );
  fprintf ( 1, '    C is the number of components of the solution \n' );
  fprintf ( 1, '      at any single point,\n' );
  fprintf ( 1, '    P is the number of points where a solution is given,\n' );
  fprintf ( 1, '    N is the number of solution vectors,\n' );
  fprintf ( 1, '    L is the number of modes to be extracted.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Then we let M = C*P be the abstract spatial dimension.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  If requested, we compute the average solution,\n' );
  fprintf ( 1, '  subtract it from each solution, and save that\n' );
  fprintf ( 1, '  as mode #0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Set up A, the M by N matrix of solution vectors,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Get A = U * S * V'', the singular value decomposition.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The first L columns of U are our dominant modes.\n' );
  fprintf ( 1, '\n' );
%
%  What is the basis size?
%
  basis_num = input ( '  How many basis vectors (L) are to be extracted?' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L = %d\n', basis_num );
%
%  Gather one or more "base" file names.
%
  data_file_base_num = 0;
  data_file_base_start = [];
  data_file_base_last = [];
  data_file_base = [];
  data_file_base_size = 0;

  while ( 1 )
%
%  Get the next base file name.
%
    fprintf ( 1, '\n' );
    fprintf ( 1, '  You specify a consecutive sequence of file names\n' );
    fprintf ( 1, '  by giving the first "base" file name.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  If there are no more sequences to enter,\n' );
    fprintf ( 1, '  just hit RETURN.\n' );

    file_name = input ( '  Enter a new base file name, or RETURN.' );

    if ( s_len_trim ( file_name ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  RETURN was entered.\n' );
      fprintf ( 1, '  Presumably, there are no more file sequences.\n' );
      break
    end

    file_name_length = length ( file_name );
    
    data_file_base_num = data_file_base_num + 1;
    data_file_base_start(data_file_base_num) = data_file_base_size + 1;
    data_file_base_last(data_file_base_num) = data_file_base_size + file_name_length;
    data_file_base = strcat ( data_file_base, file_name );
    data_file_base_size = data_file_base_size + file_name_length;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %6d:  "%s"\n', data_file_base_num, file_name );
%
%  For the very first base file, get the data sizes.
%
    if ( data_file_base_num == 1 )

      [ comp_num, node_num ] = dtable_header_read ( file_name );

      dim_num = comp_num * node_num;

      fprintf ( 1, '\n' );
      fprintf ( 1, '  According to the first base file,\n' );
      fprintf ( 1, '  The number of solution components C =   %d\n', comp_num );
      fprintf ( 1, '  The number of solution points P =       %d\n', node_num );
      fprintf ( 1, '  The "size" of each solution M = (C*P) = %d\n', dim_num );
%
%  Idiocy check.  L must be less than or equal to M.
%
      if ( dim_num < basis_num )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SVD_BASIS - Fatal error!\n' );
        fprintf ( 1, '\n' );
        fprintf ( 1, '  M < L.\n' );
        fprintf ( 1, '\n' );
        fprintf ( 1, '  That is, the number of modes requested (L) is greater\n' );
        fprintf ( 1, '  than the spatial dimension (M).\n' );
        fprintf ( 1, '  Technically, the program could pad out the answer\n' );
        fprintf ( 1, '  with L-M zero vectors, but instead, we will stop\n' );
        fprintf ( 1, '  assuming you made an error, or a misapprehension.\n' );
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SVD_BASIS:\n' );
        fprintf ( 1, '  Abnormal end of execution.\n' );
        fprintf ( 1, '\n' );
        timestamp;
        return
      end

    end

  end
%
%  Count all the data files.
%
  data_file_num = 0;

  for i = 1 : data_file_base_num

    start = data_file_base_start(i);
    last = data_file_base_last(i);
    data_file = data_file_base(start:last);

    while ( 1 )
      
      if ( ~file_exist ( data_file ) )
        break
      end

      data_file_num = data_file_num + 1;

      data_file = file_name_inc ( data_file );

    end

  end

  if ( data_file_num == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SVD_BASIS - Fatal error!\n' );
    fprintf ( 1, '  There do not seem to be any solution files;\n' );
    fprintf ( 1, '  that is, files whose names are "incremented"\n' );
    fprintf ( 1, '  versions of the first file name.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The first file we looked for was "%s"\n', data_file );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SVD_BASIS:\n' );
    fprintf ( 1, '  Abnormal end of execution.\n' );
    fprintf ( 1, '\n' );
    timestamp;
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data files N = %d\n', data_file_num );
%
%  Set up an array to hold all the data.
%
  point_num = data_file_num;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data is stored in an M by N matrix A.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The "spatial" dimension M is   %d\n', dim_num );
  fprintf ( 1, '  The number of data points N is %d\n', point_num );
%
%  Read the data.
%
  l = 0;

  for ii = 1 : data_file_base_num

    start = data_file_base_start(i);
    last = data_file_base_last(i);
    data_file = data_file_base(start:last);

    while ( 1 )

      if ( ~file_exist ( data_file ) )
        break
      end

      l = l + 1;

      table = dtable_data_read ( data_file, comp_num, node_num );

      k = 0;
      for j = 1 : node_num
        for i = 1 : comp_num
          k = k + 1;
          point(k,l) = table(i,j);
        end
      end

      data_file = file_name_inc ( data_file );

    end

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data has been read into the matrix A.\n' );
%
%----------------------------------------------------------------------------
%  Optionally, average the data, subtract the average from each entry,
%  and later save the average as vector #0.
%----------------------------------------------------------------------------
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_BASIS:\n' );
  fprintf ( 1, '  Averaging of data is optional.\n' );
  fprintf ( 1, '  The program can average the data vectors,\n' );
  fprintf ( 1, '  subtract it from each data vector,\n' );
  fprintf ( 1, '  and write out the data average vector as an\n' );
  fprintf ( 1, '  extra "mode 0" vector.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Do you want to compute and use the average? (''Y''/''N'')\n' );

  average_char = input ( '  Enter ''Y'' or ''N'':' );

  if ( average_char == 'Y' | average_char == 'y' )
    average_normalization = 1;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The user has requested the average vector.\n' );
  else
    average_normalization = 0;
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The user does not want the average vector.\n' );
  end

  if ( average_normalization )

    for i = 1 : dim_num
      point_average(i,1) = sum ( point(i,1:point_num) );
    end

    point_average(1:dim_num,1) = point_average(1:dim_num,1) / point_num;

    for i = 1 : dim_num
      point(i,1:point_num) = point(i,1:point_num) - point_average(i,1);
    end

  end
%
%----------------------------------------------------------------------------
%
%  Compute the SVD of A.
%
%----------------------------------------------------------------------------
%
  [ point, sval ] = singular_vectors ( dim_num, point_num, basis_num, point );
%
%----------------------------------------------------------------------------
%
%  "Clean" the output data.  We are having problems with some vectors
%  containing a few very tiny (and meaningless) values.
%
%----------------------------------------------------------------------------
%
  if ( clean )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Because the CLEAN option is on,\n' );
    fprintf ( 1, '  we will set very tiny vector entries to 0.\n' );

    tol = r8_epsilon ( 'DUMMY' );

    for j = 1 : basis_num
      for i = 1 : dim_num
        if ( abs ( point(i,j) ) < tol )
          point(i,j) = 0.0;
        end
      end
    end
  end
%
%----------------------------------------------------------------------------
%
%  Write the first L left singular vectors (columns of U) to files.
%
%----------------------------------------------------------------------------
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_BASIS:\n' );
  fprintf ( 1, '  Ready to write the left singular vectors to files.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Do you want comments in the header of the file?\n' );
  fprintf ( 1, '  (These begin with the "#" character.) (''Y''/''N'')\n' );

  comment_char = input ( '  Enter ''Y'' or ''N'':' );

  if ( comment_char == 'Y' | comment_char == 'y' )
    comment = 1;
  else
    comment = 0;
  end

  basis_file = 'svd_000.txt';

  if ( average_normalization )

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Writing average file "%s".\n', basis_file );

    average_value = 0.0;
    
    basis_write ( basis_file, comp_num, node_num, average_value, ...
      point_average(1:dim_num,1), comment );

  end

  for j = 1 : basis_num

    basis_file = file_name_inc ( basis_file );

    if ( j == 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Writing first file "%s".\n', basis_file );
    end

    if ( j == basis_num )
      fprintf ( 1, '  Writing last file  "%s"\n', basis_file );
    end

    basis_write ( basis_file, comp_num, node_num, sval(j), ...
      point(1:dim_num,j), comment );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SVD_BASIS\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp;