Code covered by the BSD License

Highlights fromslatec

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

[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
```