[abe,lda,n,ml,mu,ipvt,b,job]=cnbsl(abe,lda,n,ml,mu,ipvt,b,job);

function [abe,lda,n,ml,mu,ipvt,b,job]=cnbsl(abe,lda,n,ml,mu,ipvt,b,job);
%***BEGIN PROLOGUE  CNBSL
%***PURPOSE  Solve a complex band system using the factors computed by
%            CNBCO or CNBFA.
%***LIBRARY   SLATEC
%***CATEGORY  D2C2
%***TYPE      COMPLEX (SNBSL-S, DNBSL-D, CNBSL-C)
%***KEYWORDS  BANDED, LINEAR EQUATIONS, NONSYMMETRIC, SOLVE
%***AUTHOR  Voorhees, E. A., (LANL)
%***DESCRIPTION
%
%     CNBSL solves the complex band system
%     A * X = B  or  CTRANS(A) * X = B
%     using the factors computed by CNBCO or CNBFA.
%
%     On Entry
%
%        ABE     COMPLEX(LDA, NC)
%                the output from CNBCO or CNBFA.
%                NC must be .GE. 2*ML+MU+1 .
%
%        LDA     INTEGER
%                the leading dimension of the array  ABE .
%
%        N       INTEGER
%                the order of the original matrix.
%
%        ML      INTEGER
%                number of diagonals below the main diagonal.
%
%        MU      INTEGER
%                number of diagonals above the main diagonal.
%
%        IPVT    INTEGER(N)
%                the pivot vector from CNBCO or CNBFA.
%
%        B       COMPLEX(N)
%                the right hand side vector.
%
%        JOB     INTEGER
%                = 0         to solve  A*X = B .
%                = nonzero   to solve  CTRANS(A)*X = B , where
%                            CTRANS(A)  is the conjugate transpose.
%
%     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 often caused by improper arguments or improper
%        setting of LDA.  It will not occur if the subroutines are
%        called correctly and if CNBCO has set RCOND .GT. 0.0
%        or CNBFA has set INFO .EQ. 0 .
%
%     To compute  INVERSE(A) * C  where  C  is a matrix
%     with  P  columns
%           CALL CNBCO(ABE,LDA,N,ML,MU,IPVT,RCOND,Z)
%           IF (RCOND is too small) GO TO ...
%           DO 10 J = 1, P
%             CALL CNBSL(ABE,LDA,N,ML,MU,IPVT,C(1,J),0)
%        10 CONTINUE
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  CAXPY, CDOTC
%***REVISION HISTORY  (YYMMDD)
%   800730  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)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  CNBSL
persistent k kb l lb ldb lm m mlm nm1 t ; 

ipvt_shape=size(ipvt);ipvt=reshape(ipvt,1,[]);
abe_shape=size(abe);abe=reshape([abe(:).',zeros(1,ceil(numel(abe)./prod([lda])).*prod([lda])-numel(abe))],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(l), l=0; end;
if isempty(lb), lb=0; end;
if isempty(ldb), ldb=0; end;
if isempty(lm), lm=0; end;
if isempty(m), m=0; end;
if isempty(mlm), mlm=0; end;
if isempty(nm1), nm1=0; end;
%***FIRST EXECUTABLE STATEMENT  CNBSL
m = fix(mu + ml + 1);
nm1 = fix(n - 1);
ldb = fix(1 - lda);
if( job~=0 )
%
%       JOB = NONZERO, SOLVE CTRANS(A) * X = B
%       FIRST SOLVE  CTRANS(U)*Y = B
%
for k = 1 : n;
lm = fix(min(k,m) - 1);
lb = fix(k - lm);
[t ,lm,abe(sub2ind(size(abe),k-1,ml+2):end),ldb,b(sub2ind(size(b),max(lb,1)):end)]=cdotc(lm,abe(sub2ind(size(abe),k-1,ml+2):end),ldb,b(sub2ind(size(b),max(lb,1)):end),1);
b(k) =(b(k)-t)./conj(abe(k,ml+1));
end; k = fix(n+1);
%
%       NOW SOLVE CTRANS(L)*X = Y
%
if( ml~=0 )
if( nm1>=1 )
for kb = 1 : nm1;
k = fix(n - kb);
lm = fix(min(ml,n-k));
mlm = fix(ml -(lm-1));
b(k) = b(k) + cdotc(lm,abe(sub2ind(size(abe),k+lm,mlm):end),ldb,b(sub2ind(size(b),max(k+1,1)):end),1);
l = fix(ipvt(k));
if( l~=k )
t = b(l);
b(l) = b(k);
b(k) = t;
end;
end; kb = fix(nm1+1);
end;
end;
else;
%
%       JOB = 0 , SOLVE  A * X = B
%       FIRST SOLVE L*Y = B
%
if( ml~=0 )
if( nm1>=1 )
for k = 1 : nm1;
lm = fix(min(ml,n-k));
l = fix(ipvt(k));
t = b(l);
if( l~=k )
b(l) = b(k);
b(k) = t;
end;
mlm = fix(ml -(lm-1));
[lm,t,abe(sub2ind(size(abe),k+lm,mlm):end),ldb,b(sub2ind(size(b),max(k+1,1)):end)]=caxpy(lm,t,abe(sub2ind(size(abe),k+lm,mlm):end),ldb,b(sub2ind(size(b),max(k+1,1)):end),1);
end; k = fix(nm1+1);
end;
end;
%
%       NOW SOLVE  U*X = Y
%
for kb = 1 : n;
k = fix(n + 1 - kb);
b(k) = b(k)./abe(k,ml+1);
lm = fix(min(k,m) - 1);
lb = fix(k - lm);
t = -b(k);
%barrowes, modified to have max(,1) because 0-indexed when k=1
[lm,t,abe(sub2ind(size(abe),max(k-1,1),ml+2):end),ldb,b(sub2ind(size(b),max(lb,1)):end)]=caxpy(lm,t,abe(sub2ind(size(abe),max(k-1,1),ml+2):end),ldb,b(sub2ind(size(b),max(lb,1)):end),1);
end; kb = fix(n+1);
end;
ipvt_shape=zeros(ipvt_shape);ipvt_shape(:)=ipvt(1:numel(ipvt_shape));ipvt=ipvt_shape;
abe_shape=zeros(abe_shape);abe_shape(:)=abe(1:numel(abe_shape));abe=abe_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK COMBAK