function [ a_lu, info ] = r8cb_np_fa ( n, ml, mu, a )
%*****************************************************************************80
%
%% R8CB_NP_FA factors a R8CB matrix by Gaussian elimination.
%
% Discussion:
%
% The R8CB storage format is appropriate for a compact banded matrix.
% It is assumed that the matrix has lower and upper bandwidths ML and MU,
% respectively. The matrix is stored in a way similar to that used
% by LINPACK and LAPACK for a general banded matrix, except that in
% this mode, no extra rows are set aside for possible fillin during pivoting.
% Thus, this storage mode is suitable if you do not intend to factor
% the matrix, or if you can guarantee that the matrix can be factored
% without pivoting.
%
% R8CB_NP_FA is a version of the LINPACK routine R8GBFA, modifed to use
% no pivoting, and to be applied to the R8CB compressed band matrix storage
% format. It will fail if the matrix is singular, or if any zero
% pivot is encountered.
%
% If R8CB_NP_FA successfully factors the matrix, R8CB_NP_SL may be called
% to solve linear systems involving the matrix.
%
% The matrix is stored in a compact version of LINPACK general
% band storage, which does not include the fill-in entires.
% The following program segment will store the entries of a banded
% matrix in the compact format used by this routine:
%
% m = mu+1
% do j = 1, n
% i1 = max ( 1, j-mu )
% i2 = min ( n, j+ml )
% do i = i1, i2
% k = i-j+m
% a(k,j) = afull(i,j)
% end do
% end do
%
% Licensing:
%
% This code is distributed under the MIT license.
%
% Modified:
%
% 23 February 2004
%
% Author:
%
% John Burkardt
%
% Parameters:
%
% Input, integer N, the order of the matrix.
% N must be positive.
%
% Input, integer ML, MU, the lower and upper bandwidths.
% ML and MU must be nonnegative, and no greater than N-1.
%
% Input, real A(ML+MU+1,N), the compact band matrix.
%
% Output, real A_LU(ML+MU+1,N), the LU factors of the matrix.
%
% Output, integer INFO, singularity flag.
% 0, no singularity detected.
% nonzero, the factorization failed on the INFO-th step.
%
a_lu(1:ml+mu+1,1:n) = a(1:ml+mu+1,1:n);
%
% The value of M is MU + 1 rather than ML + MU + 1.
%
m = mu + 1;
info = 0;
ju = 0;
for k = 1 : n-1
%
% If our pivot entry A(MU+1,K) is zero, then we must give up.
%
if ( a_lu(m,k) == 0.0 )
info = k;
fprintf ( 1, '\n' );
fprintf ( 1, 'R8CB_NP_FA - Fatal error!\n' );
fprintf ( 1, ' Zero pivot on step %d\n', info );
error ( 'R8CB_NP_FA - Fatal error!' );
return;
end
%
% LM counts the number of nonzero elements that lie below the current
% diagonal entry, A(K,K).
%
% Multiply the LM entries below the diagonal by -1/A(K,K), turning
% them into the appropriate "multiplier" terms in the L matrix.
%
lm = min ( ml, n-k );
a_lu(m+1:m+lm,k) = -a_lu(m+1:m+lm,k) / a_lu(m,k);
%
% MM points to the row in which the next entry of the K-th row is, A(K,J).
% We then add L(I,K)*A(K,J) to A(I,J) for rows I = K+1 to K+LM.
%
ju = max ( ju, mu + k );
ju = min ( ju, n );
mm = m;
for j = k+1 : ju
mm = mm - 1;
a_lu(mm+1:mm+lm,j) = a_lu(mm+1:mm+lm,j) + a_lu(mm,j) * a_lu(m+1:m+lm,k);
end
end
if ( a_lu(m,n) == 0.0 )
info = n;
fprintf ( 1, '\n' );
fprintf ( 1, 'R8CB_NP_FA - Fatal error!\n' );
fprintf ( 1, ' Zero pivot on step %d\n', info );
error ( 'R8CB_NP_FA - Fatal error!' );
return;
end
return
end