| [nm,n,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr]=comqr2(nm,n,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr); |
function [nm,n,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr]=comqr2(nm,n,low,igh,ortr,orti,hr,hi,wr,wi,zr,zi,ierr);
%***BEGIN PROLOGUE COMQR2
%***PURPOSE Compute the eigenvalues and eigenvectors of a complex upper
% Hessenberg matrix.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C2B
%***TYPE COMPLEX (HQR2-S, COMQR2-C)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation of a unitary analogue of the
% ALGOL procedure COMLR2, NUM. MATH. 16, 181-204(1970) by Peters
% and Wilkinson.
% HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
% The unitary analogue substitutes the QR algorithm of Francis
% (COMP. JOUR. 4, 332-345(1962)) for the LR algorithm.
%
% This subroutine finds the eigenvalues and eigenvectors
% of a COMPLEX UPPER Hessenberg matrix by the QR
% method. The eigenvectors of a COMPLEX GENERAL matrix
% can also be found if CORTH has been used to reduce
% this general matrix to Hessenberg form.
%
% On INPUT
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, HR, HI, ZR, and ZI, as declared in the
% calling program dimension statement. NM is an INTEGER
% variable.
%
% N is the order of the matrix H=(HR,HI). 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 CBAL. If CBAL has not been used,
% set LOW=1 and IGH equal to the order of the matrix, N.
%
% ORTR and ORTI contain information about the unitary trans-
% formations used in the reduction by CORTH, if performed.
% Only elements LOW through IGH are used. If the eigenvectors
% of the Hessenberg matrix are desired, set ORTR(J) and
% ORTI(J) to 0.0E0 for these elements. ORTR and ORTI are
% one-dimensional REAL arrays, dimensioned ORTR(IGH) and
% ORTI(IGH).
%
% HR and HI contain the real and imaginary parts, respectively,
% of the complex upper Hessenberg matrix. Their lower
% triangles below the subdiagonal contain information about
% the unitary transformations used in the reduction by CORTH,
% if performed. If the eigenvectors of the Hessenberg matrix
% are desired, these elements may be arbitrary. HR and HI
% are two-dimensional REAL arrays, dimensioned HR(NM,N) and
% HI(NM,N).
%
% On OUTPUT
%
% ORTR, ORTI, and the upper Hessenberg portions of HR and HI
% have been destroyed.
%
% WR and WI contain the real and imaginary parts, respectively,
% of the eigenvalues of the upper Hessenberg matrix. 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).
%
% ZR and ZI contain the real and imaginary parts, respectively,
% of the eigenvectors. The eigenvectors are unnormalized.
% If an error exit is made, none of the eigenvectors has been
% found. ZR and ZI are two-dimensional REAL arrays,
% dimensioned ZR(NM,N) and ZI(NM,N).
%
% 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 CSROOT for complex square root.
% 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, CSROOT, 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 COMQR2
%
persistent en enm1 i iend ii ip1 itn its j jj k l ll lp1 m nn norm s1 s2 si sr ti tr xi xr yi yr zzi zzr ;
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(nn), nn=0; end;
if isempty(ip1), ip1=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(lp1), lp1=0; end;
if isempty(enm1), enm1=0; end;
if isempty(iend), iend=0; end;
hr_shape=size(hr);hr=reshape([hr(:).',zeros(1,ceil(numel(hr)./prod([nm])).*prod([nm])-numel(hr))],nm,[]);
hi_shape=size(hi);hi=reshape([hi(:).',zeros(1,ceil(numel(hi)./prod([nm])).*prod([nm])-numel(hi))],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,[]);
ortr_shape=size(ortr);ortr=reshape(ortr,1,[]);
orti_shape=size(orti);orti=reshape(orti,1,[]);
if isempty(si), si=0; end;
if isempty(sr), sr=0; end;
if isempty(ti), ti=0; end;
if isempty(tr), tr=0; end;
if isempty(xi), xi=0; end;
if isempty(xr), xr=0; end;
if isempty(yi), yi=0; end;
if isempty(yr), yr=0; end;
if isempty(zzi), zzi=0; end;
if isempty(zzr), zzr=0; end;
if isempty(norm), norm=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
%
%***FIRST EXECUTABLE STATEMENT COMQR2
ierr = 0;
% .......... INITIALIZE EIGENVECTOR MATRIX ..........
for i = 1 : n;
%
for j = 1 : n;
zr(i,j) = 0.0e0;
zi(i,j) = 0.0e0;
if( i==j )
zr(i,j) = 1.0e0;
end;
end; j = fix(n+1);
end; i = fix(n+1);
% .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
% FROM THE INFORMATION LEFT BY CORTH ..........
iend = fix(igh - low - 1);
if( iend>=0 )
if( iend~=0 )
% .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
for ii = 1 : iend;
i = fix(igh - ii);
if( ortr(i)~=0.0e0 || orti(i)~=0.0e0 )
%
if( hr(i,i-1)~=0.0e0 || hi(i,i-1)~=0.0e0 )
% .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
norm = hr(i,i-1).*ortr(i) + hi(i,i-1).*orti(i);
ip1 = fix(i + 1);
%
for k = ip1 : igh;
ortr(k) = hr(k,i-1);
orti(k) = hi(k,i-1);
end; k = fix(igh+1);
%
for j = i : igh;
sr = 0.0e0;
si = 0.0e0;
%
for k = i : igh;
sr = sr + ortr(k).*zr(k,j) + orti(k).*zi(k,j);
si = si + ortr(k).*zi(k,j) - orti(k).*zr(k,j);
end; k = fix(igh+1);
%
sr = sr./norm;
si = si./norm;
%
for k = i : igh;
%
zr(k,j) = zr(k,j) + sr.*ortr(k) - si.*orti(k);
zi(k,j) = zi(k,j) + sr.*orti(k) + si.*ortr(k);
end; k = fix(igh+1);
end; j = fix(igh+1);
end;
end;
end; ii = fix(iend+1);
end;
% .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
l = fix(low + 1);
%
for i = l : igh;
ll = fix(min(i+1,igh));
if( hi(i,i-1)~=0.0e0 )
%
[norm ,hr(i,i-1),hi(i,i-1)]=pythag(hr(i,i-1),hi(i,i-1));
yr = hr(i,i-1)./norm;
yi = hi(i,i-1)./norm;
hr(i,i-1) = norm;
hi(i,i-1) = 0.0e0;
%
for j = i : n;
si = yr.*hi(i,j) - yi.*hr(i,j);
hr(i,j) = yr.*hr(i,j) + yi.*hi(i,j);
hi(i,j) = si;
end; j = fix(n+1);
%
for j = 1 : ll;
si = yr.*hi(j,i) + yi.*hr(j,i);
hr(j,i) = yr.*hr(j,i) - yi.*hi(j,i);
hi(j,i) = si;
end; j = fix(ll+1);
%
for j = low : igh;
si = yr.*zi(j,i) + yi.*zr(j,i);
zr(j,i) = yr.*zr(j,i) - yi.*zi(j,i);
zi(j,i) = si;
end; j = fix(igh+1);
end;
end; i = fix(igh+1);
end;
% .......... STORE ROOTS ISOLATED BY CBAL ..........
for i = 1 : n;
if( i<low || i>igh )
wr(i) = hr(i,i);
wi(i) = hi(i,i);
end;
end; i = fix(n+1);
%
en = fix(igh);
tr = 0.0e0;
ti = 0.0e0;
itn = fix(30.*n);
% .......... SEARCH FOR NEXT EIGENVALUE ..........
while( en>=low );
its = 0;
enm1 = fix(en - 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;
s1 = abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) + abs(hr(l,l))+ abs(hi(l,l));
s2 = s1 + abs(hr(l,l-1));
if( s2==s1 )
break;
end;
end;
% .......... FORM SHIFT ..........
if( l==en )
break;
end;
if( itn==0 )
% .......... SET ERROR -- NO CONVERGENCE TO AN
% EIGENVALUE AFTER 30*N ITERATIONS ..........
ierr = fix(en);
hr_shape=zeros(hr_shape);hr_shape(:)=hr(1:numel(hr_shape));hr=hr_shape;
hi_shape=zeros(hi_shape);hi_shape(:)=hi(1:numel(hi_shape));hi=hi_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;
ortr_shape=zeros(ortr_shape);ortr_shape(:)=ortr(1:numel(ortr_shape));ortr=ortr_shape;
orti_shape=zeros(orti_shape);orti_shape(:)=orti(1:numel(orti_shape));orti=orti_shape;
return;
else;
if( its==10 || its==20 )
% .......... FORM EXCEPTIONAL SHIFT ..........
sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2));
si = 0.0e0;
else;
sr = hr(en,en);
si = hi(en,en);
xr = hr(enm1,en).*hr(en,enm1);
xi = hi(enm1,en).*hr(en,enm1);
if( xr~=0.0e0 || xi~=0.0e0 )
yr =(hr(enm1,enm1)-sr)./2.0e0;
yi =(hi(enm1,enm1)-si)./2.0e0;
[dumvar1,dumvar2,zzr,zzi]=csroot(yr.^2-yi.^2+xr,2.0e0.*yr.*yi+xi,zzr,zzi);
if( yr.*zzr+yi.*zzi<0.0e0 )
zzr = -zzr;
zzi = -zzi;
end;
xr_orig=xr; xi_orig=xi; [xr,xi,dumvar3,dumvar4,dumvar5,dumvar6]=cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi); xi(dumvar6~=xi_orig)=dumvar6(dumvar6~=xi_orig); xr(dumvar5~=xr_orig)=dumvar5(dumvar5~=xr_orig);
sr = sr - xr;
si = si - xi;
end;
end;
%
for i = low : en;
hr(i,i) = hr(i,i) - sr;
hi(i,i) = hi(i,i) - si;
end; i = fix(en+1);
%
tr = tr + sr;
ti = ti + si;
its = fix(its + 1);
itn = fix(itn - 1);
% .......... REDUCE TO TRIANGLE (ROWS) ..........
lp1 = fix(l + 1);
%
for i = lp1 : en;
sr = hr(i,i-1);
hr(i,i-1) = 0.0e0;
[norm ,dumvar2,sr]=pythag(pythag(hr(i-1,i-1),hi(i-1,i-1)),sr);
xr = hr(i-1,i-1)./norm;
wr(i-1) = xr;
xi = hi(i-1,i-1)./norm;
wi(i-1) = xi;
hr(i-1,i-1) = norm;
hi(i-1,i-1) = 0.0e0;
hi(i,i-1) = sr./norm;
%
for j = i : n;
%
yr = hr(i-1,j);
yi = hi(i-1,j);
zzr = hr(i,j);
zzi = hi(i,j);
hr(i-1,j) = xr.*yr + xi.*yi + hi(i,i-1).*zzr;
hi(i-1,j) = xr.*yi - xi.*yr + hi(i,i-1).*zzi;
hr(i,j) = xr.*zzr - xi.*zzi - hi(i,i-1).*yr;
hi(i,j) = xr.*zzi + xi.*zzr - hi(i,i-1).*yi;
end; j = fix(n+1);
end; i = fix(en+1);
%
si = hi(en,en);
if( si~=0.0e0 )
[norm ,hr(en,en),si]=pythag(hr(en,en),si);
sr = hr(en,en)./norm;
si = si./norm;
hr(en,en) = norm;
hi(en,en) = 0.0e0;
if( en~=n )
ip1 = fix(en + 1);
%
for j = ip1 : n;
yr = hr(en,j);
yi = hi(en,j);
hr(en,j) = sr.*yr + si.*yi;
hi(en,j) = sr.*yi - si.*yr;
end; j = fix(n+1);
end;
end;
% .......... INVERSE OPERATION (COLUMNS) ..........
for j = lp1 : en;
xr = wr(j-1);
xi = wi(j-1);
%
for i = 1 : j;
yr = hr(i,j-1);
yi = 0.0e0;
zzr = hr(i,j);
zzi = hi(i,j);
if( i~=j )
yi = hi(i,j-1);
hi(i,j-1) = xr.*yi + xi.*yr + hi(j,j-1).*zzi;
end;
hr(i,j-1) = xr.*yr - xi.*yi + hi(j,j-1).*zzr;
hr(i,j) = xr.*zzr + xi.*zzi - hi(j,j-1).*yr;
hi(i,j) = xr.*zzi - xi.*zzr - hi(j,j-1).*yi;
end; i = fix(j+1);
%
for i = low : igh;
%
yr = zr(i,j-1);
yi = zi(i,j-1);
zzr = zr(i,j);
zzi = zi(i,j);
zr(i,j-1) = xr.*yr - xi.*yi + hi(j,j-1).*zzr;
zi(i,j-1) = xr.*yi + xi.*yr + hi(j,j-1).*zzi;
zr(i,j) = xr.*zzr + xi.*zzi - hi(j,j-1).*yr;
zi(i,j) = xr.*zzi - xi.*zzr - hi(j,j-1).*yi;
end; i = fix(igh+1);
end; j = fix(en+1);
%
if( si~=0.0e0 )
%
for i = 1 : en;
yr = hr(i,en);
yi = hi(i,en);
hr(i,en) = sr.*yr - si.*yi;
hi(i,en) = sr.*yi + si.*yr;
end; i = fix(en+1);
%
for i = low : igh;
yr = zr(i,en);
yi = zi(i,en);
zr(i,en) = sr.*yr - si.*yi;
zi(i,en) = sr.*yi + si.*yr;
%
end; i = fix(igh+1);
end;
end;
end;
% .......... A ROOT FOUND ..........
hr(en,en) = hr(en,en) + tr;
wr(en) = hr(en,en);
hi(en,en) = hi(en,en) + ti;
wi(en) = hi(en,en);
en = fix(enm1);
end;
% .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
% VECTORS OF UPPER TRIANGULAR FORM ..........
norm = 0.0e0;
%
for i = 1 : n;
%
for j = i : n;
norm = norm + abs(hr(i,j)) + abs(hi(i,j));
end; j = fix(n+1);
end; i = fix(n+1);
%
if( n~=1 && norm~=0.0e0 )
% .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
for nn = 2 : n;
en = fix(n + 2 - nn);
xr = wr(en);
xi = wi(en);
enm1 = fix(en - 1);
% .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : enm1;
%
i = fix(en - ii);
zzr = hr(i,en);
zzi = hi(i,en);
if( i~=enm1 )
ip1 = fix(i + 1);
%
for j = ip1 : enm1;
zzr = zzr + hr(i,j).*hr(j,en) - hi(i,j).*hi(j,en);
zzi = zzi + hr(i,j).*hi(j,en) + hi(i,j).*hr(j,en);
end; j = fix(enm1+1);
end;
%
yr = xr - wr(i);
yi = xi - wi(i);
if( yr==0.0e0 && yi==0.0e0 )
yr = norm;
while( true );
yr = 0.5e0.*yr;
if( norm+yr<=norm )
break;
end;
end;
yr = 2.0e0.*yr;
end;
[zzr,zzi,yr,yi,hr(i,en),hi(i,en)]=cdiv(zzr,zzi,yr,yi,hr(i,en),hi(i,en));
end;
end;
% .......... end BACKSUBSTITUTION ..........
enm1 = fix(n - 1);
% .......... VECTORS OF ISOLATED ROOTS ..........
for i = 1 : enm1;
if( i<low || i>igh )
%
ip1 = fix(i + 1);
%
for j = ip1 : n;
zr(i,j) = hr(i,j);
zi(i,j) = hi(i,j);
end; j = fix(n+1);
end;
end; i = fix(enm1+1);
% .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
% VECTORS OF ORIGINAL FULL MATRIX.
% FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
for jj = low : enm1;
j = fix(n + low - jj);
m = fix(min(j-1,igh));
%
for i = low : igh;
%
zzr = zr(i,j);
zzi = zi(i,j);
%
for k = low : m;
zzr = zzr + zr(i,k).*hr(k,j) - zi(i,k).*hi(k,j);
zzi = zzi + zr(i,k).*hi(k,j) + zi(i,k).*hr(k,j);
end; k = fix(m+1);
%
zr(i,j) = zzr;
zi(i,j) = zzi;
end; i = fix(igh+1);
end; jj = fix(enm1+1);
end;
hr_shape=zeros(hr_shape);hr_shape(:)=hr(1:numel(hr_shape));hr=hr_shape;
hi_shape=zeros(hi_shape);hi_shape(:)=hi(1:numel(hi_shape));hi=hi_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;
ortr_shape=zeros(ortr_shape);ortr_shape(:)=ortr(1:numel(ortr_shape));ortr=ortr_shape;
orti_shape=zeros(orti_shape);orti_shape(:)=orti(1:numel(orti_shape));orti=orti_shape;
end %subroutine comqr2
%DECK COMQR
|
|
