[ap,n,kpvt,b]=chpsl(ap,n,kpvt,b);

function [ap,n,kpvt,b]=chpsl(ap,n,kpvt,b);
%***BEGIN PROLOGUE  CHPSL
%***PURPOSE  Solve a complex Hermitian system using factors obtained
%            from CHPFA.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2D1A
%***TYPE      COMPLEX (SSPSL-S, DSPSL-D, CHPSL-C, CSPSL-C)
%***KEYWORDS  HERMITIAN, LINEAR ALGEBRA, LINPACK, MATRIX, PACKED, SOLVE
%***AUTHOR  Bunch, J., (UCSD)
%***DESCRIPTION
%
%     CHISL solves the complex Hermitian system
%     A * X = B
%     using the factors computed by CHPFA.
%
%     On Entry
%
%        AP      COMPLEX(N*(N+1)/2)
%                the output from CHPFA.
%
%        N       INTEGER
%                the order of the matrix  A .
%
%        KVPT    INTEGER(N)
%                the pivot vector from CHPFA.
%
%        B       COMPLEX(N)
%                the right hand side vector.
%
%     On Return
%
%        B       the solution vector  X .
%
%     Error Condition
%
%        A division by zero may occur if  CHPCO  has set RCOND .EQ. 0.0
%        or  CHPFA  has set INFO .NE. 0  .
%
%     To compute  INVERSE(A) * C  where  C  is a matrix
%     with  P  columns
%           CALL CHPFA(AP,N,KVPT,INFO)
%           IF (INFO .NE. 0) GO TO ...
%           DO 10 J = 1, P
%              CALL CHPSL(AP,N,KVPT,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  CAXPY, CDOTC
%***REVISION HISTORY  (YYMMDD)
%   780814  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890831  Modified array declarations.  (WRB)
%   891107  Modified routine equivalence list.  (WRB)
%   891107  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  CHPSL
persistent ak akm1 bk bkm1 denom ik ikm1 ikp1 k kk km1k km1km1 kp temp ; 

kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
ap_shape=size(ap);ap=reshape(ap,1,[]);
b_shape=size(b);b=reshape(b,1,[]);
%
if isempty(ak), ak=0; end;
if isempty(akm1), akm1=0; end;
if isempty(bk), bk=0; end;
if isempty(bkm1), bkm1=0; end;
if isempty(denom), denom=0; end;
if isempty(temp), temp=0; end;
if isempty(ik), ik=0; end;
if isempty(ikm1), ikm1=0; end;
if isempty(ikp1), ikp1=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(km1k), km1k=0; end;
if isempty(km1km1), km1km1=0; end;
if isempty(kp), kp=0; end;
%
%     LOOP BACKWARD APPLYING THE TRANSFORMATIONS AND
%     D INVERSE TO B.
%
%***FIRST EXECUTABLE STATEMENT  CHPSL
k = fix(n);
ik =fix(fix((n.*(n-1))./2));
while( k~=0 );
kk = fix(ik + k);
if( kpvt(k)<0 )
%
%           2 X 2 PIVOT BLOCK.
%
ikm1 = fix(ik -(k-1));
if( k~=2 )
kp = fix(abs(kpvt(k)));
if( kp~=k-1 )
%
%                 INTERCHANGE.
%
temp = b(k-1);
b(k-1) = b(kp);
b(kp) = temp;
end;
%
%              APPLY THE TRANSFORMATION.
%
[dumvar1,dumvar2,ap(sub2ind(size(ap),max(ik+1,1)):end),dumvar4,dumvar5]=caxpy(k-2,b(k),ap(sub2ind(size(ap),max(ik+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);   dumvar2i=find((b(k))~=(dumvar2));dumvar5i=find((b(sub2ind(size(b),max(1,1)):end))~=(dumvar5));   b(k-1+dumvar2i)=dumvar2(dumvar2i); b(1-1+dumvar5i)=dumvar5(dumvar5i); 
[dumvar1,dumvar2,ap(sub2ind(size(ap),max(ikm1+1,1)):end),dumvar4,dumvar5]=caxpy(k-2,b(k-1),ap(sub2ind(size(ap),max(ikm1+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);   dumvar2i=find((b(k-1))~=(dumvar2));dumvar5i=find((b(sub2ind(size(b),max(1,1)):end))~=(dumvar5));   b(k-1-1+dumvar2i)=dumvar2(dumvar2i); b(1-1+dumvar5i)=dumvar5(dumvar5i); 
end;
%
%           APPLY D INVERSE.
%
km1k = fix(ik + k - 1);
kk = fix(ik + k);
ak = ap(kk)./conj(ap(km1k));
km1km1 = fix(ikm1 + k - 1);
akm1 = ap(km1km1)./ap(km1k);
bk = b(k)./conj(ap(km1k));
bkm1 = b(k-1)./ap(km1k);
denom = ak.*akm1 - 1.0e0;
b(k) =(akm1.*bk-bkm1)./denom;
b(k-1) =(ak.*bkm1-bk)./denom;
k = fix(k - 2);
ik = fix(ik -(k+1) - k);
else;
%
%           1 X 1 PIVOT BLOCK.
%
if( k~=1 )
kp = fix(kpvt(k));
if( kp~=k )
%
%                 INTERCHANGE.
%
temp = b(k);
b(k) = b(kp);
b(kp) = temp;
end;
%
%              APPLY THE TRANSFORMATION.
%
[dumvar1,dumvar2,ap(sub2ind(size(ap),max(ik+1,1)):end),dumvar4,dumvar5]=caxpy(k-1,b(k),ap(sub2ind(size(ap),max(ik+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);   dumvar2i=find((b(k))~=(dumvar2));dumvar5i=find((b(sub2ind(size(b),max(1,1)):end))~=(dumvar5));   b(k-1+dumvar2i)=dumvar2(dumvar2i); b(1-1+dumvar5i)=dumvar5(dumvar5i); 
end;
%
%           APPLY D INVERSE.
%
b(k) = b(k)./ap(kk);
k = fix(k - 1);
ik = fix(ik - k);
end;
end;
%
%     LOOP FORWARD APPLYING THE TRANSFORMATIONS.
%
k = 1;
ik = 0;
while( k<=n );
if( kpvt(k)<0 )
%
%           2 X 2 PIVOT BLOCK.
%
if( k~=1 )
%
%              APPLY THE TRANSFORMATION.
%
b(k) = b(k) + cdotc(k-1,ap(sub2ind(size(ap),max(ik+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);
ikp1 = fix(ik + k);
b(k+1) = b(k+1) + cdotc(k-1,ap(sub2ind(size(ap),max(ikp1+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);
kp = fix(abs(kpvt(k)));
if( kp~=k )
%
%                 INTERCHANGE.
%
temp = b(k);
b(k) = b(kp);
b(kp) = temp;
end;
end;
ik = fix(ik + k + k + 1);
k = fix(k + 2);
else;
%
%           1 X 1 PIVOT BLOCK.
%
if( k~=1 )
%
%              APPLY THE TRANSFORMATION.
%
b(k) = b(k) + cdotc(k-1,ap(sub2ind(size(ap),max(ik+1,1)):end),1,b(sub2ind(size(b),max(1,1)):end),1);
kp = fix(kpvt(k));
if( kp~=k )
%
%                 INTERCHANGE.
%
temp = b(k);
b(k) = b(kp);
b(kp) = temp;
end;
end;
ik = fix(ik + k);
k = fix(k + 1);
end;
end;
kpvt_shape=zeros(kpvt_shape);kpvt_shape(:)=kpvt(1:numel(kpvt_shape));kpvt=kpvt_shape;
ap_shape=zeros(ap_shape);ap_shape(:)=ap(1:numel(ap_shape));ap=ap_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
end %subroutine chpsl
%DECK CHU