Page 117 - Compact Numerical Methods For Computers

106 Compact numerical methods for computers
quite satisfactorily. Indeed I have performed a number of successful computations
on a non-symmetric generalised eigenvalue problem wherein both A and B were
singular by means of this artifice! However, it must be admitted that this
particular problem arose in a context in which a great deal was known about the
properties of the solutions.

Algorithm 10. Inverse iteration via Gauss elimination
The algorithm requires a working array W, n by 2 * n and two vectors x and y of order n.
procedure gii(nRow : integer; {order of problem}
var A : rmatrix; {matrices of interest packed
into array as ( A | B) }
var Y : rvector; {the eigenvector}
var shift : real; {the shift -- on input, the value
nearest which an eigenvalue is wanted. On output, the
eigenvalue approximation found.}
var itcount: integer); {On input a limit to the number
of iterations allowed to find the eigensolution. On
output, the number of iterations used. The returned
value is negative if the limit is exceeded.}
{alg10.pas == Inverse iteration to find matrix an eigensolution of
A Y = e v * B * Y
for real matrices A and B of order n. The solution found corresponds to
one of the eigensolutions having eigenvalue, ev, closest to the value
shift. Y on input contains a starting vector approximation.
Copyright 1988 J.C.Nash
}
var
i, itlimit, j, k, m, msame, nRHS :integer;
ev, s, t, to1 : real; {eigenvalue approximation}
X : rvector;
begin
itlimit:=itcount; {to get the iteration limit from the call}
nRHS:=nRow; {the number of right hand sides is also n since we will
store matrix B in array A in the right hand nRow columns}
tol:=Calceps;
s:=0.0; {to initialize a matrix norm}
for i:=1 to nRow do
begin
X[i]:=Y[i]; {copy initial vector to X for starting iteration}
Y[i]:=0.0; {to initialize iteration vector Y}
for j:=1 to nRow do
begin
A[i,j]:=A[i,j]-shift*A[i,j+nRow];
s:=s+abs(A[i,j]);
end;
end;
tol:=tol*s; {to set a reasonable tolerance for zero pivots}
gelim(nRow, nRHS, A, tol); {Gauss elimination STEPS 2-10}
itcount:=0,
msame :=0; {msame counts the number of eigenvector elements which
are unchanged since the last iteration}