| [nm,n,a,b,alfr,alfi,beta,z]=qzvec(nm,n,a,b,alfr,alfi,beta,z); |
function [nm,n,a,b,alfr,alfi,beta,z]=qzvec(nm,n,a,b,alfr,alfi,beta,z);
%***BEGIN PROLOGUE QZVEC
%***PURPOSE The optional fourth step of the QZ algorithm for
% generalized eigenproblems. Accepts a matrix in
% quasi-triangular form and another in upper triangular
% and computes the eigenvectors of the triangular problem
% and transforms them back to the original coordinates
% Usually preceded by QZHES, QZIT, and QZVAL.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C3
%***TYPE SINGLE PRECISION (QZVEC-S)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is the optional fourth step of the QZ algorithm
% for solving generalized matrix eigenvalue problems,
% SIAM J. NUMER. ANAL. 10, 241-256(1973) by MOLER and STEWART.
%
% This subroutine accepts a pair of REAL matrices, one of them in
% quasi-triangular form (in which each 2-by-2 block corresponds to
% a pair of complex eigenvalues) and the other in upper triangular
% form. It computes the eigenvectors of the triangular problem and
% transforms the results back to the original coordinate system.
% It is usually preceded by QZHES, QZIT, and QZVAL.
%
% On Input
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, A, B, and Z, as declared in the calling
% program dimension statement. NM is an INTEGER variable.
%
% N is the order of the matrices A and B. N is an INTEGER
% variable. N must be less than or equal to NM.
%
% A contains a real upper quasi-triangular matrix. A is a two-
% dimensional REAL array, dimensioned A(NM,N).
%
% B contains a real upper triangular matrix. In addition,
% location B(N,1) contains the tolerance quantity (EPSB)
% computed and saved in QZIT. B is a two-dimensional REAL
% array, dimensioned B(NM,N).
%
% ALFR, ALFI, and BETA are one-dimensional REAL arrays with
% components whose ratios ((ALFR+I*ALFI)/BETA) are the
% generalized eigenvalues. They are usually obtained from
% QZVAL. They are dimensioned ALFR(N), ALFI(N), and BETA(N).
%
% Z contains the transformation matrix produced in the reductions
% by QZHES, QZIT, and QZVAL, if performed. If the
% eigenvectors of the triangular problem are desired, Z must
% contain the identity matrix. Z is a two-dimensional REAL
% array, dimensioned Z(NM,N).
%
% On Output
%
% A is unaltered. Its subdiagonal elements provide information
% about the storage of the complex eigenvectors.
%
% B has been destroyed.
%
% ALFR, ALFI, and BETA are unaltered.
%
% Z contains the real and imaginary parts of the eigenvectors.
% If ALFI(J) .EQ. 0.0, the J-th eigenvalue is real and
% the J-th column of Z contains its eigenvector.
% If ALFI(J) .NE. 0.0, the J-th eigenvalue is complex.
% If ALFI(J) .GT. 0.0, the eigenvalue is the first of
% a complex pair and the J-th and (J+1)-th columns
% of Z contain its eigenvector.
% If ALFI(J) .LT. 0.0, the eigenvalue is the second of
% a complex pair and the (J-1)-th and J-th columns
% of Z contain the conjugate of its eigenvector.
% Each eigenvector is normalized so that the modulus
% of its largest component is 1.0 .
%
% 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 (NONE)
%***REVISION HISTORY (YYMMDD)
% 760101 DATE WRITTEN
% 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 QZVEC
%
persistent alfm almi almr betm d di dr en enm2 epsb gt i ii isw j jj k m na nn q r ra rr s sa t t1 t2 ti tr w w1 x x1 y z1 zz ;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(m), m=0; end;
if isempty(en), en=0; end;
if isempty(ii), ii=0; end;
if isempty(jj), jj=0; end;
if isempty(na), na=0; end;
if isempty(nn), nn=0; end;
if isempty(isw), isw=0; end;
if isempty(enm2), enm2=0; end;
if isempty(gt), gt=zeros(1,8); end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([nm])).*prod([nm])-numel(b))],nm,[]);
alfr_shape=size(alfr);alfr=reshape(alfr,1,[]);
alfi_shape=size(alfi);alfi=reshape(alfi,1,[]);
beta_shape=size(beta);beta=reshape(beta,1,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(d), d=0; end;
if isempty(q), q=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
if isempty(w), w=0; end;
if isempty(x), x=0; end;
if isempty(y), y=0; end;
if isempty(di), di=0; end;
if isempty(dr), dr=0; end;
if isempty(ra), ra=0; end;
if isempty(rr), rr=0; end;
if isempty(sa), sa=0; end;
if isempty(ti), ti=0; end;
if isempty(tr), tr=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(w1), w1=0; end;
if isempty(x1), x1=0; end;
if isempty(zz), zz=0; end;
if isempty(z1), z1=0; end;
if isempty(alfm), alfm=0; end;
if isempty(almi), almi=0; end;
if isempty(almr), almr=0; end;
if isempty(betm), betm=0; end;
if isempty(epsb), epsb=0; end;
%
%***FIRST EXECUTABLE STATEMENT QZVEC
epsb = b(n,1);
isw = 1;
% .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
for nn = 1 : n;
en = fix(n + 1 - nn);
na = fix(en - 1);
if( isw~=2 )
if( alfi(en)~=0.0e0 )
% .......... COMPLEX VECTOR ..........
m = fix(na);
almr = alfr(m);
almi = alfi(m);
betm = beta(m);
% .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
% EIGENVECTOR MATRIX IS TRIANGULAR ..........
y = betm.*a(en,na);
b(na,na) = -almi.*b(en,en)./y;
b(na,en) =(almr.*b(en,en)-betm.*a(en,en))./y;
b(en,na) = 0.0e0;
b(en,en) = 1.0e0;
enm2 = fix(na - 1);
if( enm2~=0 )
% .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : enm2;
i = fix(na - ii);
w = betm.*a(i,i) - almr.*b(i,i);
w1 = -almi.*b(i,i);
ra = 0.0e0;
sa = 0.0e0;
%
for j = m : en;
x = betm.*a(i,j) - almr.*b(i,j);
x1 = -almi.*b(i,j);
ra = ra + x.*b(j,na) - x1.*b(j,en);
sa = sa + x.*b(j,en) + x1.*b(j,na);
end; j = fix(en+1);
%
if( i~=1 && isw~=2 )
if( betm.*a(i,i-1)~=0.0e0 )
zz = w;
z1 = w1;
r = ra;
s = sa;
isw = 2;
continue;
end;
end;
m = fix(i);
gt(:)=0;
while (1);
if(gt(8)==0)
if(gt(6)==0)
if(gt(4)==0)
if(gt(2)==0)
if( isw==2 )
gt(6)=1;
continue;
end;
% .......... COMPLEX 1-BY-1 BLOCK ..........
tr = -ra;
ti = -sa;
end;
gt(2)=0;
dr = w;
di = w1;
% .......... COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR,DI) ..........
end;
gt(4)=0;
if( abs(di)<=abs(dr) )
rr = di./dr;
d = dr + di.*rr;
t1 =(tr+ti.*rr)./d;
t2 =(ti-tr.*rr)./d;
if( isw==1 )
break;
end;
if( isw==2 )
gt(8)=1;
continue;
end;
end;
rr = dr./di;
d = dr.*rr + di;
t1 =(tr.*rr+ti)./d;
t2 =(ti.*rr-tr)./d;
if( isw==1 )
break;
end;
if( isw==2 )
gt(8)=1;
continue;
end;
% .......... COMPLEX 2-BY-2 BLOCK ..........
end;
gt(6)=0;
x = betm.*a(i,i+1) - almr.*b(i,i+1);
x1 = -almi.*b(i,i+1);
y = betm.*a(i+1,i);
tr = y.*ra - w.*r + w1.*s;
ti = y.*sa - w.*s - w1.*r;
dr = w.*zz - w1.*z1 - x.*y;
di = w.*z1 + w1.*zz - x1.*y;
if( dr==0.0e0 && di==0.0e0 )
dr = epsb;
end;
gt(4)=1;
continue;
end;
gt(8)=0;
b(i+1,na) = t1;
b(i+1,en) = t2;
isw = 1;
if( abs(y)>abs(w)+abs(w1) )
t1 =(-r-zz.*b(i+1,na)+z1.*b(i+1,en))./y;
t2 =(-s-zz.*b(i+1,en)-z1.*b(i+1,na))./y;
else;
tr = -ra - x.*b(i+1,na) + x1.*b(i+1,en);
ti = -sa - x.*b(i+1,en) - x1.*b(i+1,na);
gt(2)=1;
continue;
end;
break;
end;
b(i,na) = t1;
b(i,en) = t2;
end;
end;
else;
% .......... REAL VECTOR ..........
m = fix(en);
b(en,en) = 1.0e0;
if( na~=0 )
alfm = alfr(m);
betm = beta(m);
% .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : na;
i = fix(en - ii);
w = betm.*a(i,i) - alfm.*b(i,i);
r = 0.0e0;
%
for j = m : en;
r = r +(betm.*a(i,j)-alfm.*b(i,j)).*b(j,en);
end; j = fix(en+1);
%
while (1);
if( i~=1 && isw~=2 )
if( betm.*a(i,i-1)~=0.0e0 )
zz = w;
s = r;
break;
end;
end;
m = fix(i);
if( isw==2 )
% .......... REAL 2-BY-2 BLOCK ..........
x = betm.*a(i,i+1) - alfm.*b(i,i+1);
y = betm.*a(i+1,i);
q = w.*zz - x.*y;
t =(x.*s-zz.*r)./q;
b(i,en) = t;
if( abs(x)<=abs(zz) )
b(i+1,en) =(-s-y.*t)./zz;
else;
b(i+1,en) =(-r-w.*t)./x;
end;
else;
% .......... REAL 1-BY-1 BLOCK ..........
t = w;
if( w==0.0e0 )
t = epsb;
end;
b(i,en) = -r./t;
continue;
end;
break;
end;
isw = fix(3 - isw);
% .......... end REAL VECTOR ..........
end;
end;
continue;
end;
end;
% .......... end COMPLEX VECTOR ..........
isw = fix(3 - isw);
end;
% .......... end BACK SUBSTITUTION.
% TRANSFORM TO ORIGINAL COORDINATE SYSTEM.
% FOR J=N STEP -1 UNTIL 1 DO -- ..........
for jj = 1 : n;
j = fix(n + 1 - jj);
%
for i = 1 : n;
zz = 0.0e0;
%
for k = 1 : j;
zz = zz + z(i,k).*b(k,j);
end; k = fix(j+1);
%
z(i,j) = zz;
end; i = fix(n+1);
end; jj = fix(n+1);
% .......... NORMALIZE SO THAT MODULUS OF LARGEST
% COMPONENT OF EACH VECTOR IS 1.
% (ISW IS 1 INITIALLY FROM BEFORE) ..........
for j = 1 : n;
d = 0.0e0;
if( isw==2 )
%
for i = 1 : n;
r = abs(z(i,j-1)) + abs(z(i,j));
if( r~=0.0e0 )
r = r.*sqrt((z(i,j-1)./r).^2+(z(i,j)./r).^2);
end;
if( r>d )
d = r;
end;
end; i = fix(n+1);
%
for i = 1 : n;
z(i,j-1) = z(i,j-1)./d;
z(i,j) = z(i,j)./d;
end; i = fix(n+1);
elseif( alfi(j)==0.0e0 ) ;
%
for i = 1 : n;
if( abs(z(i,j))>d )
d = abs(z(i,j));
end;
end; i = fix(n+1);
%
for i = 1 : n;
z(i,j) = z(i,j)./d;
end; i = fix(n+1);
%
continue;
end;
%
isw = fix(3 - isw);
end; j = fix(n+1);
%
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
alfr_shape=zeros(alfr_shape);alfr_shape(:)=alfr(1:numel(alfr_shape));alfr=alfr_shape;
alfi_shape=zeros(alfi_shape);alfi_shape(:)=alfi(1:numel(alfi_shape));alfi=alfi_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine qzvec
|
|
