| [ap,n,kpvt,info]=sspfa(ap,n,kpvt,info); |
function [ap,n,kpvt,info]=sspfa(ap,n,kpvt,info);
%***BEGIN PROLOGUE SSPFA
%***PURPOSE Factor a real symmetric matrix stored in packed form by
% elimination with symmetric pivoting.
%***LIBRARY SLATEC (LINPACK)
%***CATEGORY D2B1A
%***TYPE SINGLE PRECISION (SSPFA-S, DSPFA-D, CHPFA-C, CSPFA-C)
%***KEYWORDS LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION, PACKED,
% SYMMETRIC
%***AUTHOR Bunch, J., (UCSD)
%***DESCRIPTION
%
% SSPFA factors a real symmetric matrix stored in
% packed form by elimination with symmetric pivoting.
%
% To solve A*X = B , follow SSPFA by SSPSL.
% To compute INVERSE(A)*C , follow SSPFA by SSPSL.
% To compute DETERMINANT(A) , follow SSPFA by SSPDI.
% To compute INERTIA(A) , follow SSPFA by SSPDI.
% To compute INVERSE(A) , follow SSPFA by SSPDI.
%
% On Entry
%
% AP REAL (N*(N+1)/2)
% the packed form of a symmetric matrix A . The
% columns of the upper triangle are stored sequentially
% in a one-dimensional array of length N*(N+1)/2 .
% See comments below for details.
%
% N INTEGER
% the order of the matrix A .
%
% Output
%
% AP a block diagonal matrix and the multipliers which
% were used to obtain it stored in packed form.
% The factorization can be written A = U*D*TRANS(U)
% where U is a product of permutation and unit
% upper triangular matrices , TRANS(U) is the
% transpose of U , and D is block diagonal
% with 1 by 1 and 2 by 2 blocks.
%
% KPVT INTEGER(N)
% an integer vector of pivot indices.
%
% INFO INTEGER
% = 0 normal value.
% = K if the K-th pivot block is singular. This is
% not an error condition for this subroutine,
% but it does indicate that SSPSL or SSPDI may
% divide by zero if called.
%
% Packed Storage
%
% The following program segment will pack the upper
% triangle of a symmetric matrix.
%
% K = 0
% DO 20 J = 1, N
% DO 10 I = 1, J
% K = K + 1
% AP(K) = A(I,J)
% 10 CONTINUE
% 20 CONTINUE
%
%***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
% Stewart, LINPACK Users' Guide, SIAM, 1979.
%***ROUTINES CALLED ISAMAX, SAXPY, SSWAP
%***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 SSPFA
persistent absakk ak akm1 alpha bk bkm1 colmax denom ij ik ikm1 im imax imaxp1 imim imj imk j jj jk jkm1 jmax jmim k kk km1 km1k km1km1 km2 kstep mulk mulkm1 rowmax swap t ;
kpvt_shape=size(kpvt);kpvt=reshape(kpvt,1,[]);
ap_shape=size(ap);ap=reshape(ap,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(mulk), mulk=0; end;
if isempty(mulkm1), mulkm1=0; end;
if isempty(t), t=0; end;
if isempty(absakk), absakk=0; end;
if isempty(alpha), alpha=0; end;
if isempty(colmax), colmax=0; end;
if isempty(rowmax), rowmax=0; end;
if isempty(ij), ij=0; end;
if isempty(ik), ik=0; end;
if isempty(ikm1), ikm1=0; end;
if isempty(im), im=0; end;
if isempty(imax), imax=0; end;
if isempty(imaxp1), imaxp1=0; end;
if isempty(imim), imim=0; end;
if isempty(imj), imj=0; end;
if isempty(imk), imk=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(jk), jk=0; end;
if isempty(jkm1), jkm1=0; end;
if isempty(jmax), jmax=0; end;
if isempty(jmim), jmim=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(km1), km1=0; end;
if isempty(km1k), km1k=0; end;
if isempty(km1km1), km1km1=0; end;
if isempty(km2), km2=0; end;
if isempty(kstep), kstep=0; end;
if isempty(swap), swap=false; end;
%***FIRST EXECUTABLE STATEMENT SSPFA
%
% INITIALIZE
%
% ALPHA IS USED IN CHOOSING PIVOT BLOCK SIZE.
%
alpha =(1.0e0+sqrt(17.0e0))./8.0e0;
%
info = 0;
%
% MAIN LOOP ON K, WHICH GOES FROM N TO 1.
%
k = fix(n);
ik =fix(fix((n.*(n-1))./2));
%
% LEAVE THE LOOP IF K=0 OR K=1.
%
while( k~=0 );
if( k<=1 )
kpvt(1) = 1;
if( ap(1)==0.0e0 )
info = 1;
end;
break;
else;
%
% THIS SECTION OF CODE DETERMINES THE KIND OF
% ELIMINATION TO BE PERFORMED. WHEN IT IS COMPLETED,
% KSTEP WILL BE SET TO THE SIZE OF THE PIVOT BLOCK, AND
% SWAP WILL BE SET TO true IF AN INTERCHANGE IS
% REQUIRED.
%
km1 = fix(k - 1);
kk = fix(ik + k);
absakk = abs(ap(kk));
%
% DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
% COLUMN K.
%
[imax ,dumvar2,ap(sub2ind(size(ap),max(ik+1,1)):end)]=isamax(k-1,ap(sub2ind(size(ap),max(ik+1,1)):end),1);
imk = fix(ik + imax);
colmax = abs(ap(imk));
if( absakk<alpha.*colmax )
%
% DETERMINE THE LARGEST OFF-DIAGONAL ELEMENT IN
% ROW IMAX.
%
rowmax = 0.0e0;
imaxp1 = fix(imax + 1);
im = fix(fix((imax.*(imax-1))./2));
imj = fix(im + 2.*imax);
for j = imaxp1 : k;
rowmax = max(rowmax,abs(ap(imj)));
imj = fix(imj + j);
end; j = fix(k+1);
if( imax~=1 )
[jmax ,dumvar2,ap(sub2ind(size(ap),max(im+1,1)):end)]=isamax(imax-1,ap(sub2ind(size(ap),max(im+1,1)):end),1);
jmim = fix(jmax + im);
rowmax = max(rowmax,abs(ap(jmim)));
end;
imim = fix(imax + im);
if( abs(ap(imim))>=alpha.*rowmax )
kstep = 1;
swap = true;
elseif( absakk<alpha.*colmax.*(colmax./rowmax) ) ;
kstep = 2;
swap = imax~=km1;
else;
kstep = 1;
swap = false;
end;
else;
kstep = 1;
swap = false;
end;
if( max(absakk,colmax)==0.0e0 )
%
% COLUMN K IS ZERO. SET INFO AND ITERATE THE LOOP.
%
kpvt(k) = fix(k);
info = fix(k);
elseif( kstep==2 ) ;
%
% 2 X 2 PIVOT BLOCK.
%
km1k = fix(ik + k - 1);
ikm1 = fix(ik -(k-1));
if( swap )
%
% PERFORM AN INTERCHANGE.
%
[imax,dumvar2,dumvar3,dumvar4]=sswap(imax,ap(sub2ind(size(ap),max(im+1,1)):end),1,ap(sub2ind(size(ap),max(ikm1+1,1)):end),1); dumvar2i=find((ap(sub2ind(size(ap),max(im+1,1)):end))~=(dumvar2));dumvar4i=find((ap(sub2ind(size(ap),max(ikm1+1,1)):end))~=(dumvar4)); ap(im+1-1+dumvar2i)=dumvar2(dumvar2i); ap(ikm1+1-1+dumvar4i)=dumvar4(dumvar4i);
imj = fix(ikm1 + imax);
for jj = imax : km1;
j = fix(km1 + imax - jj);
jkm1 = fix(ikm1 + j);
t = ap(jkm1);
ap(jkm1) = ap(imj);
ap(imj) = t;
imj = fix(imj -(j-1));
end; jj = fix(km1+1);
t = ap(km1k);
ap(km1k) = ap(imk);
ap(imk) = t;
end;
%
% PERFORM THE ELIMINATION.
%
km2 = fix(k - 2);
if( km2~=0 )
ak = ap(kk)./ap(km1k);
km1km1 = fix(ikm1 + k - 1);
akm1 = ap(km1km1)./ap(km1k);
denom = 1.0e0 - ak.*akm1;
ij = fix(ik -(k-1) -(k-2));
for jj = 1 : km2;
j = fix(km1 - jj);
jk = fix(ik + j);
bk = ap(jk)./ap(km1k);
jkm1 = fix(ikm1 + j);
bkm1 = ap(jkm1)./ap(km1k);
mulk =(akm1.*bk-bkm1)./denom;
mulkm1 =(ak.*bkm1-bk)./denom;
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ik+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1); dumvar3i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5)); ap(ik+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
t = mulkm1;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ikm1+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1); dumvar3i=find((ap(sub2ind(size(ap),max(ikm1+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5)); ap(ikm1+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
ap(jk) = mulk;
ap(jkm1) = mulkm1;
ij = fix(ij -(j-1));
end; jj = fix(km2+1);
end;
%
% SET THE PIVOT ARRAY.
%
kpvt(k) = fix(1 - k);
if( swap )
kpvt(k) = fix(-imax);
end;
kpvt(k-1) = fix(kpvt(k));
else;
%
% 1 X 1 PIVOT BLOCK.
%
if( swap )
%
% PERFORM AN INTERCHANGE.
%
[imax,dumvar2,dumvar3,dumvar4]=sswap(imax,ap(sub2ind(size(ap),max(im+1,1)):end),1,ap(sub2ind(size(ap),max(ik+1,1)):end),1); dumvar2i=find((ap(sub2ind(size(ap),max(im+1,1)):end))~=(dumvar2));dumvar4i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar4)); ap(im+1-1+dumvar2i)=dumvar2(dumvar2i); ap(ik+1-1+dumvar4i)=dumvar4(dumvar4i);
imj = fix(ik + imax);
for jj = imax : k;
j = fix(k + imax - jj);
jk = fix(ik + j);
t = ap(jk);
ap(jk) = ap(imj);
ap(imj) = t;
imj = fix(imj -(j-1));
end; jj = fix(k+1);
end;
%
% PERFORM THE ELIMINATION.
%
ij = fix(ik -(k-1));
for jj = 1 : km1;
j = fix(k - jj);
jk = fix(ik + j);
mulk = -ap(jk)./ap(kk);
t = mulk;
[j,t,dumvar3,dumvar4,dumvar5]=saxpy(j,t,ap(sub2ind(size(ap),max(ik+1,1)):end),1,ap(sub2ind(size(ap),max(ij+1,1)):end),1); dumvar3i=find((ap(sub2ind(size(ap),max(ik+1,1)):end))~=(dumvar3));dumvar5i=find((ap(sub2ind(size(ap),max(ij+1,1)):end))~=(dumvar5)); ap(ik+1-1+dumvar3i)=dumvar3(dumvar3i); ap(ij+1-1+dumvar5i)=dumvar5(dumvar5i);
ap(jk) = mulk;
ij = fix(ij -(j-1));
end; jj = fix(km1+1);
%
% SET THE PIVOT ARRAY.
%
kpvt(k) = fix(k);
if( swap )
kpvt(k) = fix(imax);
end;
end;
ik = fix(ik -(k-1));
if( kstep==2 )
ik = fix(ik -(k-2));
end;
k = fix(k - kstep);
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;
end %subroutine sspfa
%DECK SSPMV
|
|
