[abd,lda,n,m,b]=spbsl(abd,lda,n,m,b);

function [abd,lda,n,m,b]=spbsl(abd,lda,n,m,b);
%***BEGIN PROLOGUE  SPBSL
%***PURPOSE  Solve a real symmetric positive definite band system
%            using the factors computed by SPBCO or SPBFA.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B2
%***TYPE      SINGLE PRECISION (SPBSL-S, DPBSL-D, CPBSL-C)
%***KEYWORDS  BANDED, LINEAR ALGEBRA, LINPACK, MATRIX,
%             POSITIVE DEFINITE, SOLVE
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     SPBSL solves the real symmetric positive definite band
%     system  A*X = B
%     using the factors computed by SPBCO or SPBFA.
%
%     On Entry
%
%        ABD     REAL(LDA, N)
%                the output from SPBCO or SPBFA.
%
%        LDA     INTEGER
%                the leading dimension of the array  ABD .
%
%        N       INTEGER
%                the order of the matrix  A .
%
%        M       INTEGER
%                the number of diagonals above the main diagonal.
%
%        B       REAL(N)
%                the right hand side vector.
%
%     On Return
%
%        B       the solution vector  X .
%
%     Error Condition
%
%        A division by zero will occur if the input factor contains
%        a zero on the diagonal.  Technically, this indicates
%        singularity, but it is usually caused by improper subroutine
%        arguments.  It will not occur if the subroutines are called
%        correctly and  INFO .EQ. 0 .
%
%     To compute  INVERSE(A) * C  where  C  is a matrix
%     with  P  columns
%           CALL SPBCO(ABD,LDA,N,RCOND,Z,INFO)
%           IF (RCOND is too small .OR. INFO .NE. 0) GO TO ...
%           DO 10 J = 1, P
%              CALL SPBSL(ABD,LDA,N,C(1,J))
%        10 CONTINUE
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  SAXPY, 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  SPBSL
persistent k kb la lb lm t ; 

abd_shape=size(abd);abd=reshape([abd(:).',zeros(1,ceil(numel(abd)./prod([lda])).*prod([lda])-numel(abd))],lda,[]);
b_shape=size(b);b=reshape(b,1,[]);
%
if isempty(t), t=0; end;
if isempty(k), k=0; end;
if isempty(kb), kb=0; end;
if isempty(la), la=0; end;
if isempty(lb), lb=0; end;
if isempty(lm), lm=0; end;
%
%     SOLVE TRANS(R)*Y = B
%
%***FIRST EXECUTABLE STATEMENT  SPBSL
for k = 1 : n;
lm = fix(min(k-1,m));
la = fix(m + 1 - lm);
lb = fix(k - lm);
[t ,lm,abd(sub2ind(size(abd),la,k):end),dumvar4,b(sub2ind(size(b),max(lb,1)):end)]=sdot(lm,abd(sub2ind(size(abd),la,k):end),1,b(sub2ind(size(b),max(lb,1)):end),1);
b(k) =(b(k)-t)./abd(m+1,k);
end; k = fix(n+1);
%
%     SOLVE R*X = Y
%
for kb = 1 : n;
k = fix(n + 1 - kb);
lm = fix(min(k-1,m));
la = fix(m + 1 - lm);
lb = fix(k - lm);
b(k) = b(k)./abd(m+1,k);
t = -b(k);
[lm,t,abd(sub2ind(size(abd),la,k):end),dumvar4,b(sub2ind(size(b),max(lb,1)):end)]=saxpy(lm,t,abd(sub2ind(size(abd),la,k):end),1,b(sub2ind(size(b),max(lb,1)):end),1);
end; kb = fix(n+1);
abd_shape=zeros(abd_shape);abd_shape(:)=abd(1:numel(abd_shape));abd=abd_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK SPELI4