function x = r83v_jac_sl ( n, a, b, c, ax, x, it_max )

%*****************************************************************************80
%
%% R83V_JAC_SL solves a R83V system A*x=b using Jacobi iteration.
%
%  Discussion:
%
%    The R83V storage format is used for a tridiagonal matrix.
%    The subdiagonal is in A(min(M-1,N)).
%    The diagonal is in B(min(M,N)).
%    The superdiagonal is in C(min(M,N-1)).
%
%  Example:
%
%    An R83V matrix of order 3x5 would be stored:
%
%      B1  C1  **  **  **
%      A1  B2  C2  **  **
%      **  A2  B3  C3  **
%
%    An R83 matrix of order 5x5 would be stored:
%
%      B1  C1  **  **  **
%      A1  B2  C2  **  **
%      **  A2  B3  C3  **
%      **  **  A3  B4  C4
%      **  **  **  A4  B5
%
%    An R83 matrix of order 5x3 would be stored:
%
%      B1  C1  **
%      A1  B2  C2
%      **  A2  B3
%      **  **  A3
%      **  **  **
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 February 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be at least 2.
%
%    Input, real A(min(M-1,N)), B(min(M,N)), C(min(M,N-1)), the R83V matrix.
%
%    Input, real AX(N), the right hand side of the linear system.
%
%    Input, real X(N), an approximate solution to the system.
%
%    Input, integer IT_MAX, the maximum number of iterations.
%
%    Output, real X(N), an updated approximate solution to the system.
%
  ax = ax(:);
  x = x(:); 
%
%  No diagonal matrix entry can be zero.
%
  for i = 1 : n
    if ( b(i) == 0.0E+00 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R83V_JAC_SL - Fatal error!\n' );
      fprintf ( 1, '  Zero diagonal entry, index = %d\n', i );
      error ( 'R83V_JAC_SL - Fatal error!' );
    end
  end
%
%  Iterate IT_MAX times.
%
  x_new = zeros ( n, 1 );

  for it_num = 1 : it_max

    x_new(1:n)   = ax(1:n);
    x_new(1:n-1) = x_new(1:n-1) - c(1:n-1) .* x(2:n);
    x_new(2:n)   = x_new(2:n)   - a(1:n-1) .* x(1:n-1);
%
%  Divide by diagonal terms.
%
    x_new(1:n) = x_new(1:n) ./ b(1:n);
%
%  Update.
%
    x(1:n) = x_new(1:n);

  end

  return
end

