function x = r83_cr_sls ( n, a_cr, nb, b )

%*****************************************************************************80
%
%% r83_cr_sls() solves several real linear systems factored by R83_CR_FA.
%
%  Discussion:
%
%    The R83 storage format is used for a tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:min(M+1,N)).
%    The diagonal in entries (2,1:min(M,N)).
%    The subdiagonal in (3,min(M-1,N)).
%    R8GE A(I,J) = R83 A(I-J+2,J).   
%
%    The matrix A must be tridiagonal.  R83_CR_FA is called to compute the
%    LU factors of A.  It does so using a form of cyclic reduction.  If
%    the factors computed by R83_CR_FA are passed to R83_CR_SLS, then one or many
%    linear systems involving the matrix A may be solved.
%
%    Note that R83_CR_FA does not perform pivoting, and so the solution 
%    produced by R83_CR_SLS may be less accurate than a solution produced 
%    by a standard Gauss algorithm.  However, such problems can be 
%    guaranteed not to occur if the matrix A is strictly diagonally 
%    dominant, that is, if the absolute value of the diagonal coefficient 
%    is greater than the sum of the absolute values of the two off diagonal 
%    coefficients, for each row of the matrix.
%
%  Example:
%
%    An R83 matrix of order 3x5 would be stored:
%
%       *  A12 A23 A34  *
%      A11 A22 A33  *   *
%      A21 A32  *   *   *
%
%    An R83 matrix of order 5x5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%    An R83 matrix of order 5x3 would be stored:
%
%       *  A12 A23
%      A11 A22 A33
%      A21 A32 A43
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 May 2010
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Hockney,
%    A fast direct solution of Poisson's equation using Fourier Analysis,
%    Journal of the ACM,
%    Volume 12, Number 1, pages 95-113, January 1965.
%
%  Input:
%
%    integer N, the order of the matrix.
%    N must be positive.
%
%    real A_CR(3,2*N+1), factorization information computed by R83_CR_FA.
%
%    integer NB, the number of right hand sides.
%
%    real B(N,NB), the right hand side vectors.
%
%  Output:
%
%    real X(N,NB), the solutions of the linear system.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83_CR_SLS - Fatal error!\n' );
    fprintf ( 1, '  Nonpositive N = %d\n', n );
    return
  end

  if ( n == 1 )
    x(1,1:nb) = a_cr(2,2) * b(1,1:nb);
    return
  end
%
%  Set up RHS.
%
  rhs(1,1:nb) = 0.0;
  rhs(2:n+1,1:nb) = b(1:n,1:nb);
  rhs(n+2:2*n+1,1:nb) = 0.0;

  il = n;
  ndiv = 1;
  ipntp = 0;

  while ( 1 < il )

    ipnt = ipntp;
    ipntp = ipntp + il;
    il = floor ( il / 2 );
    ndiv = ndiv * 2;
    ihaf = ipntp;

    for iful = ipnt + 2 : 2 : ipntp
      ihaf = ihaf + 1;
      rhs(ihaf+1,1:nb) = rhs(iful+1,1:nb) ...
        - a_cr(3,iful) * rhs(iful,1:nb) ...
        - a_cr(1,iful+1) * rhs(iful+2,1:nb);
    end

  end

  rhs(ihaf+1,1:nb) = rhs(ihaf+1,1:nb) * a_cr(2,ihaf+1);
  ipnt = ipntp;

  while ( 0 < ipnt )

    ipntp = ipnt;
    ndiv = floor ( ndiv / 2 );
    il = floor ( n / ndiv );
    ipnt = ipnt - il;
    ihaf = ipntp;

    for ifulm = ipnt + 1 : 2 : ipntp
      iful = ifulm + 1;
      ihaf = ihaf + 1;
      rhs(iful+1,1:nb) = rhs(ihaf+1,1:nb);
      rhs(ifulm+1,1:nb) = a_cr(2,ifulm+1) * ( rhs(ifulm+1,1:nb) ...
        - a_cr(3,ifulm) * rhs(ifulm,1:nb) ...
        - a_cr(1,ifulm+1) * rhs(iful+1,1:nb) );
    end

  end

  x(1:n,1:nb) = rhs(2:n+1,1:nb);

  return
end
