Page 126 - Compact Numerical Methods For Computers

The algebraic eigenvalue problem 115
Algorithm 26. Eigensolutions of a complex matrix by Eberlein’s method (cont.)
begin
hr := hr+A[k,i]*A[m,i]+Z[k,i]*Z[m,i];
hr := hr-A[i,k]*A[i,m]-Z[i,k]*Z[i,m];
hi := hi+Z[k,i]*A[m,i]-A[k,i]*Z[m,i];
hi := hi-A[i,k]*Z[i,m]+Z[i,k]*A[i,m];
te := A[i,k]*A[i,k]+Z[i,k]*Z[i,k]+A[m,i]*A[m,i]+Z[m,i]*Z[m,i];
tee := A[i,m]*A[i,m]+Z[i,m]*Z[i,m]+A[k,i]*A[k,i]+Z[k,i]*Z[k,i];
g := g+te+tee; hj := hj-te+tee;
end; {if i<>k and i<>m}
end; {loop on i}
br := A[k,m]+A[m,k]; bi := Z[k,m]+Z[m,k]; er := A[k,m]-A[m,k];
ei := Z[k,m]-Z[m,k]; dr := A[k,k]-A[m,m]; di := Z[k,k]-Z[m,m];
te := br*br+ei*ei+dr*dr; tee := bi*bi+er*er+di*di;
if te>=tee then
begin
isw := 1.0; c := br; s := ei; d := dr, de := di;
root2 := sqrt(te);
e n d
else {te<tee}
begin
isw := -1.0; c := bi; s := -er; d := di; de := dr,
root2 := sqrt(tee);
end;
root1 := sqrt(s*s+c*c); sig := -1.0; if d>=0.0 then sig := 1.0;
sa := 0.0; ca := -1.0; if c>=0.0 then ca := 1.0;
if root1<=eps then
begin
sx := 0.0; sa := 0.0; cx := 1.0; ca := 1.0;
if isw<=0.0 then
begin
e := ei; bv := -br;
end
e l s e
begin
e := er; bv := bi;
end; {if isw<=0.0}
nd := d*d+de*de;
end
else {root1>eps}
begin
if abs(s)>eps then
begin
ca := c/root1; sa := s/root1;
end; {abs(s)>eps}
cot2x := d/rootl; cotx := cot2x+(sig*sqrt(1.0+cot2x*cot2x));
sx := sig/sqrt(1.0+cotx*cotx); cx := sx*cotx;
{find rotated elements}
eta := (er*br+ei*bi)/root1; tse := (br*bi-er*ei)/root1;
te := sig*(tse*d-de*root1)/root2; tee := (d*de+root1*tse)/root2;
nd := root2*root2+tee*tee; tee := hj*cx*sx; cos2a := ca*ca-sa*sa;
sin2a := 2.0*ca*sa; tem := hr*cos2a+hi*sin2a;
tep := hi*cos2a-hr*sin2a; hr := hr*cx*cx-tem*sx*sx-ca*tee;