| [nm,n,a,b,eps1,matz,z,ierr]=qzit(nm,n,a,b,eps1,matz,z,ierr); |
function [nm,n,a,b,eps1,matz,z,ierr]=qzit(nm,n,a,b,eps1,matz,z,ierr);
%***BEGIN PROLOGUE QZIT
%***PURPOSE The second step of the QZ algorithm for generalized
% eigenproblems. Accepts an upper Hessenberg and an upper
% triangular matrix and reduces the former to
% quasi-triangular form while preserving the form of the
% latter. Usually preceded by QZHES and followed by QZVAL
% and QZVEC.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C1B3
%***TYPE SINGLE PRECISION (QZIT-S)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is the second step of the QZ algorithm
% for solving generalized matrix eigenvalue problems,
% SIAM J. NUMER. ANAL. 10, 241-256(1973) by MOLER and STEWART,
% as modified in technical note NASA TN D-7305(1973) by WARD.
%
% This subroutine accepts a pair of REAL matrices, one of them
% in upper Hessenberg form and the other in upper triangular form.
% It reduces the Hessenberg matrix to quasi-triangular form using
% orthogonal transformations while maintaining the triangular form
% of the other matrix. It is usually preceded by QZHES and
% followed by QZVAL and, possibly, QZVEC.
%
% 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 Hessenberg matrix. A is a two-
% dimensional REAL array, dimensioned A(NM,N).
%
% B contains a real upper triangular matrix. B is a two-
% dimensional REAL array, dimensioned B(NM,N).
%
% EPS1 is a tolerance used to determine negligible elements.
% EPS1 = 0.0 (or negative) may be input, in which case an
% element will be neglected only if it is less than roundoff
% error times the norm of its matrix. If the input EPS1 is
% positive, then an element will be considered negligible
% if it is less than EPS1 times the norm of its matrix. A
% positive value of EPS1 may result in faster execution,
% but less accurate results. EPS1 is a REAL variable.
%
% MATZ should be set to true if the right hand transformations
% are to be accumulated for later use in computing
% eigenvectors, and to false otherwise. MATZ is a LOGICAL
% variable.
%
% Z contains, if MATZ has been set to true, the transformation
% matrix produced in the reduction by QZHES, if performed, or
% else the identity matrix. If MATZ has been set to false,
% Z is not referenced. Z is a two-dimensional REAL array,
% dimensioned Z(NM,N).
%
% On Output
%
% A has been reduced to quasi-triangular form. The elements
% below the first subdiagonal are still zero, and no two
% consecutive subdiagonal elements are nonzero.
%
% B is still in upper triangular form, although its elements
% have been altered. The location B(N,1) is used to store
% EPS1 times the norm of B for later use by QZVAL and QZVEC.
%
% Z contains the product of the right hand transformations
% (for both steps) if MATZ has been set to true
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% J if neither A(J,J-1) nor A(J-1,J-2) has become
% zero after a total of 30*N iterations.
%
% 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
% 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 QZIT
%
persistent a1 a11 a12 a2 a21 a22 a3 a33 a34 a43 a44 ani anorm b11 b12 b22 b33 b34 b44 bni bnorm en enm2 enorn ep epsa epsb gt i ish itn its j k k1 k2 km1 l l1 ld ll lm1 lor1 na notlas r s sh t u1 u2 u3 v1 v2 v3 ;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(en), en=0; end;
if isempty(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(ld), ld=0; end;
if isempty(ll), ll=0; end;
if isempty(l1), l1=0; end;
if isempty(na), na=0; end;
if isempty(ish), ish=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(km1), km1=0; end;
if isempty(lm1), lm1=0; end;
if isempty(gt), gt=zeros(1,5); end;
if isempty(enm2), enm2=0; end;
if isempty(lor1), lor1=0; end;
if isempty(enorn), enorn=0; 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,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(t), t=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(a3), a3=0; end;
if isempty(ep), ep=0; end;
if isempty(sh), sh=0; end;
if isempty(u1), u1=0; end;
if isempty(u2), u2=0; end;
if isempty(u3), u3=0; end;
if isempty(v1), v1=0; end;
if isempty(v2), v2=0; end;
if isempty(v3), v3=0; end;
if isempty(ani), ani=0; end;
if isempty(a11), a11=0; end;
if isempty(a12), a12=0; end;
if isempty(a21), a21=0; end;
if isempty(a22), a22=0; end;
if isempty(a33), a33=0; end;
if isempty(a34), a34=0; end;
if isempty(a43), a43=0; end;
if isempty(a44), a44=0; end;
if isempty(bni), bni=0; end;
if isempty(b11), b11=0; end;
if isempty(b12), b12=0; end;
if isempty(b22), b22=0; end;
if isempty(b33), b33=0; end;
if isempty(b34), b34=0; end;
if isempty(b44), b44=0; end;
if isempty(epsa), epsa=0; end;
if isempty(epsb), epsb=0; end;
if isempty(anorm), anorm=0; end;
if isempty(bnorm), bnorm=0; end;
if isempty(notlas), notlas=false; end;
%
%***FIRST EXECUTABLE STATEMENT QZIT
ierr = 0;
% .......... COMPUTE EPSA,EPSB ..........
anorm = 0.0e0;
bnorm = 0.0e0;
%
for i = 1 : n;
ani = 0.0e0;
if( i~=1 )
ani = abs(a(i,i-1));
end;
bni = 0.0e0;
%
for j = i : n;
ani = ani + abs(a(i,j));
bni = bni + abs(b(i,j));
end; j = fix(n+1);
%
if( ani>anorm )
anorm = ani;
end;
if( bni>bnorm )
bnorm = bni;
end;
end; i = fix(n+1);
%
if( anorm==0.0e0 )
anorm = 1.0e0;
end;
if( bnorm==0.0e0 )
bnorm = 1.0e0;
end;
ep = eps1;
if( ep<=0.0e0 )
% .......... COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO ..........
ep = 1.0e0;
while (1);
ep = ep./2.0e0;
if( 1.0e0+ep>1.0e0 )
continue;
end;
break;
end;
end;
epsa = ep.*anorm;
epsb = ep.*bnorm;
% .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE
% KEEPING B TRIANGULAR ..........
lor1 = 1;
enorn = fix(n);
en = fix(n);
itn = fix(30.*n);
% .......... BEGIN QZ STEP ..........
gt(:)=0;
while (1);
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
if( en<=2 )
% .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC ..........
if( n>1 )
b(n,1) = epsb;
end;
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
if( ~matz )
enorn = fix(en);
end;
its = 0;
na = fix(en - 1);
enm2 = fix(na - 1);
end;
end;
gt(2)=0;
ish = 2;
% .......... CHECK FOR CONVERGENCE OR REDUCIBILITY.
% FOR L=EN STEP -1 UNTIL 1 DO -- ..........
for ll = 1 : en;
lm1 = fix(en - ll);
l = fix(lm1 + 1);
if( l==1 )
% .......... CHECK FOR SMALL TOP OF B ..........
ld = fix(l);
gt(4)=1;
break;
elseif( abs(a(l,lm1))<=epsa ) ;
break;
end;
end;
if(gt(4)==1)
continue;
end;
%
end;
gt(3)=0;
a(l,lm1) = 0.0e0;
if( l<na )
ld = fix(l);
else;
% .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ..........
en = fix(lm1);
continue;
end;
end;
gt(4)=0;
l1 = fix(l + 1);
b11 = b(l,l);
if( abs(b11)>epsb )
a11 = a(l,l)./b11;
a21 = a(l1,l)./b11;
if( ish~=1 )
% .......... ITERATION STRATEGY ..........
if( itn==0 )
% .......... SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT
% HAS BECOME NEGLIGIBLE AFTER 30*N ITERATIONS ..........
ierr = fix(en);
if( n>1 )
b(n,1) = epsb;
end;
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
elseif( its==10 ) ;
% .......... AD HOC SHIFT ..........
a1 = 0.0e0;
a2 = 1.0e0;
a3 = 1.1605e0;
gt(5)=1;
continue;
else;
% .......... DETERMINE TYPE OF SHIFT ..........
b22 = b(l1,l1);
if( abs(b22)<epsb )
b22 = epsb;
end;
b33 = b(na,na);
if( abs(b33)<epsb )
b33 = epsb;
end;
b44 = b(en,en);
if( abs(b44)<epsb )
b44 = epsb;
end;
a33 = a(na,na)./b33;
a34 = a(na,en)./b44;
a43 = a(en,na)./b33;
a44 = a(en,en)./b44;
b34 = b(na,en)./b44;
t = 0.5e0.*(a43.*b34-a33-a44);
r = t.*t + a34.*a43 - a33.*a44;
if( r<0.0e0 )
% .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ..........
a12 = a(l,l1)./b22;
a22 = a(l1,l1)./b22;
b12 = b(l,l1)./b22;
a1 =((a33-a11).*(a44-a11)-a34.*a43+a43.*b34.*a11)./a21 + a12 -a11.*b12;
a2 =(a22-a11) - a21.*b12 -(a33-a11) -(a44-a11) + a43.*b34;
a3 = a(l1+1,l1)./b22;
gt(5)=1;
continue;
else;
% .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ..........
ish = 1;
r = sqrt(r);
sh = -t + r;
s = -t - r;
if( abs(s-a44)<abs(sh-a44) )
sh = s;
end;
% .......... LOOK FOR TWO CONSECUTIVE SMALL
% SUB-DIAGONAL ELEMENTS OF A.
% FOR L=EN-2 STEP -1 UNTIL LD DO -- ..........
for ll = ld : enm2;
l = fix(enm2 + ld - ll);
if( l==ld )
break;
end;
lm1 = fix(l - 1);
l1 = fix(l + 1);
t = a(l,l);
if( abs(b(l,l))>epsb )
t = t - sh.*b(l,l);
end;
if( abs(a(l,lm1))<=abs(t./a(l1,l)).*epsa )
gt(5)=1;
continue;
end;
end;
end;
end;
end;
%
a1 = a11 - sh;
a2 = a21;
if( l~=ld )
a(l,lm1) = -a(l,lm1);
end;
else;
b(l,l) = 0.0e0;
s = abs(a(l,l)) + abs(a(l1,l));
u1 = a(l,l)./s;
u2 = a(l1,l)./s;
r = (abs(sqrt(u1.*u1+u2.*u2)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
u2 = v2./v1;
%
for j = l : enorn;
t = a(l,j) + u2.*a(l1,j);
a(l,j) = a(l,j) + t.*v1;
a(l1,j) = a(l1,j) + t.*v2;
t = b(l,j) + u2.*b(l1,j);
b(l,j) = b(l,j) + t.*v1;
b(l1,j) = b(l1,j) + t.*v2;
end; j = fix(enorn+1);
%
if( l~=1 )
a(l,lm1) = -a(l,lm1);
end;
lm1 = fix(l);
l = fix(l1);
gt(3)=1;
continue;
end;
end;
gt(5)=0;
its = fix(its + 1);
itn = fix(itn - 1);
if( ~matz )
lor1 = fix(ld);
end;
% .......... MAIN LOOP ..........
for k = l : na;
notlas = k~=na & ish==2;
k1 = fix(k + 1);
k2 = fix(k + 2);
km1 = fix(max(k-1,l));
ll = fix(min(en,k1+ish));
if( notlas )
% .......... ZERO A(K+1,K-1) AND A(K+2,K-1) ..........
if( k~=l )
a1 = a(k,km1);
a2 = a(k1,km1);
a3 = a(k2,km1);
end;
s = abs(a1) + abs(a2) + abs(a3);
if( s==0.0e0 )
continue;
end;
u1 = a1./s;
u2 = a2./s;
u3 = a3./s;
r = (abs(sqrt(u1.*u1+u2.*u2+u3.*u3)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
v3 = -u3./r;
u2 = v2./v1;
u3 = v3./v1;
%
for j = km1 : enorn;
t = a(k,j) + u2.*a(k1,j) + u3.*a(k2,j);
a(k,j) = a(k,j) + t.*v1;
a(k1,j) = a(k1,j) + t.*v2;
a(k2,j) = a(k2,j) + t.*v3;
t = b(k,j) + u2.*b(k1,j) + u3.*b(k2,j);
b(k,j) = b(k,j) + t.*v1;
b(k1,j) = b(k1,j) + t.*v2;
b(k2,j) = b(k2,j) + t.*v3;
end; j = fix(enorn+1);
%
if( k~=l )
a(k1,km1) = 0.0e0;
a(k2,km1) = 0.0e0;
end;
% .......... ZERO B(K+2,K+1) AND B(K+2,K) ..........
s = abs(b(k2,k2)) + abs(b(k2,k1)) + abs(b(k2,k));
if( s~=0.0e0 )
u1 = b(k2,k2)./s;
u2 = b(k2,k1)./s;
u3 = b(k2,k)./s;
r = (abs(sqrt(u1.*u1+u2.*u2+u3.*u3)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
v3 = -u3./r;
u2 = v2./v1;
u3 = v3./v1;
%
for i = lor1 : ll;
t = a(i,k2) + u2.*a(i,k1) + u3.*a(i,k);
a(i,k2) = a(i,k2) + t.*v1;
a(i,k1) = a(i,k1) + t.*v2;
a(i,k) = a(i,k) + t.*v3;
t = b(i,k2) + u2.*b(i,k1) + u3.*b(i,k);
b(i,k2) = b(i,k2) + t.*v1;
b(i,k1) = b(i,k1) + t.*v2;
b(i,k) = b(i,k) + t.*v3;
end; i = fix(ll+1);
%
b(k2,k) = 0.0e0;
b(k2,k1) = 0.0e0;
if( matz )
%
for i = 1 : n;
t = z(i,k2) + u2.*z(i,k1) + u3.*z(i,k);
z(i,k2) = z(i,k2) + t.*v1;
z(i,k1) = z(i,k1) + t.*v2;
z(i,k) = z(i,k) + t.*v3;
end; i = fix(n+1);
end;
end;
else;
% .......... ZERO A(K+1,K-1) ..........
if( k~=l )
a1 = a(k,km1);
a2 = a(k1,km1);
end;
s = abs(a1) + abs(a2);
if( s==0.0e0 )
break;
end;
u1 = a1./s;
u2 = a2./s;
r = (abs(sqrt(u1.*u1+u2.*u2)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
u2 = v2./v1;
%
for j = km1 : enorn;
t = a(k,j) + u2.*a(k1,j);
a(k,j) = a(k,j) + t.*v1;
a(k1,j) = a(k1,j) + t.*v2;
t = b(k,j) + u2.*b(k1,j);
b(k,j) = b(k,j) + t.*v1;
b(k1,j) = b(k1,j) + t.*v2;
end; j = fix(enorn+1);
%
if( k~=l )
a(k1,km1) = 0.0e0;
end;
end;
% .......... ZERO B(K+1,K) ..........
s = abs(b(k1,k1)) + abs(b(k1,k));
if( s~=0.0e0 )
u1 = b(k1,k1)./s;
u2 = b(k1,k)./s;
r = (abs(sqrt(u1.*u1+u2.*u2)).*sign(u1));
v1 = -(u1+r)./r;
v2 = -u2./r;
u2 = v2./v1;
%
for i = lor1 : ll;
t = a(i,k1) + u2.*a(i,k);
a(i,k1) = a(i,k1) + t.*v1;
a(i,k) = a(i,k) + t.*v2;
t = b(i,k1) + u2.*b(i,k);
b(i,k1) = b(i,k1) + t.*v1;
b(i,k) = b(i,k) + t.*v2;
end; i = fix(ll+1);
%
b(k1,k) = 0.0e0;
if( matz )
%
for i = 1 : n;
t = z(i,k1) + u2.*z(i,k);
z(i,k1) = z(i,k1) + t.*v1;
z(i,k) = z(i,k) + t.*v2;
end; i = fix(n+1);
end;
end;
%
end;
% .......... end QZ STEP ..........
gt(2)=1;
continue;
end;
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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine qzit
%DECK QZVAL
|
|
