Page 127 - Compact Numerical Methods For Computers

116 Compact numerical methods for computers

Algorithm 26. Eigensolutions of a complex matrix by Eberlein’s method (cont.)
hi := hi*cx*cx+tep*sx*sx-sa*tee;
bv := isw*te*ca+eta*sa; e := ca*eta-isw*te*sa;
end; {root1>eps}
{label ‘enter1’ is here in Algol version}
s := hr-sig*root2*e; c := hi-sig*root2*bv; root := sqrt(c*c+s*s);
if root<eps then
begin
cb := 1.0; ch := 1.0; sb := 0.0; sh := 0.0;
end {if root<eps}
else {root>=eps}
begin
cb := -c/root; sb := s/root; tee := cb*bv-e*sb; nc := tee*tee;
tanh :=root/(g+2.0*(nc+nd)); ch := 1.0/sqrt(1.0-tanh*tanh);
sh := ch*tanh;
end; {root>=eps}
tem := sx*sh*(sa*cb-sb*ca); c1r := cx*ch-tem; c2r := cx*ch+tem;
c1i := -sx*sh*(ca*cb+sa*sb); c2i := c1i; tep := sx*ch*ca;
tem := cx*sh*sb* s1r := tep-tem; s2r := -tep-tem; tep := sx*ch*sa;
tem := cx*sh*cb: s1i := tep+tem; s2i := tep-tem;
tem := sqrt(s1r*s1r+s1i*s1i); tep := sqrt(s2r*s2r+s2i*s2i);
if tep>eps then mark := false;
if (tep>eps) and (tem>eps) then
begin
for i := 1 to n do
begin
aki := A[k,i]; ami := A[m,i]; zki := Z[k,i]; zmi := Z[m,i];
A[k,i] := c1r*aki-c1i*zki+s1r*ami-s1i*zmi;
Z[k,i] := c1r*zki+c1i*aki+s1r*zmi+s1i*ami;
A[m,i] := s2r*aki-s2i*zki+c2r*ami-c2i*zmi;
Z[m,i] := s2r*zki+s2i*aki+c2r*zmi+c2i*ami;
end; {loop on i}
for i := l to n do
begin
aki := A[i,k]; ami := A[i,m]; zki := Z[i,k]; zmi := Z[i,m];
A[i,k] := c2r*aki-c2i*zki-s2r*ami+s2i*zmi;
Z[i,k] := c2r*zki+c2i*aki-s2r*zmi-s2i*ami;
A[i,m] := -s1r*aki+s1i*zki+c1r*ami-c1i*zmi;
Z[i,m] := -s1r*zki-s1i*aki+c1r*zmi+c1i*ami;
aki := T[i,k]; ami := T[i,m]; zki := U[i,k]; zmi := U[i,m];
T[i,k] := c2r*aki-c2i*zki-s2r*ami+s2i*zmi;
U[i,k] := c2r*zki+c2i*aki-s2r*zmi-s2i*ami;
T[i,m] := -s1r*aki+s1i*zki+c1r*ami-c1i*zmi;
U[i,m] := -s1r*zki-s1i*aki+c1r*zmi+c1i*ami;
end; {loop on i}
end; {if tep and tern>eps}
end; {loop on m}
end; {loop on k}
end {if tau>=l00*eps}
else mark := true; {to indicate convergence}
end; {while itcount<=itlimit}
if itcount>itlimit then itcount := -itcount; {negative iteration count
means process has not converged properly}
end; {alg26.pas == comeig}