function mgmres_test02 ( )

%*****************************************************************************80
%
%% mgmres_test02() tests mgmres_st() on a 9 by 9 matrix.
%
%  Discussion:
%
%    A =
%      2  0  0 -1  0  0  0  0  0
%      0  2 -1  0  0  0  0  0  0
%      0 -1  2  0  0  0  0  0  0
%     -1  0  0  2 -1  0  0  0  0
%      0  0  0 -1  2 -1  0  0  0
%      0  0  0  0 -1  2 -1  0  0
%      0  0  0  0  0 -1  2 -1  0
%      0  0  0  0  0  0 -1  2 -1
%      0  0  0  0  0  0  0 -1  2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 July 2007
%
%  Author:
%
%    John Burkardt
%
  n = 9;
  nz_num = 23;

  a = [ ...
    2.0, -1.0, ...
    2.0, -1.0, ...
    -1.0, 2.0, ...
    -1.0, 2.0, -1.0, ...
    -1.0, 2.0, -1.0, ...
    -1.0, 2.0, -1.0, ...
    -1.0, 2.0, -1.0, ...
    -1.0, 2.0, -1.0, ...
    -1.0, 2.0 ];
  ia = [ ...
    1, 1, ...
    2, 2, ...
    3, 3, ...
    4, 4, 4, ...
    5, 5, 5, ...
    6, 6, 6, ...
    7, 7, 7, ...
    8, 8, 8, ...
    9, 9 ];
  ja = [ ...
    1, 4, ...
    2, 3, ...
    2, 3, ...
    1, 4, 5, ...
    4, 5, 6, ...
    5, 6, 7, ...
    6, 7, 8, ...
    7, 8, 9, ...
    8, 9 ];
  rhs = [ ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0, ...
    1.0 ];
  x_exact = [ ...
    3.5, ...
    1.0, ...
    1.0, ...
    6.0, ...
    7.5, ...
    8.0, ...
    7.5, ...
    6.0, ...
    3.5 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'mgmres_test02():\n' );
  fprintf ( 1, '  Test mgmres_st() on matrix that is not quite \n' );
  fprintf ( 1, '  the -1,2,-1 matrix, of order N = %d\n', n );

  for test = 1 : 2

    if ( test == 1 )

      fprintf ( 1, '\n' );
      fprintf ( 1, '  First try, set initial estimate zero.\n' );
      fprintf ( 1, '\n' );

      x_estimate(1:n) = 0;

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Second try, set initial estimate random.\n' );
      fprintf ( 1, '\n' );
      x_estimate = rand ( n, 1 );

    end

    itr_max = 20;
    mr = n - 1;
    tol_abs = 1.0E-08;
    tol_rel = 1.0E-08;

    x_error = sqrt ( sum ( ( x_exact(1:n) - x_estimate(1:n) ).^2 ) );

    fprintf ( 1, '  Before solving, X_ERROR = %f\n', x_error );

    x_estimate = mgmres_st ( n, nz_num, ia, ja, a, x_estimate, rhs, ...
      itr_max, mr, tol_abs, tol_rel );

    x_error = sqrt ( sum ( ( x_exact(1:n) - x_estimate(1:n) ).^2 ) );

    fprintf ( 1, '  After solving, X_ERROR = %f\n', x_error );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Final solution estimate:\n' );
    fprintf ( 1, '\n' );
    for i = 1 : n
      fprintf ( 1, '  %8d  %14f\n', i, x_estimate(i) );
    end

  end

  return
end
