function x = gauss ( A, b, wait )

%*****************************************************************************80
%
%% gauss() uses Gauss elimination to solve a linear system A*x=b.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(m,n), the system matrix.
%
%    real b(m,k), the right hand side vectors.
%
%    option logical wait:
%    0: no extra printout (default)
%    1: extra printout, but no pause.
%    2: extra printout, and pausing.
%
%  Output:
%
%    global real x(m,k), the solution vectors.
%
  if ( nargin < 3 )
    wait = 0;
  end

  Ab = init ( A, b );

  [ m, n ] = size ( A );

  if ( 0 < wait )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Augmented matrix Ab, step 0\n' );
    fprintf ( 1, '\n' );
    disp ( Ab );
    if ( 1 < wait )
      pause ( );
    end
  end

  for j = 1 : n

    [ Ab, p ] = pivot ( j );

    Ab = scale ( j );

    for i = 1 : n
      if ( i ~= j )
        Ab = elim ( i, j );
      end
    end

    if ( 0 < wait )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Augmented matrix Ab, after step %d\n', j );
    fprintf ( 1, '\n' );
      disp ( Ab );
      if ( 1 < wait )
        pause ( );
      end
    end

  end

  x = solution ( );

  if ( 0 < wait )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Computed solution x\n' );
    fprintf ( 1, '\n' );
    disp ( x );
    if ( 1 < wait )
      pause ( );
    end
  end

  return
end
function Ab = elim ( i, j )

%*****************************************************************************80
%
%% elim uses Ab(j,j) to eliminate Ab(i,j).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer i, j, the row and column indices of the entry to eliminate.
%
%  Output:
%
%    real Ab(m,n+k), the augmented matrix has been modified.
%
  global Ab
  global k
  global m
  global n

  if ( i == j )
    fprintf ( 1, 'elim: Ignoring this command\n' );
    fprintf ( 1, '  I and J must be distinct.\n' );
    return
  end

  if ( Ab(i,j) == 0.0 )
    return
  end

  s = - Ab(i,j) / Ab(j,j);
  Ab(i,:) = Ab(i,:) + Ab(j,:) * s;
  Ab(i,j) = 0.0;

  return
end
function Ab = init ( A, b )

%*****************************************************************************80
%
%% init initializes the data.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real Ab(m,n+k), the augmented matrix has been set up.
%
  global Ab
  global k
  global m
  global n

  if ( nargin < 2 )
    b = [];
  end

  [ m, n ] = size ( A );
  [ m, k ] = size ( b );

  Ab = [ A, b ];

  return
end
function [ Ab, p ] = pivot ( j )

%*****************************************************************************80
%
%% pivot finds the pivot row for variable j.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer j, the variable for which a pivot is sought.
%
%  Output:
%
%    real Ab(m,n+k), the augmented matrix after pivoting.
%
%    integer p, the row that was chosen for the pivot.
%
  global Ab
  global k
  global m
  global n

  [ value, p ] = max ( abs ( Ab(j:n,j) ) );

  if ( value == 0.0 )
    fprintf ( 'pivot - Fatal error!\n' );
    fprintf ( '  No nonzero pivot could be found.\n' );
    error ( 'pivot - Fatal error!' );
  end

  p = p + j - 1;
  if ( p ~= j )
    t       = Ab(j,:);
    Ab(j,:) = Ab(p,:);
    Ab(p,:) = t;
  end

  return
end
function Ab = scale ( i, j )

%*****************************************************************************80
%
%% scale scales row i so that Ab(i,j) = 1.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer i, j, the row and column indices of the entry whose value
%    is used to normalize row i.  Ab(i,j) must not be zero.
%    j defaults to i.
%
%  Output:
%
%    real Ab(m,n+k), the augmented matrix has been modified.
%    Row i was divided by Ab(i,j).
%
  global Ab
  global k
  global m
  global n

  if ( nargin == 1 )
    j = i;
  end

  if ( i < 1 || m < i )
    fprintf ( 1, 'scale - Fatal error!\n' );
    fprintf ( 1, '  1 <= I <= %d is required.\n', m );
    error ( 'scale: Fatal error!' );
  end

  if ( j < 1 || n < j )
    fprintf ( 1, 'scale - Fatal error!\n' );
    fprintf ( 1, '  1 <= J <= %d is required.\n', n );
    error ( 'scale: Fatal error!' );
  end

  if ( Ab(i,j) == 0.0 )
    fprintf ( 1, 'scale - Fatal error!\n' );
    fprintf ( 1, '  Ab(%d,%d) = 0.\n', i, j );
    error ( 'scale: Fatal error!' );
  end

  Ab(i,:) = Ab(i,:) / Ab(i,j);

  return
end
function x = solution ( )

%*****************************************************************************80
%
%% solution returns the solution.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    03 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Output:
%
%    real x(m,k), the solution vectors.
%
  global Ab
  global k
  global m
  global n

  I = eye ( m, n );

  norm_fro = norm ( Ab(1:m,1:n) - I, 'fro' );

  if ( eps < norm_fro )
    fprintf ( 1, 'solution - Fatal error!\n' );
    fprintf ( 1, '  Augmented matrix Ab does not start with the identity.\n' );
    return
  end

  x = Ab(1:m,n+1:n+k);

  return
end

