function x = r83v_fs ( n, a1, a2, a3, b )

%*****************************************************************************80
%
%% R83V_FS solves a linear system with R83V matrix.
%
%  Discussion:
%
%    This function is based on the LINPACK SGTSL routine.
%
%    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:
%
%    20 February 2016
%
%  Author:
%
%    John Burkardt, based on the LINPACK SGTSL function.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, integer N, the order of the tridiagonal matrix.
%
%    Input, real A1(N-1), A2(N), A3(N-1), the R83V matrix.
%
%    Input, real B(N), the right hand side.
%
%    Output, real X(N), the solution.
%

%
%  Copy the input data.
%
  c(1) = 0.0;
  c(2:n) = a1(1:n-1);
  d(1:n) = a2(1:n);
  e(1:n-1) = a3(1:n-1);
  e(n) = 0.0;

  x(1:n) = b(1:n);
%
%  Factor.
%
  c(1) = a2(1);

  if ( 2 <= n )

    d(1) = e(1);
    e(1) = 0.0;
    e(n) = 0.0;

    for k = 2 : n
%
%  Find the larger of the two rows and interchange if necessary.
%
      if ( abs ( c(k-1) ) <= abs ( c(k) ) )

        t = c(k);
        c(k) = c(k-1);
        c(k-1) = t;

        t = d(k);
        d(k) = d(k-1);
        d(k-1) = t;

        t = e(k);
        e(k) = e(k-1);
        e(k-1) = t;

        t = x(k);
        x(k) = x(k-1);
        x(k-1) = t;

      end
%
%  Zero elements.
%
      if ( c(k) == 0.0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'R83V_FS - Fatal error!\n' )
        fprintf ( 1, '  Zero pivot on step K = %d\n', k );
        error ( 'R83V_FS - Fatal error!' );
      end

      t = - c(k) / c(k-1);
      c(k) = d(k) + t * d(k-1);
      d(k) = e(k) + t * e(k-1);
      e(k) = 0.0;
      x(k) = x(k) + t * x(k-1);

    end

  end

  if ( c(n) == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83V_FS - Fatal error!\n' )
    fprintf ( 1, '  Zero pivot on step K = %d\n', n );
    error ( 'R83V_FS - Fatal error!' );
  end
%
%  Back solve.
%
  x(n) = x(n) / c(n);

  if ( 1 < n )

    x(n-1) = ( x(n-1) - d(n-1) * x(n) ) / c(n-1);

    for k = n - 2 : -1 : 1
      x(k) = ( x(k) - d(k) * x(k+1) - e(k) * x(k+2) ) / c(k);
    end

  end

  return
end