Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[a,lda,n,kpvt,rcond,z]=chico(a,lda,n,kpvt,rcond,z);
```function [a,lda,n,kpvt,rcond,z]=chico(a,lda,n,kpvt,rcond,z);
%***BEGIN PROLOGUE  CHICO
%***PURPOSE  Factor a complex Hermitian matrix by elimination with sym-
%            metric pivoting and estimate the condition of the matrix.
%***LIBRARY   SLATEC (LINPACK)
%***CATEGORY  D2D1A
%***TYPE      COMPLEX (SSICO-S, DSICO-D, CHICO-C, CSICO-C)
%***KEYWORDS  CONDITION NUMBER, HERMITIAN, LINEAR ALGEBRA, LINPACK,
%             MATRIX FACTORIZATION
%***AUTHOR  Moler, C. B., (U. of New Mexico)
%***DESCRIPTION
%
%     CHICO factors a complex Hermitian matrix by elimination with
%     symmetric pivoting and estimates the condition of the matrix.
%
%     If  RCOND  is not needed, CHIFA is slightly faster.
%     To solve  A*X = B , follow CHICO by CHISL.
%     To compute  INVERSE(A)*C , follow CHICO by CHISL.
%     To compute  INVERSE(A) , follow CHICO by CHIDI.
%     To compute  DETERMINANT(A) , follow CHICO by CHIDI.
%     To compute  INERTIA(A), follow CHICO by CHIDI.
%
%     On Entry
%
%        A       COMPLEX(LDA, N)
%                the Hermitian matrix to be factored.
%                Only the diagonal and upper triangle are used.
%
%        LDA     INTEGER
%                the leading dimension of the array  A .
%
%        N       INTEGER
%                the order of the matrix  A .
%
%     Output
%
%        A       a block diagonal matrix and the multipliers which
%                were used to obtain it.
%                The factorization can be written  A = U*D*CTRANS(U)
%                where  U  is a product of permutation and unit
%                upper triangular matrices , CTRANS(U) is the
%                conjugate transpose of  U , and  D  is block diagonal
%                with 1 by 1 and 2 by 2 blocks.
%
%        KVPT    INTEGER(N)
%                an integer vector of pivot indices.
%
%        RCOND   REAL
%                an estimate of the reciprocal condition of  A .
%                For the system  A*X = B , relative perturbations
%                in  A  and  B  of size  EPSILON  may cause
%                relative perturbations in  X  of size  EPSILON/RCOND .
%                If  RCOND  is so small that the logical expression
%                           1.0 + RCOND .EQ. 1.0
%                is truemlv, then  A  may be singular to working
%                precision.  In particular,  RCOND  is zero  if
%                exact singularity is detected or the estimate
%                underflows.
%
%        Z       COMPLEX(N)
%                a work vector whose contents are usually unimportant.
%                If  A  is close to a singular matrix, then  Z  is
%                an approximate null vector in the sense that
%                NORM(A*Z) = RCOND*NORM(A)*NORM(Z) .
%
%***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
%                 Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED  CAXPY, CDOTC, CHIFA, CSSCAL, SCASUM
%***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  CHICO
persistent ak akm1 anorm bk bkm1 denom ek i info j jm1 k kp kps ks s t ynorm zdum zdum2 ;

kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([lda])).*prod([lda])-numel(a))],lda,[]);
z_shape=size(z);z=reshape(z,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(ek), ek=0; end;
if isempty(t), t=0; end;
if isempty(anorm), anorm=0; end;
if isempty(s), s=0; end;
if isempty(ynorm), ynorm=0; end;
if isempty(i), i=0; end;
if isempty(info), info=0; end;
if isempty(j), j=0; end;
if isempty(jm1), jm1=0; end;
if isempty(k), k=0; end;
if isempty(kp), kp=0; end;
if isempty(kps), kps=0; end;
if isempty(ks), ks=0; end;
if isempty(zdum), zdum=0; end;
if isempty(zdum2), zdum2=0; end;
% csign1= @(zdum,zdum2)  cabs1(zdum)*(zdum2/cabs1(zdum2));complex :: csign1;
% cabs1= @(zdum)  abs(real(zdum)) + abs(aimag(zdum));real :: cabs1;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
csign1= @(zdum,zdum2)  cabs1(zdum).*(zdum2./cabs1(zdum2));
%
%     FIND NORM OF A USING ONLY UPPER HALF
%
%***FIRST EXECUTABLE STATEMENT  CHICO
for j = 1 : n;
z(j) = complex(scasum(j,a(sub2ind(size(a),1,j):end),1),0.0e0);
jm1 = fix(j - 1);
if( jm1>=1 )
for i = 1 : jm1;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
z(i) = complex(real(z(i))+cabs1(a(i,j)),0.0e0);
end; i = fix(jm1+1);
end;
end; j = fix(n+1);
anorm = 0.0e0;
for j = 1 : n;
anorm = max(anorm,real(z(j)));
end; j = fix(n+1);
%
%     FACTOR
%
[a,lda,n,kpvt,info]=chifa(a,lda,n,kpvt,info);
%
%     RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) .
%     ESTIMATE = NORM(Z)/NORM(Y) WHERE  A*Z = Y  AND  A*Y = E .
%     THE COMPONENTS OF  E  ARE CHOSEN TO CAUSE MAXIMUM LOCAL
%     GROWTH IN THE ELEMENTS OF W  WHERE  U*D*W = E .
%     THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW.
%
%     SOLVE U*D*W = E
%
ek =complex(1.0e0,0.0e0);
for j = 1 : n;
z(j) =complex(0.0e0,0.0e0);
end; j = fix(n+1);
k = fix(n);
while( k~=0 );
ks = 1;
if( kpvt(k)<0 )
ks = 2;
end;
kp = fix(abs(kpvt(k)));
kps = fix(k + 1 - ks);
if( kp~=kps )
t = z(kps);
z(kps) = z(kp);
z(kp) = t;
end;
csign1= @(zdum,zdum2)  cabs1(zdum).*(zdum2./cabs1(zdum2));
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(z(k))~=0.0e0 )
ek = csign1(ek,z(k));
end;
z(k) = z(k) + ek;
[dumvar1,dumvar2,a(sub2ind(size(a),1,k):end),dumvar4,dumvar5]=caxpy(k-ks,z(k),a(sub2ind(size(a),1,k):end),1,z(sub2ind(size(z),max(1,1)):end),1);   dumvar2i=find((z(k))~=(dumvar2));dumvar5i=find((z(sub2ind(size(z),max(1,1)):end))~=(dumvar5));   z(k-1+dumvar2i)=dumvar2(dumvar2i); z(1-1+dumvar5i)=dumvar5(dumvar5i);
if( ks~=1 )
csign1= @(zdum,zdum2)  cabs1(zdum).*(zdum2./cabs1(zdum2));
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(z(k-1))~=0.0e0 )
ek = csign1(ek,z(k-1));
end;
z(k-1) = z(k-1) + ek;
[dumvar1,dumvar2,a(sub2ind(size(a),1,k-1):end),dumvar4,dumvar5]=caxpy(k-ks,z(k-1),a(sub2ind(size(a),1,k-1):end),1,z(sub2ind(size(z),max(1,1)):end),1);   dumvar2i=find((z(k-1))~=(dumvar2));dumvar5i=find((z(sub2ind(size(z),max(1,1)):end))~=(dumvar5));   z(k-1-1+dumvar2i)=dumvar2(dumvar2i); z(1-1+dumvar5i)=dumvar5(dumvar5i);
end;
if( ks==2 )
ak = a(k,k)./conj(a(k-1,k));
akm1 = a(k-1,k-1)./a(k-1,k);
bk = z(k)./conj(a(k-1,k));
bkm1 = z(k-1)./a(k-1,k);
denom = ak.*akm1 - 1.0e0;
z(k) =(akm1.*bk-bkm1)./denom;
z(k-1) =(ak.*bkm1-bk)./denom;
else;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(z(k))>cabs1(a(k,k)) )
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
s = cabs1(a(k,k))./cabs1(z(k));
[n,s,z]=csscal(n,s,z,1);
ek = complex(s,0.0e0).*ek;
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(a(k,k))~=0.0e0 )
z(k) = z(k)./a(k,k);
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(a(k,k))==0.0e0 )
z(k) =complex(1.0e0,0.0e0);
end;
end;
k = fix(k - ks);
end;
s = 1.0e0./scasum(n,z,1);
[n,s,z]=csscal(n,s,z,1);
%
%     SOLVE CTRANS(U)*Y = W
%
k = 1;
while( k<=n );
ks = 1;
if( kpvt(k)<0 )
ks = 2;
end;
if( k~=1 )
z(k) = z(k) + cdotc(k-1,a(sub2ind(size(a),1,k):end),1,z(sub2ind(size(z),max(1,1)):end),1);
if( ks==2 )
z(k+1) = z(k+1) + cdotc(k-1,a(sub2ind(size(a),1,k+1):end),1,z(sub2ind(size(z),max(1,1)):end),1);
end;
kp = fix(abs(kpvt(k)));
if( kp~=k )
t = z(k);
z(k) = z(kp);
z(kp) = t;
end;
end;
k = fix(k + ks);
end;
s = 1.0e0./scasum(n,z,1);
[n,s,z]=csscal(n,s,z,1);
%
ynorm = 1.0e0;
%
%     SOLVE U*D*V = Y
%
k = fix(n);
while( k~=0 );
ks = 1;
if( kpvt(k)<0 )
ks = 2;
end;
if( k~=ks )
kp = fix(abs(kpvt(k)));
kps = fix(k + 1 - ks);
if( kp~=kps )
t = z(kps);
z(kps) = z(kp);
z(kp) = t;
end;
[dumvar1,dumvar2,a(sub2ind(size(a),1,k):end),dumvar4,dumvar5]=caxpy(k-ks,z(k),a(sub2ind(size(a),1,k):end),1,z(sub2ind(size(z),max(1,1)):end),1);   dumvar2i=find((z(k))~=(dumvar2));dumvar5i=find((z(sub2ind(size(z),max(1,1)):end))~=(dumvar5));   z(k-1+dumvar2i)=dumvar2(dumvar2i); z(1-1+dumvar5i)=dumvar5(dumvar5i);
if( ks==2 )
[dumvar1,dumvar2,a(sub2ind(size(a),1,k-1):end),dumvar4,dumvar5]=caxpy(k-ks,z(k-1),a(sub2ind(size(a),1,k-1):end),1,z(sub2ind(size(z),max(1,1)):end),1);   dumvar2i=find((z(k-1))~=(dumvar2));dumvar5i=find((z(sub2ind(size(z),max(1,1)):end))~=(dumvar5));   z(k-1-1+dumvar2i)=dumvar2(dumvar2i); z(1-1+dumvar5i)=dumvar5(dumvar5i);
end;
end;
if( ks==2 )
ak = a(k,k)./conj(a(k-1,k));
akm1 = a(k-1,k-1)./a(k-1,k);
bk = z(k)./conj(a(k-1,k));
bkm1 = z(k-1)./a(k-1,k);
denom = ak.*akm1 - 1.0e0;
z(k) =(akm1.*bk-bkm1)./denom;
z(k-1) =(ak.*bkm1-bk)./denom;
else;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(z(k))>cabs1(a(k,k)) )
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
s = cabs1(a(k,k))./cabs1(z(k));
[n,s,z]=csscal(n,s,z,1);
ynorm = s.*ynorm;
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(a(k,k))~=0.0e0 )
z(k) = z(k)./a(k,k);
end;
cabs1= @(zdum)  abs(real(zdum)) + abs(imag(zdum));
if( cabs1(a(k,k))==0.0e0 )
z(k) =complex(1.0e0,0.0e0);
end;
end;
k = fix(k - ks);
end;
s = 1.0e0./scasum(n,z,1);
[n,s,z]=csscal(n,s,z,1);
ynorm = s.*ynorm;
%
%     SOLVE CTRANS(U)*Z = V
%
k = 1;
while( k<=n );
ks = 1;
if( kpvt(k)<0 )
ks = 2;
end;
if( k~=1 )
z(k) = z(k) + cdotc(k-1,a(sub2ind(size(a),1,k):end),1,z(sub2ind(size(z),max(1,1)):end),1);
if( ks==2 )
z(k+1) = z(k+1) + cdotc(k-1,a(sub2ind(size(a),1,k+1):end),1,z(sub2ind(size(z),max(1,1)):end),1);
end;
kp = fix(abs(kpvt(k)));
if( kp~=k )
t = z(k);
z(k) = z(kp);
z(kp) = t;
end;
end;
k = fix(k + ks);
end;
%     MAKE ZNORM = 1.0
s = 1.0e0./scasum(n,z,1);
[n,s,z]=csscal(n,s,z,1);
ynorm = s.*ynorm;
%
if( anorm~=0.0e0 )
rcond = ynorm./anorm;
end;
if( anorm==0.0e0 )
rcond = 0.0e0;
end;
kpvt_shape=zeros(kpvt_shape);kpvt_shape(:)=kpvt(1:numel(kpvt_shape));kpvt=kpvt_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine chico
%DECK CHIDI
```