function x = r83v_cg ( n, a, b, c, ax, x )

%*****************************************************************************80
%
%% R83V_CG uses the conjugate gradient method on an R83V system.
%
%  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
%
%  Reference:
%
%    Frank Beckman,
%    The Solution of Linear Equations by the Conjugate Gradient Method,
%    in Mathematical Methods for Digital Computers,
%    edited by John Ralston, Herbert Wilf,
%    Wiley, 1967,
%    ISBN: 0471706892,
%    LC: QA76.5.R3.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    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 vector.
%
%    Input/output, real X(N).
%    On input, an estimate for the solution, which may be 0.
%    On output, the approximate solution vector.
%
  ax = ax(:);
  x = x(:);
%
%  Initialize
%    AP = A * x,
%    R  = b - A * x,
%    P  = b - A * x.
%
  ap = r83v_mv ( n, n, a, b, c, x );

  r(1:n,1) = ax(1:n,1) - ap(1:n,1);
  p(1:n,1) = ax(1:n,1) - ap(1:n,1);
%
%  Do the N steps of the conjugate gradient method.
%
  for it = 1 : n
%
%  Compute the matrix*vector product AP=A*P.
%
    ap = r83v_mv ( n, n, a, b, c, p );
%
%  Compute the dot products
%    PAP = P*AP,
%    PR  = P*R
%  Set
%    ALPHA = PR / PAP.
%
    pap = p' * ap;
    pr = p' * r;

    if ( pap == 0.0 )
      return
    end

    alpha = pr / pap;
%
%  Set
%    X = X + ALPHA * P
%    R = R - ALPHA * AP.
%
    x(1:n,1) = x(1:n,1) + alpha * p(1:n,1);
    r(1:n,1) = r(1:n,1) - alpha * ap(1:n,1);
%
%  Compute the vector dot product
%    RAP = R*AP
%  Set
%    BETA = - RAP / PAP.
%
    rap = r' * ap;

    beta = - rap / pap;
%
%  Update the perturbation vector
%    P = R + BETA * P.
%
    p(1:n,1) = r(1:n,1) + beta * p(1:n,1);

  end

  return
end
