| [nm,n,ar,ai,wr,wi,select,mm,m,zr,zi,ierr,rm1,rm2,rv1,rv2]=cinvit(nm,n,ar,ai,wr,wi,select,mm,m,zr,zi,ierr,rm1,rm2,rv1,rv2); |
function [nm,n,ar,ai,wr,wi,select,mm,m,zr,zi,ierr,rm1,rm2,rv1,rv2]=cinvit(nm,n,ar,ai,wr,wi,select,mm,m,zr,zi,ierr,rm1,rm2,rv1,rv2);
%***BEGIN PROLOGUE CINVIT
%***PURPOSE Compute the eigenvectors of a complex upper Hessenberg
% associated with specified eigenvalues using inverse
% iteration.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C2B
%***TYPE COMPLEX (INVIT-S, CINVIT-C)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure CXINVIT
% by Peters and Wilkinson.
% HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971).
%
% This subroutine finds those eigenvectors of A COMPLEX UPPER
% Hessenberg matrix corresponding to specified eigenvalues,
% using inverse iteration.
%
% On INPUT
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, AR, AI, ZR and ZI, as declared in the
% calling program dimension statement. NM is an INTEGER
% variable.
%
% N is the order of the matrix A=(AR,AI). N is an INTEGER
% variable. N must be less than or equal to NM.
%
% AR and AI contain the real and imaginary parts, respectively,
% of the complex upper Hessenberg matrix. AR and AI are
% two-dimensional REAL arrays, dimensioned AR(NM,N)
% and AI(NM,N).
%
% WR and WI contain the real and imaginary parts, respectively,
% of the eigenvalues of the matrix. The eigenvalues must be
% stored in a manner identical to that of subroutine COMLR,
% which recognizes possible splitting of the matrix. WR and
% WI are one-dimensional REAL arrays, dimensioned WR(N) and
% WI(N).
%
% SELECT specifies the eigenvectors to be found. The
% eigenvector corresponding to the J-th eigenvalue is
% specified by setting SELECT(J) to true SELECT is a
% one-dimensional LOGICAL array, dimensioned SELECT(N).
%
% MM should be set to an upper bound for the number of
% eigenvectors to be found. MM is an INTEGER variable.
%
% On OUTPUT
%
% AR, AI, WI, and SELECT are unaltered.
%
% WR may have been altered since close eigenvalues are perturbed
% slightly in searching for independent eigenvectors.
%
% M is the number of eigenvectors actually found. M is an
% INTEGER variable.
%
% ZR and ZI contain the real and imaginary parts, respectively,
% of the eigenvectors corresponding to the flagged eigenvalues.
% The eigenvectors are normalized so that the component of
% largest magnitude is 1. Any vector which fails the
% acceptance test is set to zero. ZR and ZI are
% two-dimensional REAL arrays, dimensioned ZR(NM,MM) and
% ZI(NM,MM).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% -(2*N+1) if more than MM eigenvectors have been requested
% (the MM eigenvectors calculated to this point are
% in ZR and ZI),
% -K if the iteration corresponding to the K-th
% value fails (if this occurs more than once, K
% is the index of the last occurrence); the
% corresponding columns of ZR and ZI are set to
% zero vectors,
% -(N+K) if both error situations occur.
%
% RV1 and RV2 are one-dimensional REAL arrays used for
% temporary storage, dimensioned RV1(N) and RV2(N).
% They hold the approximate eigenvectors during the inverse
% iteration process.
%
% RM1 and RM2 are two-dimensional REAL arrays used for
% temporary storage, dimensioned RM1(N,N) and RM2(N,N).
% These arrays hold the triangularized form of the upper
% Hessenberg matrix used in the inverse iteration process.
%
% The ALGOL procedure GUESSVEC appears in CINVIT in-line.
%
% Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
% Calls CDIV for complex division.
%
% Questions and comments should be directed to B. S. Garbow,
% APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
% ------------------------------------------------------------------
%
%***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
% Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
% system Routines - EISPACK Guide, Springer-Verlag,
% 1976.
%***ROUTINES CALLED CDIV, PYTHAG
%***REVISION HISTORY (YYMMDD)
% 760101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890831 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE CINVIT
%
persistent eps3 growto gt10 i igo ii ilambd ip1 its j k km1 mp norm normv rlambd s uk ukroot x y ;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(s), s=0; end;
if isempty(ii), ii=0; end;
if isempty(mp), mp=0; end;
if isempty(uk), uk=0; end;
if isempty(ip1), ip1=0; end;
if isempty(its), its=0; end;
if isempty(km1), km1=0; end;
if isempty(igo), igo=0; end;
if isempty(gt10), gt10=0; end;
ar_shape=size(ar);ar=reshape([ar(:).',zeros(1,ceil(numel(ar)./prod([nm])).*prod([nm])-numel(ar))],nm,[]);
ai_shape=size(ai);ai=reshape([ai(:).',zeros(1,ceil(numel(ai)./prod([nm])).*prod([nm])-numel(ai))],nm,[]);
wr_shape=size(wr);wr=reshape(wr,1,[]);
wi_shape=size(wi);wi=reshape(wi,1,[]);
zr_shape=size(zr);zr=reshape([zr(:).',zeros(1,ceil(numel(zr)./prod([nm])).*prod([nm])-numel(zr))],nm,[]);
zi_shape=size(zi);zi=reshape([zi(:).',zeros(1,ceil(numel(zi)./prod([nm])).*prod([nm])-numel(zi))],nm,[]);
rm1_shape=size(rm1);rm1=reshape([rm1(:).',zeros(1,ceil(numel(rm1)./prod([n])).*prod([n])-numel(rm1))],n,[]);
rm2_shape=size(rm2);rm2=reshape([rm2(:).',zeros(1,ceil(numel(rm2)./prod([n])).*prod([n])-numel(rm2))],n,[]);
rv1_shape=size(rv1);rv1=reshape(rv1,1,[]);
rv2_shape=size(rv2);rv2=reshape(rv2,1,[]);
if isempty(x), x=0; end;
if isempty(y), y=0; end;
if isempty(eps3), eps3=0; end;
if isempty(norm), norm=0; end;
if isempty(normv), normv=0; end;
if isempty(growto), growto=0; end;
if isempty(ilambd), ilambd=0; end;
if isempty(rlambd), rlambd=0; end;
if isempty(ukroot), ukroot=0; end;
%
%***FIRST EXECUTABLE STATEMENT CINVIT
ierr = 0;
uk = 0;
s = 1;
%
for k = 1 : n;
if( select(k) )
%
if( s>mm )
% .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR
% SPACE REQUIRED ..........
if( ierr~=0 )
ierr = fix(ierr - n);
end;
if( ierr==0 )
ierr = fix(-(2.*n+1));
end;
break;
else;
if( uk<k )
% .......... CHECK FOR POSSIBLE SPLITTING ..........
for uk = k : n;
if( uk==n )
break;
end;
if( ar(uk+1,uk)==0.0e0 && ai(uk+1,uk)==0.0e0 )
break;
end;
end;
% .......... COMPUTE INFINITY NORM OF LEADING UK BY UK
% (HESSENBERG) MATRIX ..........
norm = 0.0e0;
mp = 1;
%
for i = 1 : uk;
x = 0.0e0;
%
for j = mp : uk;
x = x + pythag(ar(i,j),ai(i,j));
end; j = fix(uk+1);
%
if( x>norm )
norm = x;
end;
mp = fix(i);
end; i = fix(uk+1);
% .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION
% AND CLOSE ROOTS ARE MODIFIED BY EPS3 ..........
if( norm==0.0e0 )
norm = 1.0e0;
end;
eps3 = norm;
while( true );
eps3 = 0.5e0.*eps3;
if( norm+eps3<=norm )
break;
end;
end;
eps3 = 2.0e0.*eps3;
% .......... GROWTO IS THE CRITERION FOR GROWTH ..........
ukroot = sqrt(real(uk));
growto = 0.1e0./ukroot;
end;
rlambd = wr(k);
ilambd = wi(k);
if( k~=1 )
km1 = fix(k - 1);
% .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : km1;
i = fix(k - ii);
if( select(i) && abs(wr(i)-rlambd)<eps3 && abs(wi(i)-ilambd)<eps3 )
% .......... PERTURB EIGENVALUE IF IT IS CLOSE
% TO ANY PREVIOUS EIGENVALUE ..........
rlambd = rlambd + eps3;
continue;
end;
end; ii = fix(km1+1);
%
wr(k) = rlambd;
end;
% .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I
% AND INITIAL COMPLEX VECTOR ..........
mp = 1;
%
for i = 1 : uk;
%
for j = mp : uk;
rm1(i,j) = ar(i,j);
rm2(i,j) = ai(i,j);
end; j = fix(uk+1);
%
rm1(i,i) = rm1(i,i) - rlambd;
rm2(i,i) = rm2(i,i) - ilambd;
mp = fix(i);
rv1(i) = eps3;
end; i = fix(uk+1);
% .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
% REPLACING ZERO PIVOTS BY EPS3 ..........
if( uk~=1 )
%
for i = 2 : uk;
mp = fix(i - 1);
if( pythag(rm1(i,mp),rm2(i,mp))>pythag(rm1(mp,mp),rm2(mp,mp)) )
%
for j = mp : uk;
y = rm1(i,j);
rm1(i,j) = rm1(mp,j);
rm1(mp,j) = y;
y = rm2(i,j);
rm2(i,j) = rm2(mp,j);
rm2(mp,j) = y;
end; j = fix(uk+1);
end;
%
if( rm1(mp,mp)==0.0e0 && rm2(mp,mp)==0.0e0 )
rm1(mp,mp) = eps3;
end;
[rm1(i,mp),rm2(i,mp),rm1(mp,mp),rm2(mp,mp),x,y]=cdiv(rm1(i,mp),rm2(i,mp),rm1(mp,mp),rm2(mp,mp),x,y);
if( x~=0.0e0 || y~=0.0e0 )
%
%
for j = i : uk;
rm1(i,j) = rm1(i,j) - x.*rm1(mp,j) + y.*rm2(mp,j);
rm2(i,j) = rm2(i,j) - x.*rm2(mp,j) - y.*rm1(mp,j);
end; j = fix(uk+1);
end;
end; i = fix(uk+1);
end;
%
if( rm1(uk,uk)==0.0e0 && rm2(uk,uk)==0.0e0 )
rm1(uk,uk)= eps3;
end;
its = 0;
% .......... BACK SUBSTITUTION
% FOR I=UK STEP -1 UNTIL 1 DO -- ..........
igo=0;
while( true );
for ii = 1 : uk;
i = fix(uk + 1 - ii);
x = rv1(i);
y = 0.0e0;
if( i~=uk )
ip1 = fix(i + 1);
%
for j = ip1 : uk;
x = x - rm1(i,j).*rv1(j) + rm2(i,j).*rv2(j);
y = y - rm1(i,j).*rv2(j) - rm2(i,j).*rv1(j);
end; j = fix(uk+1);
end;
%
[x,y,rm1(i,i),rm2(i,i),rv1(i),rv2(i)]=cdiv(x,y,rm1(i,i),rm2(i,i),rv1(i),rv2(i));
end; ii = fix(uk+1);
% .......... ACCEPTANCE TEST FOR EIGENVECTOR
% AND NORMALIZATION ..........
its = fix(its + 1);
norm = 0.0e0;
normv = 0.0e0;
%
for i = 1 : uk;
[x ,rv1(i),rv2(i)]=pythag(rv1(i),rv2(i));
if( normv<x )
normv = x;
j = fix(i);
end;
norm = norm + x;
end; i = fix(uk+1);
%
gt10=0;
if( norm>=growto )
% .......... ACCEPT VECTOR ..........
x = rv1(j);
y = rv2(j);
%
for i = 1 : uk;
[rv1(i),rv2(i),x,y,zr(i,s),zi(i,s)]=cdiv(rv1(i),rv2(i),x,y,zr(i,s),zi(i,s));
end; i = fix(uk+1);
%
if( uk==n )
igo=1;
break;
end;
j = fix(uk + 1);
gt10=1;
break;
else;
if( its>=uk )
break;
end;
x = ukroot;
y = eps3./(x+1.0e0);
rv1(1) = eps3;
%
for i = 2 : uk;
rv1(i) = y;
end; i = fix(uk+1);
%
j = fix(uk - its + 1);
rv1(j) = rv1(j) - eps3.*x;
end;
end;
if(gt10==0)
if(igo==1)
s = fix(s + 1);
continue;
end;
% .......... SET ERROR -- UNACCEPTED EIGENVECTOR ..........
j = 1;
ierr = fix(-k);
end;
% .......... IN-LINE PROCEDURE FOR CHOOSING
% A NEW STARTING VECTOR ..........
% .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
for i = j : n;
zr(i,s) = 0.0e0;
zi(i,s) = 0.0e0;
end; i = fix(n+1);
end;
%
s = fix(s + 1);
end;
end;
m = fix(s - 1);
ar_shape=zeros(ar_shape);ar_shape(:)=ar(1:numel(ar_shape));ar=ar_shape;
ai_shape=zeros(ai_shape);ai_shape(:)=ai(1:numel(ai_shape));ai=ai_shape;
wr_shape=zeros(wr_shape);wr_shape(:)=wr(1:numel(wr_shape));wr=wr_shape;
wi_shape=zeros(wi_shape);wi_shape(:)=wi(1:numel(wi_shape));wi=wi_shape;
zr_shape=zeros(zr_shape);zr_shape(:)=zr(1:numel(zr_shape));zr=zr_shape;
zi_shape=zeros(zi_shape);zi_shape(:)=zi(1:numel(zi_shape));zi=zi_shape;
rm1_shape=zeros(rm1_shape);rm1_shape(:)=rm1(1:numel(rm1_shape));rm1=rm1_shape;
rm2_shape=zeros(rm2_shape);rm2_shape(:)=rm2(1:numel(rm2_shape));rm2=rm2_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
rv2_shape=zeros(rv2_shape);rv2_shape(:)=rv2(1:numel(rv2_shape));rv2=rv2_shape;
end %subroutine cinvit
%DECK CKSCL
|
|
