[n,c,d,e,b,info]=cgtsl(n,c,d,e,b,info);

function [n,c,d,e,b,info]=cgtsl(n,c,d,e,b,info);
%***BEGIN PROLOGUE  CGTSL
%***PURPOSE  Solve a tridiagonal linear system.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2C2A
%***TYPE      COMPLEX (SGTSL-S, DGTSL-D, CGTSL-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, SOLVE, TRIDIAGONAL
%***AUTHOR  Dongarra, J., (ANL)
%***DESCRIPTION
%
%     CGTSL given a general tridiagonal matrix and a right hand
%     side will find the solution.
%
%     On Entry
%
%        N       INTEGER
%                is the order of the tridiagonal matrix.
%
%        C       COMPLEX(N)
%                is the subdiagonal of the tridiagonal matrix.
%                C(2) through C(N) should contain the subdiagonal.
%                On output C is destroyed.
%
%        D       COMPLEX(N)
%                is the diagonal of the tridiagonal matrix.
%                On output D is destroyed.
%
%        E       COMPLEX(N)
%                is the superdiagonal of the tridiagonal matrix.
%                E(1) through E(N-1) should contain the superdiagonal.
%                On output E is destroyed.
%
%        B       COMPLEX(N)
%                is the right hand side vector.
%
%     On Return
%
%        B       is the solution vector.
%
%        INFO    INTEGER
%                = 0 normal value.
%                = K if the K-th element of the diagonal becomes
%                    exactly zero.  The subroutine returns when
%                    this is detected.
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   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  CGTSL
persistent k kb kp1 nm1 nm2 t zdum ; 

c_shape=size(c);c=reshape(c,1,[]);
d_shape=size(d);d=reshape(d,1,[]);
e_shape=size(e);e=reshape(e,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
%
if isempty(k), k=0; end;
if isempty(kb), kb=0; end;
if isempty(kp1), kp1=0; end;
if isempty(nm1), nm1=0; end;
if isempty(nm2), nm2=0; end;
if isempty(t), t=0; end;
if isempty(zdum), zdum=0; end;
% cabs1= @(zdum)  abs(real(zdum)) + abs(aimag(zdum));real :: cabs1;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
%***FIRST EXECUTABLE STATEMENT  CGTSL
info = 0;
c(1) = d(1);
nm1 = fix(n - 1);
if( nm1>=1 )
d(1) = e(1);
e(1) =complex(0.0e0,0.0e0);
e(n) =complex(0.0e0,0.0e0);
%
for k = 1 : nm1;
kp1 = fix(k + 1);
%
%              FIND THE LARGEST OF THE TWO ROWS
%
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(c(kp1))>=cabs1(c(k)) )
%
%                 INTERCHANGE ROW
%
t = c(kp1);
c(kp1) = c(k);
c(k) = t;
t = d(kp1);
d(kp1) = d(k);
d(k) = t;
t = e(kp1);
e(kp1) = e(k);
e(k) = t;
t = b(kp1);
b(kp1) = b(k);
b(k) = t;
end;
%
%              ZERO ELEMENTS
%
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(c(k))~=0.0e0 )
t = -c(kp1)./c(k);
c(kp1) = d(kp1) + t.*d(k);
d(kp1) = e(kp1) + t.*e(k);
e(kp1) =complex(0.0e0,0.0e0);
b(kp1) = b(kp1) + t.*b(k);
else;
info = fix(k);
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
return;
end;
end; k = fix(nm1+1);
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(c(n))~=0.0e0 )
%
%           BACK SOLVE
%
nm2 = fix(n - 2);
b(n) = b(n)./c(n);
if( n~=1 )
b(nm1) =(b(nm1)-d(nm1).*b(n))./c(nm1);
if( nm2>=1 )
for kb = 1 : nm2;
k = fix(nm2 - kb + 1);
b(k) =(b(k)-d(k).*b(k+1)-e(k).*b(k+2))./c(k);
end; kb = fix(nm2+1);
end;
end;
else;
info = fix(n);
end;
%
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
d_shape=zeros(d_shape);d_shape(:)=d(1:numel(d_shape));d=d_shape;
e_shape=zeros(e_shape);e_shape(:)=e(1:numel(e_shape));e=e_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end
%DECK CHBMV