Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[abd,lda,n,ml,mu,ipvt,info]=sgbfa(abd,lda,n,ml,mu,ipvt,info);
```function [abd,lda,n,ml,mu,ipvt,info]=sgbfa(abd,lda,n,ml,mu,ipvt,info);
%***BEGIN PROLOGUE  SGBFA
%***PURPOSE  Factor a band matrix using Gaussian elimination.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2A2
%***TYPE      SINGLE PRECISION (SGBFA-S, DGBFA-D, CGBFA-C)
%***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SGBFA factors a real band matrix by elimination.
%
%     SGBFA is usually called by SBGCO, but it can be called
%     directly with a saving in time if  RCOND  is not needed.
%
%     On Entry
%
%        ABD     REAL(LDA, N)
%                contains the matrix in band storage.  The columns
%                of the matrix are stored in the columns of  ABD  and
%                the diagonals of the matrix are stored in rows
%                ML+1 through 2*ML+MU+1 of  ABD .
%                See the comments below for details.
%
%        LDA     INTEGER
%                the leading dimension of the array  ABD .
%                LDA must be .GE. 2*ML + MU + 1 .
%
%        N       INTEGER
%                the order of the original matrix.
%
%        ML      INTEGER
%                number of diagonals below the main diagonal.
%                0 .LE. ML .LT. N .
%
%        MU      INTEGER
%                number of diagonals above the main diagonal.
%                0 .LE. MU .LT. N .
%                More efficient if  ML .LE. MU .
%     On Return
%
%        ABD     an upper triangular matrix in band storage and
%                the multipliers which were used to obtain it.
%                The factorization can be written  A = L*U , where
%                L  is a product of permutation and unit lower
%                triangular matrices and  U  is upper triangular.
%
%        IPVT    INTEGER(N)
%                an integer vector of pivot indices.
%
%        INFO    INTEGER
%                = 0  normal value.
%                = K  if  U(K,K) .EQ. 0.0 .  This is not an error
%                     condition for this subroutine, but it does
%                     indicate that SGBSL will divide by zero if
%                     called.  Use  RCOND  in SBGCO for a reliable
%                     indication of singularity.
%
%     Band Storage
%
%           If  A  is a band matrix, the following program segment
%           will set up the input.
%
%                   ML = (band width below the diagonal)
%                   MU = (band width above the diagonal)
%                   M = ML + MU + 1
%                   DO 20 J = 1, N
%                      I1 = MAX(1, J-MU)
%                      I2 = MIN(N, J+ML)
%                      DO 10 I = I1, I2
%                         K = I - J + M
%                         ABD(K,J) = A(I,J)
%                10    CONTINUE
%                20 CONTINUE
%
%           This uses rows  ML+1  through  2*ML+MU+1  of  ABD .
%           In addition, the first  ML  rows in  ABD  are used for
%           elements generated during the triangularization.
%           The total number of rows needed in  ABD  is  2*ML+MU+1 .
%           The  ML+MU by ML+MU  upper left triangle and the
%           ML by ML  lower right triangle are not referenced.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  ISAMAX, SAXPY, SSCAL
%***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  SGBFA
persistent i i0 j j0 j1 ju jz k kp1 l lm m mm nm1 t ;

ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
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(i), i=0; end;
if isempty(i0), i0=0; end;
if isempty(j), j=0; end;
if isempty(ju), ju=0; end;
if isempty(jz), jz=0; end;
if isempty(j0), j0=0; end;
if isempty(j1), j1=0; end;
if isempty(k), k=0; end;
if isempty(kp1), kp1=0; end;
if isempty(l), l=0; end;
if isempty(lm), lm=0; end;
if isempty(m), m=0; end;
if isempty(mm), mm=0; end;
if isempty(nm1), nm1=0; end;
%
%***FIRST EXECUTABLE STATEMENT  SGBFA
m = fix(ml + mu + 1);
info = 0;
%
%     ZERO INITIAL FILL-IN COLUMNS
%
j0 = fix(mu + 2);
j1 = fix(min(n,m) - 1);
if( j1>=j0 )
for jz = j0 : j1;
i0 = fix(m + 1 - jz);
for i = i0 : ml;
abd(i,jz) = 0.0e0;
end; i = fix(ml+1);
end; jz = fix(j1+1);
end;
jz = fix(j1);
ju = 0;
%
%     GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
%
nm1 = fix(n - 1);
if( nm1>=1 )
for k = 1 : nm1;
kp1 = fix(k + 1);
%
%        ZERO NEXT FILL-IN COLUMN
%
jz = fix(jz + 1);
if( jz<=n )
if( ml>=1 )
for i = 1 : ml;
abd(i,jz) = 0.0e0;
end; i = fix(ml+1);
end;
end;
%
%        FIND L = PIVOT INDEX
%
lm = fix(min(ml,n-k));
l = fix(isamax(lm+1,abd(sub2ind(size(abd),m,k):end),1) + m - 1);
ipvt(k) = fix(l + k - m);
%
%        ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
%
if( abd(l,k)==0.0e0 )
info = fix(k);
else;
%
%           INTERCHANGE IF NECESSARY
%
if( l~=m )
t = abd(l,k);
abd(l,k) = abd(m,k);
abd(m,k) = t;
end;
%
%           COMPUTE MULTIPLIERS
%
t = -1.0e0./abd(m,k);
[lm,t,abd(sub2ind(size(abd),m+1,k):end)]=sscal(lm,t,abd(sub2ind(size(abd),m+1,k):end),1);
%
%           ROW ELIMINATION WITH COLUMN INDEXING
%
ju = fix(min(max(ju,mu+ipvt(k)),n));
mm = fix(m);
if( ju>=kp1 )
for j = kp1 : ju;
l = fix(l - 1);
mm = fix(mm - 1);
t = abd(l,j);
if( l~=mm )
abd(l,j) = abd(mm,j);
abd(mm,j) = t;
end;
[lm,t,dumvar3,dumvar4,dumvar5]=saxpy(lm,t,abd(sub2ind(size(abd),m+1,k):end),1,abd(sub2ind(size(abd),mm+1,j):end),1);      abd(sub2ind(size(abd),m+1,k):end)=dumvar3; abd(sub2ind(size(abd),mm+1,j):end)=dumvar5;
end; j = fix(ju+1);
end;
end;
end; k = fix(nm1+1);
end;
ipvt(n) = fix(n);
if( abd(m,n)==0.0e0 )
info = fix(n);
end;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
abd_shape=zeros(abd_shape);abd_shape(:)=abd(1:numel(abd_shape));abd=abd_shape;
end
%DECK SGBMV
```