[abd,lda,n,m,info]=spbfa(abd,lda,n,m,info);

function [abd,lda,n,m,info]=spbfa(abd,lda,n,m,info);
%***BEGIN PROLOGUE  SPBFA
%***PURPOSE  Factor a real symmetric positive definite matrix stored in
%            band form.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B2
%***TYPE      SINGLE PRECISION (SPBFA-S, DPBFA-D, CPBFA-C)
%***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION,
%             POSITIVE DEFINITE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SPBFA factors a real symmetric positive definite matrix
%     stored in band form.
%
%     SPBFA is usually called by SPBCO, but it can be called
%     directly with a saving in time if  RCOND  is not needed.
%
%     On Entry
%
%        ABD     REAL(LDA, N)
%                the matrix to be factored.  The columns of the upper
%                triangle are stored in the columns of ABD and the
%                diagonals of the upper triangle are stored in the
%                rows of ABD .  See the comments below for details.
%
%        LDA     INTEGER
%                the leading dimension of the array  ABD .
%                LDA must be .GE. M + 1 .
%
%        N       INTEGER
%                the order of the matrix  A .
%
%        M       INTEGER
%                the number of diagonals above the main diagonal.
%                0 .LE. M .LT. N .
%
%     On Return
%
%        ABD     an upper triangular matrix  R , stored in band
%                form, so that  A = TRANS(R)*R .
%
%        INFO    INTEGER
%                = 0  for normal return.
%                = K  if the leading minor of order  K  is not
%                     positive definite.
%
%     Band Storage
%
%           If  A  is a symmetric positive definite band matrix,
%           the following program segment will set up the input.
%
%                   M = (band width above diagonal)
%                   DO 20 J = 1, N
%                      I1 = MAX(1, J-M)
%                      DO 10 I = I1, J
%                         K = I-J+M+1
%                         ABD(K,J) = A(I,J)
%                10    CONTINUE
%                20 CONTINUE
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SDOT
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  SPBFA
persistent ik j jk k mu s t ; 

abd_shape=size(abd);abd=reshape([abd(:).',zeros(1,ceil(numel(abd)./prod([lda])).*prod([lda])-numel(abd))],lda,[]);
%
if isempty(t), t=0; end;
if isempty(s), s=0; end;
if isempty(ik), ik=0; end;
if isempty(j), j=0; end;
if isempty(jk), jk=0; end;
if isempty(k), k=0; end;
if isempty(mu), mu=0; end;
%***FIRST EXECUTABLE STATEMENT  SPBFA
for j = 1 : n;
info = fix(j);
s = 0.0e0;
ik = fix(m + 1);
jk = fix(max(j-m,1));
mu = fix(max(m+2-j,1));
if( m>=mu )
for k = mu : m;
t = abd(k,j) - sdot(k-mu,abd(sub2ind(size(abd),ik,jk):end),1,abd(sub2ind(size(abd),mu,j):end),1);
t = t./abd(m+1,jk);
abd(k,j) = t;
s = s + t.*t;
ik = fix(ik - 1);
jk = fix(jk + 1);
end; k = fix(m+1);
end;
s = abd(m+1,j) - s;
if( s<=0.0e0 )
abd_shape=zeros(abd_shape);abd_shape(:)=abd(1:numel(abd_shape));abd=abd_shape;
return;
end;
abd(m+1,j) = sqrt(s);
end; j = fix(n+1);
info = 0;
abd_shape=zeros(abd_shape);abd_shape(:)=abd(1:numel(abd_shape));abd=abd_shape;
end
%DECK SPBSL