Page 124 - Compact Numerical Methods For Computers

The algebraic eigenvalue problem 113

Algorithm 12. Residuals of a complex eigensolution (cont.)
I
end; {loop on j}
writeln(‘Sum of squares = ’,ss);
writeln;
end; {alg12.pas == comres}

Algorithm 26. Eigensolutions of a complex matrix by Eberlein’s method
procedure comeig( n : integer; {order of problem}
var itcount: integer; {on entry, a limit to the iteration
count; on output, the iteration count to convergence
which is set negative if the limit is exceeded}
var A, Z, T, U : rmatrix); {the matrix for which the
eigensolutions are to be found is A + sqrt(-1)*Z,
and this will be transformed so that the diagonal
elements become the eigenvalues; on exit,
T + sqrt(-1)*U has the eigenvectors in its columns}
{alg26.pas == Pascal version of Eberlein’s comeig.
Translated from comeig.bas, comeig.for and the original Algol
version in Eberlein (1971).
Copyright J C Nash 1988
}
var
Rvec : rvector;
{Rvec was called en by Eberlein, but we use upper case vector names.}
i, itlimit, j, k, kl, m, m9, n1 : integer;
aki, ami, bv, br, bi : real;
c, cli, clr, c2i, c2r, ca, cb, ch, cos2a, cot2x, cotx, cx : real;
d, de, di, diag, dr, e, ei, er, eps, eta, g, hi, hj, hr : real;
isw, max, nc, nd, root1, root2, root : real;
s, s1i, s1r, s2i, s2r, sa, sb, sh, sig, sin2a, sx : real;
tanh, tau, te, tee, tern, tep, tse, zki, zmi : real;
mark : boolean;
begin {comeig}
{Commentary in this version has been limited to notes on differences
between this implementation and the Algol original.}
writeln(‘alg26.pas -- comeig’);
eps := Calceps; {calculate machine precision}
{NOTE: we have not scaled eps here, but probably should do so to avoid
unnecessary computations.}
mark := false; n1 := n-1;
for i := l to n do
begin
for j := 1 to n do
begin {initialize eigenvector arrays}
T[i,j] := 0.0; U[i,j] := 0.0; if i=j then T[i,i] := 1.0;
end; {loop on j}
end, {loop on i}
itlimit := itcount; {use value on entry as a limit}
itcount := 0; {and then set counter to zero}
while (itcount<=itlimit) and (not mark) do