function x = r8bto_sl ( m, l, a, b )

%*****************************************************************************80
%
%% R8BTO_SL solves an R8BTO system.
%
%  Discussion:
%
%    The R8BTO storage format is for a block Toeplitz matrix. The matrix
%    can be regarded as an L by L array of blocks, each of size M by M.
%    The full matrix has order N = M * L.  The L by L matrix is Toeplitz,
%    that is, along its diagonal, the blocks repeat.
%
%    Storage for the matrix consists of the L blocks of the first row,
%    followed by the L-1 blocks of the first column (skipping the first row).
%    These items are stored in the natural way in an (M,M,2*L-1) array.
%
%  Example:
%
%    M = 2, L = 3
%
%    1 2 | 3 4 | 5 6
%    5 5 | 6 6 | 7 7
%    ----+-----+-----
%    7 8 | 1 2 | 3 4
%    8 8 | 5 5 | 6 6
%    ----+-----+-----
%    9 0 | 7 8 | 1 2
%    9 9 | 8 8 | 5 5
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    08 July 2016
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the order of the blocks of the matrix A.
%
%    Input, integer L, the number of blocks in a row or column.
%
%    Input, real A(M,M,2*L-1), the R8BTO matrix.
%
%    Input, real B(M*L), the right hand side vector.
%
%    Output, real X(M*L), the solution vector.  X and B
%    may share storage.
%
  bmat = reshape ( b, m, l );
  c1 = zeros ( m, m, l - 1 );
  c2 = zeros ( m, m, l - 1 );
  pivot = zeros ( m );
  r1 = zeros ( m, m );
  r2 = zeros ( m, m );
  r3 = zeros ( m, m );
  r5 = zeros ( m, m );
  r6 = zeros ( m, m );
  xmat = zeros ( m, l );
%
%  Solve the system with the principal minor of order M.
%
  i3 = 1;
  for j = 1 : m
    for i = 1 : m
      c1(i,j,1) = a(i,j,1);
      r1(i,j) = a(i,j,1);
      i3 = i3 + 1;
    end
  end

  r3(1:m,1:m) = r1(1:m,1:m);
  xmat(1:m,1) = b(1:m,1);

  [ r3, pivot, info ] = r8ge_fa ( m, r3 );

  job = 0;
  xmat(1:m,1) = r8ge_sl ( m, r3, pivot, xmat(1:m,1), job );

  if ( l == 1 )
    return
  end
%
%  Recurrent process for solving the system
%  with the block Toeplitz matrix for N = 2 through L.
%
  for n = 2 : l
%
%  Compute multiples of the first and last block columns of
%  the inverse of the principal minor of order M*N.
%
    i3 = 1;
    for j = 1 : m
      for i = 1 : m
        r5(i,j) = a(i,j,l+n-1);
        r6(i,j) = a(i,j,n);
        i3 = i3 + 1;
      end
    end

    if ( 2 < n )

      c1(1:m,1:m,n-1) = r2(1:m,1:m);

      for i1 = 1 : n - 2
        i2 = n - i1;
        for j = 1 : m
          i3 = 1;
%
%  My apologies, but we have an I1+L, followed in the next line by I1+1.
%
          for i = 1 : m
            r5(1:m,j) = r5(1:m,j) + c1(i,j,i2) * a(1:m,i,i1+l);
            r6(1:m,j) = r6(1:m,j) + c2(i,j,i1) * a(1:m,i,i1+1);
            i3 = i3 + m;
          end
        end
      end

    end

    for j = 1 : m
      r2(1:m,j) = - r5(1:m,j);
      job = 0;
      r2(1:m,j) = r8ge_sl ( m, r3, pivot, r2(1:m,j), job );
    end

    r3(1:m,1:m) = r6(1:m,1:m);
    r6(1:m,1:m) = - c1(1:m,1:m,1);

    for j = 1 : m
      for i = 1 : m
        c1(1:m,j,1) = c1(1:m,j,1) + r2(i,j) * r3(1:m,i);
      end
    end

    [ r6, pivot, info ] = r8ge_fa ( m, r6 );

    for j = 1 : m
      r3(1:m,j) = r8ge_sl ( m, r6, pivot, r3(1:m,j), job );
      for i = 1 : m
        r1(1:m,j) = r1(1:m,j) + r3(i,j) * r5(1:m,i);
      end
    end

    if ( 2 < n )

      r6(1:m,1:m) = c2(1:m,1:m,1);

      for i1 = 2 : n - 1

        if ( i1 ~= n - 1 )
          r5(1:m,1:m) = c2(1:m,1:m,i1);
        end

        for j = 1 : m
          c2(1:m,j,i1) = r6(1:m,j);
          for i = 1 : m
            c2(1:m,j,i1) = c2(1:m,j,i1) + r3(i,j) * c1(1:m,i,i1);
          end
        end

        for j = 1 : m
          for i = 1 : m
            c1(1:m,j,i1) = c1(1:m,j,i1) + r2(i,j) * r6(1:m,i);
          end
        end

        r6(1:m,1:m) = r5(1:m,1:m);

      end

    end

    c2(1:m,1:m,1) = r3(1:m,1:m);
%
%  Compute the solution of the system with the principal minor of order M*N.
%
    r3(1:m,1:m) = r1(1:m,1:m);
    xmat(1:m,n) = bmat(1:m,n);

    for i1 = 1 : n - 1
      i2 = n - i1;
      i3 = 1;
      for i = 1 : m
        for jack = 1 : m
          xmat(jack,n) = xmat(jack,n) - xmat(i,i2) * a(jack,i,i1+l);
        end
        i3 = i3 + m;
      end
    end

    [ r3, pivot, info ] = r8ge_fa ( m, r3 );

    xmat(1:m,n) = r8ge_sl ( m, r3, pivot, xmat(1:m,n), job );

    for i1 = 1 : n - 1
      for i = 1 : m
        xmat(1:m,i1) = xmat(1:m,i1) + xmat(i,n) * c2(1:m,i,i1);
      end
    end

  end

  x = reshape ( xmat, m * l, 1 );

  return
end
