[n,d,e,b]=sptsl(n,d,e,b);

function [n,d,e,b]=sptsl(n,d,e,b);
%***BEGIN PROLOGUE  SPTSL
%***PURPOSE  Solve a positive definite tridiagonal linear system.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2B2A
%***TYPE      SINGLE PRECISION (SPTSL-S, DPTSL-D, CPTSL-C)
%***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX, POSITIVE DEFINITE, SOLVE,
%             TRIDIAGONAL
%***AUTHOR  Dongarra, J., (ANL)
%***DESCRIPTION
%
%     SPTSL given a positive definite tridiagonal matrix and a right
%     hand side will find the solution.
%
%     On Entry
%
%        N        INTEGER
%                 is the order of the tridiagonal matrix.
%
%        D        REAL(N)
%                 is the diagonal of the tridiagonal matrix.
%                 On output, D is destroyed.
%
%        E        REAL(N)
%                 is the offdiagonal of the tridiagonal matrix.
%                 E(1) through E(N-1) should contain the
%                 offdiagonal.
%
%        B        REAL(N)
%                 is the right hand side vector.
%
%     On Return
%
%        B        contains the solution.
%
%***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
%   890505  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  SPTSL
persistent k kbm1 ke kf kp1 nm1 nm1d2 t1 t2 ; 

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(kbm1), kbm1=0; end;
if isempty(ke), ke=0; end;
if isempty(kf), kf=0; end;
if isempty(kp1), kp1=0; end;
if isempty(nm1), nm1=0; end;
if isempty(nm1d2), nm1d2=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
%
%     CHECK FOR 1 X 1 CASE
%
%***FIRST EXECUTABLE STATEMENT  SPTSL
if( n~=1 )
nm1 = fix(n - 1);
nm1d2 = fix(fix(nm1./2));
if( n~=2 )
kbm1 = fix(n - 1);
%
%           ZERO TOP HALF OF SUBDIAGONAL AND BOTTOM HALF OF
%           SUPERDIAGONAL
%
for k = 1 : nm1d2;
t1 = e(k)./d(k);
d(k+1) = d(k+1) - t1.*e(k);
b(k+1) = b(k+1) - t1.*b(k);
t2 = e(kbm1)./d(kbm1+1);
d(kbm1) = d(kbm1) - t2.*e(kbm1);
b(kbm1) = b(kbm1) - t2.*b(kbm1+1);
kbm1 = fix(kbm1 - 1);
end; k = fix(nm1d2+1);
end;
kp1 = fix(nm1d2 + 1);
%
%        CLEAN UP FOR POSSIBLE 2 X 2 BLOCK AT CENTER
%
if( rem(n,2)==0 )
t1 = e(kp1)./d(kp1);
d(kp1+1) = d(kp1+1) - t1.*e(kp1);
b(kp1+1) = b(kp1+1) - t1.*b(kp1);
kp1 = fix(kp1 + 1);
end;
%
%        BACK SOLVE STARTING AT THE CENTER, GOING TOWARDS THE TOP
%        AND BOTTOM
%
b(kp1) = b(kp1)./d(kp1);
if( n~=2 )
k = fix(kp1 - 1);
ke = fix(kp1 + nm1d2 - 1);
for kf = kp1 : ke;
b(k) =(b(k)-e(k).*b(k+1))./d(k);
b(kf+1) =(b(kf+1)-e(kf).*b(kf))./d(kf+1);
k = fix(k - 1);
end; kf = fix(ke+1);
end;
if( rem(n,2)==0 )
b(1) =(b(1)-e(1).*b(2))./d(1);
end;
else;
b(1) = b(1)./d(1);
end;
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 SQRDC