| [nm,n,low,igh,h,wr,wi,z,ierr]=hqr2(nm,n,low,igh,h,wr,wi,z,ierr); |
function [nm,n,low,igh,h,wr,wi,z,ierr]=hqr2(nm,n,low,igh,h,wr,wi,z,ierr);
%***BEGIN PROLOGUE HQR2
%***PURPOSE Compute the eigenvalues and eigenvectors of a real upper
% Hessenberg matrix using QR method.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C2B
%***TYPE SINGLE PRECISION (HQR2-S, COMQR2-C)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure HQR2,
% NUM. MATH. 16, 181-204(1970) by Peters and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
%
% This subroutine finds the eigenvalues and eigenvectors
% of a REAL UPPER Hessenberg matrix by the QR method. The
% eigenvectors of a REAL GENERAL matrix can also be found
% if ELMHES and ELTRAN or ORTHES and ORTRAN have
% been used to reduce this general matrix to Hessenberg form
% and to accumulate the similarity transformations.
%
% On INPUT
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, H and Z, as declared in the calling
% program dimension statement. NM is an INTEGER variable.
%
% N is the order of the matrix H. N is an INTEGER variable.
% N must be less than or equal to NM.
%
% LOW and IGH are two INTEGER variables determined by the
% balancing subroutine BALANC. If BALANC has not been
% used, set LOW=1 and IGH equal to the order of the matrix, N.
%
% H contains the upper Hessenberg matrix. H is a two-dimensional
% REAL array, dimensioned H(NM,N).
%
% Z contains the transformation matrix produced by ELTRAN
% after the reduction by ELMHES, or by ORTRAN after the
% reduction by ORTHES, if performed. If the eigenvectors
% of the Hessenberg matrix are desired, Z must contain the
% identity matrix. Z is a two-dimensional REAL array,
% dimensioned Z(NM,M).
%
% On OUTPUT
%
% H has been destroyed.
%
% WR and WI contain the real and imaginary parts, respectively,
% of the eigenvalues. The eigenvalues are unordered except
% that complex conjugate pairs of values appear consecutively
% with the eigenvalue having the positive imaginary part first.
% If an error exit is made, the eigenvalues should be correct
% for indices IERR+1, IERR+2, ..., N. WR and WI are one-
% dimensional REAL arrays, dimensioned WR(N) and WI(N).
%
% Z contains the real and imaginary parts of the eigenvectors.
% If the J-th eigenvalue is real, the J-th column of Z
% contains its eigenvector. If the J-th eigenvalue is complex
% with positive imaginary part, the J-th and (J+1)-th
% columns of Z contain the real and imaginary parts of its
% eigenvector. The eigenvectors are unnormalized. If an
% error exit is made, none of the eigenvectors has been found.
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% J if the J-th eigenvalue has not been
% determined after a total of 30*N iterations.
% The eigenvalues should be correct for indices
% IERR+1, IERR+2, ..., N, but no eigenvectors are
% computed.
%
% 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
%***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 HQR2
%
persistent en enm2 gt i ii itn its j jj k l ll m mm mp2 na nn norm notlas p q r ra s s1 s2 sa t vi vr w x y zz ;
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(m), m=0; end;
if isempty(en), en=0; end;
if isempty(ii), ii=0; end;
if isempty(jj), jj=0; end;
if isempty(ll), ll=0; end;
if isempty(mm), mm=0; end;
if isempty(na), na=0; end;
if isempty(nn), nn=0; end;
if isempty(gt), gt=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(mp2), mp2=0; end;
if isempty(enm2), enm2=0; end;
h_shape=size(h);h=reshape([h(:).',zeros(1,ceil(numel(h)./prod([nm])).*prod([nm])-numel(h))],nm,[]);
wr_shape=size(wr);wr=reshape(wr,1,[]);
wi_shape=size(wi);wi=reshape(wi,1,[]);
z_shape=size(z);z=reshape([z(:).',zeros(1,ceil(numel(z)./prod([nm])).*prod([nm])-numel(z))],nm,[]);
if isempty(p), p=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(ra), ra=0; end;
if isempty(sa), sa=0; end;
if isempty(vi), vi=0; end;
if isempty(vr), vr=0; end;
if isempty(zz), zz=0; end;
if isempty(norm), norm=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(notlas), notlas=false; end;
%
%***FIRST EXECUTABLE STATEMENT HQR2
gt=0;
ierr = 0;
norm = 0.0e0;
k = 1;
% .......... STORE ROOTS ISOLATED BY BALANC
% AND COMPUTE MATRIX NORM ..........
for i = 1 : n;
%
for j = k : n;
norm = norm + abs(h(i,j));
end; j = fix(n+1);
%
k = fix(i);
if( i<low || i>igh )
wr(i) = h(i,i);
wi(i) = 0.0e0;
end;
end; i = fix(n+1);
%
en = fix(igh);
t = 0.0e0;
itn = fix(30.*n);
% .......... SEARCH FOR NEXT EIGENVALUES ..........
while( en>=low );
gt=0;
its = 0;
na = fix(en - 1);
enm2 = fix(na - 1);
% .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
% FOR L=EN STEP -1 UNTIL LOW DO -- ..........
while( true );
for ll = low : en;
l = fix(en + low - ll);
if( l==low )
break;
end;
s = abs(h(l-1,l-1)) + abs(h(l,l));
if( s==0.0e0 )
s = norm;
end;
s2 = s + abs(h(l,l-1));
if( s2==s )
break;
end;
end;
% .......... FORM SHIFT ..........
x = h(en,en);
if( l==en )
% .......... ONE ROOT FOUND ..........
h(en,en) = x + t;
wr(en) = h(en,en);
wi(en) = 0.0e0;
en = fix(na);
gt=1;
break;
else;
y = h(na,na);
w = h(en,na).*h(na,en);
if( l==na )
break;
end;
if( itn==0 )
% .......... SET ERROR -- NO CONVERGENCE TO AN
% EIGENVALUE AFTER 30*N ITERATIONS ..........
ierr = fix(en);
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
return;
else;
if( its==10 || its==20 )
% .......... FORM EXCEPTIONAL SHIFT ..........
t = t + x;
%
for i = low : en;
h(i,i) = h(i,i) - x;
end; i = fix(en+1);
%
s = abs(h(en,na)) + abs(h(na,enm2));
x = 0.75e0.*s;
y = x;
w = -0.4375e0.*s.*s;
end;
its = fix(its + 1);
itn = fix(itn - 1);
% .......... LOOK FOR TWO CONSECUTIVE SMALL
% SUB-DIAGONAL ELEMENTS.
% FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
for mm = l : enm2;
m = fix(enm2 + l - mm);
zz = h(m,m);
r = x - zz;
s = y - zz;
p =(r.*s-w)./h(m+1,m) + h(m,m+1);
q = h(m+1,m+1) - zz - r - s;
r = h(m+2,m+1);
s = abs(p) + abs(q) + abs(r);
p = p./s;
q = q./s;
r = r./s;
if( m==l )
break;
end;
s1 = abs(p).*(abs(h(m-1,m-1))+abs(zz)+abs(h(m+1,m+1)));
s2 = s1 + abs(h(m,m-1)).*(abs(q)+abs(r));
if( s2==s1 )
break;
end;
end;
%
mp2 = fix(m + 2);
%
for i = mp2 : en;
h(i,i-2) = 0.0e0;
if( i~=mp2 )
h(i,i-3) = 0.0e0;
end;
end; i = fix(en+1);
% .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
% COLUMNS M TO EN ..........
for k = m : na;
notlas = k~=na;
if( k~=m )
p = h(k,k-1);
q = h(k+1,k-1);
r = 0.0e0;
if( notlas )
r = h(k+2,k-1);
end;
x = abs(p) + abs(q) + abs(r);
if( x==0.0e0 )
continue;
end;
p = p./x;
q = q./x;
r = r./x;
end;
s = (abs(sqrt(p.*p+q.*q+r.*r)).*sign(p));
if( k==m )
if( l~=m )
h(k,k-1) = -h(k,k-1);
end;
else;
h(k,k-1) = -s.*x;
end;
p = p + s;
x = p./s;
y = q./s;
zz = r./s;
q = q./p;
r = r./p;
% .......... ROW MODIFICATION ..........
for j = k : n;
p = h(k,j) + q.*h(k+1,j);
if( notlas )
p = p + r.*h(k+2,j);
h(k+2,j) = h(k+2,j) - p.*zz;
end;
h(k+1,j) = h(k+1,j) - p.*y;
h(k,j) = h(k,j) - p.*x;
end; j = fix(n+1);
%
j = fix(min(en,k+3));
% .......... COLUMN MODIFICATION ..........
for i = 1 : j;
p = x.*h(i,k) + y.*h(i,k+1);
if( notlas )
p = p + zz.*h(i,k+2);
h(i,k+2) = h(i,k+2) - p.*r;
end;
h(i,k+1) = h(i,k+1) - p.*q;
h(i,k) = h(i,k) - p;
end; i = fix(j+1);
% .......... ACCUMULATE TRANSFORMATIONS ..........
for i = low : igh;
%
p = x.*z(i,k) + y.*z(i,k+1);
if( notlas )
p = p + zz.*z(i,k+2);
z(i,k+2) = z(i,k+2) - p.*r;
end;
z(i,k+1) = z(i,k+1) - p.*q;
z(i,k) = z(i,k) - p;
%
end; i = fix(igh+1);
end; k = fix(na+1);
end;
end;
end;
if(gt==1)
continue;
end;
% .......... TWO ROOTS FOUND ..........
p =(y-x)./2.0e0;
q = p.*p + w;
zz = sqrt(abs(q));
h(en,en) = x + t;
x = h(en,en);
h(na,na) = y + t;
if( q<0.0e0 )
% .......... COMPLEX PAIR ..........
wr(na) = x + p;
wr(en) = x + p;
wi(na) = zz;
wi(en) = -zz;
else;
% .......... REAL PAIR ..........
zz = p + (abs(zz).*sign(p));
wr(na) = x + zz;
wr(en) = wr(na);
if( zz~=0.0e0 )
wr(en) = x - w./zz;
end;
wi(na) = 0.0e0;
wi(en) = 0.0e0;
x = h(en,na);
s = abs(x) + abs(zz);
p = x./s;
q = zz./s;
r = sqrt(p.*p+q.*q);
p = p./r;
q = q./r;
% .......... ROW MODIFICATION ..........
for j = na : n;
zz = h(na,j);
h(na,j) = q.*zz + p.*h(en,j);
h(en,j) = q.*h(en,j) - p.*zz;
end; j = fix(n+1);
% .......... COLUMN MODIFICATION ..........
for i = 1 : en;
zz = h(i,na);
h(i,na) = q.*zz + p.*h(i,en);
h(i,en) = q.*h(i,en) - p.*zz;
end; i = fix(en+1);
% .......... ACCUMULATE TRANSFORMATIONS ..........
for i = low : igh;
zz = z(i,na);
z(i,na) = q.*zz + p.*z(i,en);
z(i,en) = q.*z(i,en) - p.*zz;
%
end; i = fix(igh+1);
end;
en = fix(enm2);
end;
% .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
% VECTORS OF UPPER TRIANGULAR FORM ..........
if( norm~=0.0e0 )
% .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
for nn = 1 : n;
en = fix(n + 1 - nn);
p = wr(en);
q = wi(en);
na = fix(en - 1);
if( q<0 )
% .......... end COMPLEX VECTOR ..........
% .......... COMPLEX VECTOR ..........
m = fix(na);
% .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
% EIGENVECTOR MATRIX IS TRIANGULAR ..........
if( abs(h(en,na))<=abs(h(na,en)) )
[dumvar1,dumvar2,dumvar3,q,h(na,na),h(na,en)]=cdiv(0.0e0,-h(na,en),h(na,na)-p,q,h(na,na),h(na,en));
else;
h(na,na) = q./h(en,na);
h(na,en) = -(h(en,en)-p)./h(en,na);
end;
h(en,na) = 0.0e0;
h(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 = h(i,i) - p;
ra = 0.0e0;
sa = h(i,en);
%
for j = m : na;
ra = ra + h(i,j).*h(j,na);
sa = sa + h(i,j).*h(j,en);
end; j = fix(na+1);
%
if( wi(i)>=0.0e0 )
m = fix(i);
if( wi(i)==0.0e0 )
[dumvar1,dumvar2,w,q,h(i,na),h(i,en)]=cdiv(-ra,-sa,w,q,h(i,na),h(i,en));
else;
% .......... SOLVE COMPLEX EQUATIONS ..........
x = h(i,i+1);
y = h(i+1,i);
vr =(wr(i)-p).*(wr(i)-p) + wi(i).*wi(i) - q.*q;
vi =(wr(i)-p).*2.0e0.*q;
if( vr==0.0e0 && vi==0.0e0 )
s1 = norm.*(abs(w)+abs(q)+abs(x)+abs(y)+abs(zz));
vr = s1;
while( true );
vr = 0.5e0.*vr;
if( s1+vr<=s1 )
break;
end;
end;
vr = 2.0e0.*vr;
end;
[dumvar1,dumvar2,vr,vi,h(i,na),h(i,en)]=cdiv(x.*r-zz.*ra+q.*sa,x.*s-zz.*sa-q.*ra,vr,vi,h(i,na),h(i,en));
if( abs(x)<=abs(zz)+abs(q) )
[dumvar1,dumvar2,zz,q,h(i+1,na),h(i+1,en)]=cdiv(-r-y.*h(i,na),-s-y.*h(i,en),zz,q,h(i+1,na),h(i+1,en));
else;
h(i+1,na) =(-ra-w.*h(i,na)+q.*h(i,en))./x;
h(i+1,en) =(-sa-w.*h(i,en)-q.*h(i,na))./x;
end;
end;
else;
zz = w;
r = ra;
s = sa;
end;
end;
end;
elseif( q==0 ) ;
% .......... REAL VECTOR ..........
m = fix(en);
h(en,en) = 1.0e0;
if( na~=0 )
% .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : na;
i = fix(en - ii);
w = h(i,i) - p;
r = h(i,en);
if( m<=na )
%
for j = m : na;
r = r + h(i,j).*h(j,en);
end; j = fix(na+1);
end;
%
if( wi(i)>=0.0e0 )
% .......... end REAL VECTOR ..........
m = fix(i);
if( wi(i)==0.0e0 )
t = w;
if( t==0.0e0 )
t = norm;
while( true );
t = 0.5e0.*t;
if( norm+t<=norm )
break;
end;
end;
t = 2.0e0.*t;
end;
h(i,en) = -r./t;
else;
% .......... SOLVE REAL EQUATIONS ..........
x = h(i,i+1);
y = h(i+1,i);
q =(wr(i)-p).*(wr(i)-p) + wi(i).*wi(i);
t =(x.*s-zz.*r)./q;
h(i,en) = t;
if( abs(x)<=abs(zz) )
h(i+1,en) =(-s-y.*t)./zz;
else;
h(i+1,en) =(-r-w.*t)./x;
end;
end;
else;
zz = w;
s = r;
end;
end;
end;
end;
end;
% .......... end BACK SUBSTITUTION.
% VECTORS OF ISOLATED ROOTS ..........
for i = 1 : n;
if( i<low || i>igh )
%
%
for j = i : n;
z(i,j) = h(i,j);
end; j = fix(n+1);
end;
end; i = fix(n+1);
% .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
% VECTORS OF ORIGINAL FULL MATRIX.
% FOR J=N STEP -1 UNTIL LOW DO -- ..........
for jj = low : n;
j = fix(n + low - jj);
m = fix(min(j,igh));
%
for i = low : igh;
%
zz = 0.0e0;
%
for k = low : m;
zz = zz + z(i,k).*h(k,j);
end; k = fix(m+1);
%
z(i,j) = zz;
end; i = fix(igh+1);
end; jj = fix(n+1);
end;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_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;
z_shape=zeros(z_shape);z_shape(:)=z(1:numel(z_shape));z=z_shape;
end %subroutine hqr2
%DECK HQR
|
|
