| [nm,n,low,igh,intmlv,hr,hi,wr,wi,zr,zi,ierr]=comlr2(nm,n,low,igh,intmlv,hr,hi,wr,wi,zr,zi,ierr); |
function [nm,n,low,igh,intmlv,hr,hi,wr,wi,zr,zi,ierr]=comlr2(nm,n,low,igh,intmlv,hr,hi,wr,wi,zr,zi,ierr);
%***BEGIN PROLOGUE COMLR2
%***PURPOSE Compute the eigenvalues and eigenvectors of a complex upper
% Hessenberg matrix using the modified LR method.
%***LIBRARY SLATEC (EISPACK)
%***CATEGORY D4C2B
%***TYPE COMPLEX (COMLR2-C)
%***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK, LR METHOD
%***AUTHOR Smith, B. T., et al.
%***DESCRIPTION
%
% This subroutine is a translation 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).
%
% This subroutine finds the eigenvalues and eigenvectors
% of a COMPLEX UPPER Hessenberg matrix by the modified LR
% method. The eigenvectors of a COMPLEX GENERAL matrix
% can also be found if COMHES 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.
%
% INT contains information on the rows and columns
% interchanged in the reduction by COMHES, if performed.
% Only elements LOW through IGH are used. If you want the
% eigenvectors of a complex general matrix, leave INT as it
% came from COMHES. If the eigenvectors of the Hessenberg
% matrix are desired, set INT(J)=J for these elements. INT
% is a one-dimensional INTEGER array, dimensioned INT(IGH).
%
% HR and HI contain the real and imaginary parts, respectively,
% of the complex upper Hessenberg matrix. Their lower
% triangles below the subdiagonal contain the multipliers
% which were used in the reduction by COMHES, if performed.
% If the eigenvectors of a complex general matrix are
% desired, leave these multipliers in the lower triangles.
% If the eigenvectors of the Hessenberg matrix are desired,
% these elements must be set to zero. HR and HI are
% two-dimensional REAL arrays, dimensioned HR(NM,N) and
% HI(NM,N).
%
% On OUTPUT
%
% The upper Hessenberg portions of HR and HI have been
% destroyed, but the location HR(1,1) contains the norm
% of the triangularized matrix.
%
% 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 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
%***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 COMLR2
%
persistent en enm1 i iend ii im1 ip1 itn its j jj k l ll m mm mp1 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(mm), mm=0; end;
if isempty(nn), nn=0; end;
if isempty(im1), im1=0; end;
if isempty(ip1), ip1=0; end;
if isempty(itn), itn=0; end;
if isempty(its), its=0; end;
if isempty(mp1), mp1=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,[]);
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;
int_shape=size(intmlv);intmlv=reshape(intmlv,1,[]);
%
%***FIRST EXECUTABLE STATEMENT COMLR2
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 COMHES ..........
iend = fix(igh - low - 1);
if( iend>0 )
% .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
for ii = 1 : iend;
i = fix(igh - ii);
ip1 = fix(i + 1);
%
for k = ip1 : igh;
zr(k,i) = hr(k,i-1);
zi(k,i) = hi(k,i-1);
end; k = fix(igh+1);
%
j = fix(intmlv(i));
if( i~=j )
%
for k = i : igh;
zr(i,k) = zr(j,k);
zi(i,k) = zi(j,k);
zr(j,k) = 0.0e0;
zi(j,k) = 0.0e0;
end; k = fix(igh+1);
%
zr(j,i) = 1.0e0;
end;
end; ii = fix(iend+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)) + abs(hi(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;
int_shape=zeros(int_shape);int_shape(:)=intmlv(1:numel(int_shape));intmlv=int_shape;
return;
else;
if( its==10 || its==20 )
% .......... FORM EXCEPTIONAL SHIFT ..........
sr = abs(hr(en,enm1)) + abs(hr(enm1,en-2));
si = abs(hi(en,enm1)) + abs(hi(enm1,en-2));
else;
sr = hr(en,en);
si = hi(en,en);
xr = hr(enm1,en).*hr(en,enm1) - hi(enm1,en).*hi(en,enm1);
xi = hr(enm1,en).*hi(en,enm1) + 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);
% .......... LOOK FOR TWO CONSECUTIVE SMALL
% SUB-DIAGONAL ELEMENTS ..........
xr = abs(hr(enm1,enm1)) + abs(hi(enm1,enm1));
yr = abs(hr(en,enm1)) + abs(hi(en,enm1));
zzr = abs(hr(en,en)) + abs(hi(en,en));
% .......... FOR M=EN-1 STEP -1 UNTIL L DO -- ..........
for mm = l : enm1;
m = fix(enm1 + l - mm);
if( m==l )
break;
end;
yi = yr;
yr = abs(hr(m,m-1)) + abs(hi(m,m-1));
xi = zzr;
zzr = xr;
xr = abs(hr(m-1,m-1)) + abs(hi(m-1,m-1));
s1 = zzr./yi.*(zzr+xr+xi);
s2 = s1 + yr;
if( s2==s1 )
break;
end;
end;
% .......... TRIANGULAR DECOMPOSITION H=L*R ..........
mp1 = fix(m + 1);
%
for i = mp1 : en;
im1 = fix(i - 1);
xr = hr(im1,im1);
xi = hi(im1,im1);
yr = hr(i,im1);
yi = hi(i,im1);
if( abs(xr)+abs(xi)>=abs(yr)+abs(yi) )
[yr,yi,xr,xi,zzr,zzi]=cdiv(yr,yi,xr,xi,zzr,zzi);
wr(i) = -1.0e0;
else;
% .......... INTERCHANGE ROWS OF HR AND HI ..........
for j = im1 : n;
zzr = hr(im1,j);
hr(im1,j) = hr(i,j);
hr(i,j) = zzr;
zzi = hi(im1,j);
hi(im1,j) = hi(i,j);
hi(i,j) = zzi;
end; j = fix(n+1);
%
[xr,xi,yr,yi,zzr,zzi]=cdiv(xr,xi,yr,yi,zzr,zzi);
wr(i) = 1.0e0;
end;
hr(i,im1) = zzr;
hi(i,im1) = zzi;
%
for j = i : n;
%
hr(i,j) = hr(i,j) - zzr.*hr(im1,j) + zzi.*hi(im1,j);
hi(i,j) = hi(i,j) - zzr.*hi(im1,j) - zzi.*hr(im1,j);
end; j = fix(n+1);
end; i = fix(en+1);
% .......... COMPOSITION R*L=H ..........
for j = mp1 : en;
xr = hr(j,j-1);
xi = hi(j,j-1);
hr(j,j-1) = 0.0e0;
hi(j,j-1) = 0.0e0;
% .......... INTERCHANGE COLUMNS OF HR, HI, ZR, AND ZI,
% IF NECESSARY ..........
if( wr(j)>0.0e0 )
%
for i = 1 : j;
zzr = hr(i,j-1);
hr(i,j-1) = hr(i,j);
hr(i,j) = zzr;
zzi = hi(i,j-1);
hi(i,j-1) = hi(i,j);
hi(i,j) = zzi;
end; i = fix(j+1);
%
for i = low : igh;
zzr = zr(i,j-1);
zr(i,j-1) = zr(i,j);
zr(i,j) = zzr;
zzi = zi(i,j-1);
zi(i,j-1) = zi(i,j);
zi(i,j) = zzi;
end; i = fix(igh+1);
end;
%
for i = 1 : j;
hr(i,j-1) = hr(i,j-1) + xr.*hr(i,j) - xi.*hi(i,j);
hi(i,j-1) = hi(i,j-1) + xr.*hi(i,j) + xi.*hr(i,j);
end; i = fix(j+1);
% .......... ACCUMULATE TRANSFORMATIONS ..........
for i = low : igh;
%
zr(i,j-1) = zr(i,j-1) + xr.*zr(i,j) - xi.*zi(i,j);
zi(i,j-1) = zi(i,j-1) + xr.*zi(i,j) + xi.*zr(i,j);
%
end; i = fix(igh+1);
end; j = fix(en+1);
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);
%
hr(1,1) = norm;
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;
int_shape=zeros(int_shape);int_shape(:)=intmlv(1:numel(int_shape));intmlv=int_shape;
end %subroutine comlr2
%DECK COMLR
|
|
