/* ACM Algorithm 343. EIGENVALUES AND EIGENVECTORS OF A REAL GENERAL MATRIX */
/* J. Grad and M. A. Brebner, Computer Services, University of Birmingham, */
/* Birmingham 15, England. */
/* Translated to PL/I by R. A. Vowels, 18-24 February 2006. */
/* Gives the same results for the three matrixes given in the article
in the Communications of the ACM in which this algorithm was
first reported (within the limits of single precision used).
It also gave the same results for the 8 x 8 matrix in Knoble's
certification. */
/* This version runs under Digital Research PL/I. */
/*
This subroutine finds all the eigenvalues and the
eigenvectors of a real general matrix A of order N.
First in the subroutine SCALE the matrix is scaled so that
the corresponding rows and columns are approximately
balanced and then the matrix A is normalised so that the
value of the Euclidian norm of the matrix is equal to 1.
The eigenvalues are computed by the double-step method
in the subroutine HESQR.
The eigenvectors are computed by inverse iteration in
the subroutine REAL_VECTOR, for the real eigenvalues,
or in the subroutine COMPPLEX_VECTOR, for the complex eigenvalues.
The real parts of the N computed eigenvalues will be found
in the array EVR and the imaginary parts in the array EVI.
The real components of the normalised eigenvector I
(I = 1, 2, ... , N) corresponding to the eigenvalue stored in
EVR(I) and EVI(I) will be found in the column I of the
two-dimensional array VECR and the imaginary components
in the column I of the two-dimensional array VECI.
The real eigenvector is normalised so that the sum of the
squares of the two components is equal to one.
The complex eigenvector is normalised so that the
component with the largest value in modulus has its real
part equal to one and the imaginary part equal to zero.
The array INDIC indicates the success of the subroutine
EIGENP as follows:
Value of INDIC(I) Eigenvalue I Eigenvector I
0 Not found Not found
1 Found Not found
2 Found Found
*/
eigenp:
procedure (n, a, t, evr, evi, vecr, veci, indic);
declare n fixed binary;
declare (a(8,8), t, evr(20), evi(20), vecr(8,8), veci(8,8) ) float;
declare indic(20) fixed binary;
declare (iwork(20), local(100) ) fixed binary;
declare (prfact(20), subdia(20), work1(20), work2(20), work(20)) float;
declare (enorm, eps, ex, r, r1) float;
declare (d1, d2, d3) float,
(i, ivec, j, k, k1, kon, L, L1, m) fixed binary;
if n = 1 then
do;
evr(1) = a(1,1); evi(1) = 0; vecr(1,1) = 1;
veci(1,1) = 0; indic(1) = 2; return;
end;
call scale (n, a, veci, prfact, enorm);
/* The computation of the eigenvalues of the normalised matrix. */
ex = 1; r = 1;
loop: ex = ex/2;
if (r + ex) ^= r then go to Loop;
ex = exp(-t*log(2.0000000000000000e0));
call hesqr (n, a, veci, evr, evi, subdia, indic, eps, ex);
/* The possible decomposition of the upper-Hessenberg matrix
into the sub-matrices of lower order is indicated in the
array LOCAL. The decomposition occurs when some
sub-diagonal elements are in modulus less than a small
positive number EPS defined in the subroutine HESQR. The
amount of work in the eigenvector problem may be
diminished in this way.
*/
j = n; i = 1; local(1) = 1;
if j = 1 then go to L4;
L2: if abs(subdia(j-1)) > eps then go to L3;
i = i + 1;
local(i) = 0;
L3: j = j - 1;
local(i) = local(i) + 1;
if j ^= 1 then go to L2;
/* The eigenvector problem. */
L4: k = 1; kon = 0; L = local(1); m = n;
do i = 1 to n;
ivec = n - i + 1;
if i <= L then go to L5;
k = k + 1; m = n - L; L = L + local(k);
L5: if indic(ivec) = 0 then go to L10;
if evi(ivec) ^= 0 then go to L8;
/* Transfer of an upper-Hessenberg matrix of the order M from
the arrays VECI and SUBDIA into the array A. */
do k1 = 1 to m;
do L1 = k1 to m;
a(k1,L1) = veci(k1,L1);
end;
if k1 ^= 1 then
a(k1,k1-1) = subdia(k1-1);
end;
/* The computation of the real eigenvector IVEC of the upper-
Hessenberg matrix corresponding to the real eigenvalue
EVR(IVEC). */
call REAL_VECTOR (n, m, ivec, a, vecr, evr, evi, iwork, work, indic, eps, ex);
go to L10;
/* The computation of the complex eigenvector IVEC of the
upper-Hessenberg matrix corresponding to the complex
eigenvalue EVR(IVEC) + I*EVI(IVEC). If the value of KON is
not equal to zero then this complex eigenvector has
already been found from its conjugate. */
L8: if kon ^= 0 then go to L9;
kon = 1;
;
call COMPLEX_VECTOR (n, m, ivec, a, vecr, veci, evr, evi, indic,
iwork, subdia, work1, work2, work, eps, ex);
go to L10;
L9: kon = 0;
L10: end;
/* The reconstruction of the matrix used in the reduction of
the matrix A to an upper-Hessenberg form by Householder method. */
do i = 1 to n;
do j = i to n;
a(i,j) = 0;
a(j,i) = 0;
end;
a(i,i) = 1;
end;
if n <= 2 then go to L15;
m = n - 2;
do k = 1 to m;
L = k + 1;
do j = 2 to n;
d1 = 0;
do i = L to n;
d2 = veci(i,k);
d1 = d1 + d2*a(j,i);
end;
do i = L to n;
a(j,i) = a(j,i) - veci(i,k)*d1;
end;
end;
end;
/* The computation of the eigenvectors of the original
non-scaled matrix. */
L15: kon = 1;
do i = 1 to n;
L = 0;
if evi(i) = 0 then go to L16;
L = 1;
if kon = 0 then go to L16;
kon = 0;
go to L24;
L16: do j = 1 to n;
d1 = 0; d2 = 0;
do k = 1 to n;
d3 = a(j,k); d1 = d1 + d3*vecr(k,i);
if L ^= 0 then
d2 = d2 + d3*vecr(k,i-1);
end;
work(j) = d1/prfact(j);
if L ^= 0 then
subdia(j) = d2/prfact(j);
end;
/* The normalisation of the eigenvectors and the computation
of the eigenvalues of the original non-normalized matrix. */
if L = 1 then go to L21;
d1 = 0;
do m = 1 to n;
d1 = d1 + work(m)*work(m);
end;
d1 = sqrt(d1);
do m = 1 to n;
veci(m,i) = 0;
vecr(m,i) = work(m)/d1;
end;
evr(i) = evr(i)*enorm;
go to L24;
L21: kon = 1;
evr(i) = evr(i)*enorm;
evr(i-1) = evr(i);
evi(i) = evi(i)*enorm;
evi(i-1) = -evi(i);
r = 0;
do j = 1 to n;
r1 = work(j)*work(j) + subdia(j)*subdia(j);
if r < r1 then
do; r = r1; L = j; end;
end;
d3 = work(L);
r1 = subdia(L);
do j = 1 to n;
d1 = work(j); d2 = subdia(j);
vecr(j,i) = (d1*d3 + d2*r1)/r;
veci(j,i) = (d2*d3 - d1*r1)/r;
vecr(j,i-1) = vecr(j,i);
veci(j,i-1) = -veci(j,i);
end;
L24: end;
/* This subroutine stores the matrix of the order N from the
array A into the array H. Afterward the matrix in the
array A is scaled so that the quotient of the absolute sum
of the off-diagonal elements of column I and the absolute
sum of the off-diagonal elements of row I lies within the
values of BOUND1 and BOUND2.
The component I of the eigenvector obtained by using the
scaled matrix must be divided by the value found in the
PRFACT(I) of the array PRFACT. In this way the eigenvector
of the non-scaled matrix is obtained.
After the matrix is scaled it is normalised so that the
value of the Euclidean norm is equal to one.
If the process of scaling was not successful the original
matrix from the array H would be stored back into A and
the eigenproblem would be solved by using this matrix.
The eigenvalues of the normalised matrix must be
multiplied by the scalar ENORM in order that they become
the eigenvalues of the non-normalised matrix.
*/
scale: procedure (n, a, h, prfact, enorm);
declare n fixed binary,
(a(8,8), h(8,8), prfact(20), enorm) float;
declare (i, j, ncount, iter) fixed binary;
declare (column, row, factor, fnorm, q) float;
declare (bound1, bound2) float;
do i = 1 to n;
do j = 1 to n;
h(i,j) = a(i,j);
end;
prfact(i) = 1;
end;
bound1 = 0.75; bound2 = 1.33;
iter = 0;
L3: ncount = 0;
do i = 1 to n;
column = 0; row = 0;
do j = 1 to n;
if i ^= j then
do;
column = column + abs(a(j,i));
row = row + abs(a(i,j));
end;
end;
if (column ^= 0) & (row ^= 0) then
do;
q = column / row;
if q < bound1 then go to L6;
if q > bound2 then go to L6;
end;
ncount = ncount + 1;
go to L8;
L6: factor = sqrt(q);
do j = 1 to n;
if i ^= j then
do;
a(i,j) = a(i,j) * factor;
a(j,i) = a(j,i) / factor;
end;
end;
prfact(i) = prfact(i)*factor;
L8: end;
iter = iter + 1;
if iter > 30 then go to L11;
if ncount < n then go to L3;
fnorm = 0;
do i = 1 to n;
do j = 1 to n;
q = a(i,j);
fnorm = fnorm + q*q;
end;
end;
fnorm = sqrt(fnorm);
do i = 1 to n;
do j = 1 to n;
a(i,j) = a(i,j)/fnorm;
end;
end;
enorm = fnorm;
return;
/* Scaling was unsuccessful. Replace the original matrix A. */
L11: do i = 1 to n;
prfact(i) = 1;
do j = 1 to n;
a(i,j) = h(i,j);
end;
end;
enorm = 1;
end scale;
/* This subroutine finds all the eigenvalues of a real
general matrix. The original matrix A of order N is
reduced to the upper-Hessenberg form H by means of
similarity transformations (Householder method). The matrix
H is preserved in the upper half of the array H and in the
array SUBDIA. The special vectors used in the definition
of the Householder transformation matrices are stored in
the lower part of the array H.
The real parts of the N eigenvalues will be found in the
first N places of the array EVR, and
the imaginary parts in the first N places of the array EVI.
The array INDIC indicates the success of the routine as
follows:
Value of INDIC(I) Eigenvalue I
0 Not found
1 Found
EPS is a small positive number that numerically represents
zero in the program. EPS = (Euclidean norm of H)*EX , where
EX = 2**(-T). T is the number of binary digits in the
mantissa of a floating-point number.
*/
hesqr:
procedure (n, a, h, evr, evi, subdia, indic, eps, ex);
declare (n, indic(20)) fixed binary,
(a(8,8), h(8,8), evr(20), evi(20), subdia(20), eps, ex) float;
declare (r, shift, t, s, sr, sr2, x, y, z) float;
declare (i, j, k, L, maxst, m, m1, ns) fixed binary;
/*
Reduction of the matrix A to an upper-Hessenberg form H.
There are N-2 steps.
*/
if n < 2 then go to L14;
if n > 2 then go to L2;
subdia(1) = a(2,1);
go to L14;
L2: m = n - 2;
do k = 1 to m;
L = k + 1;
s = 0;
do i = L to n;
h(i,k) = a(i,k);
s = s + abs(a(i,k));
end;
if s = abs(a(k+1,k)) then
do;
subdia(k) = a(k+1,k);
h(k+1,k) = 0;
go to L12;
end;
sr2 = 0;
do i = L to n;
sr = a(i,k);
sr = sr / s;
a(i,k) = sr;
sr2 = sr2 + sr * sr;
end;
sr = sqrt(sr2);
if a(L,k) >= 0 then sr = -sr;
sr2 = sr2 - sr*a(L,k);
a(L,k) = a(L,k) - sr;
h(L,k) = h(L,k) - sr*s;
subdia(k) = sr*s;
x = s*sqrt(sr2);
do i = L to n;
h(i,k) = h(i,k) / x;
subdia(i) = a(i,k) / sr2;
end;
/* Pre-multiplication by the matrix PR. */
do j = L to n;
sr = 0;
do i = L to n;
sr = sr + a(i,k) * a(i,j);
end;
do i = L to n;
a(i,j) = a(i,j) - subdia(i)*sr;
end;
end;
/* Post-multiplication by the matrix PR. */
do j = 1 to n;
sr = 0;
do i = L to n;
sr = sr + a(j,i) * a(i,k);
end;
do i = L to n;
a(j,i) = a(j,i) - subdia(i) * sr;
end;
end;
L12: end;
do k = 1 to m;
a(k+1,k) = subdia(k);
end;
/* Transfer of the upper half of the matrix A into the
array H and the calculation of the small positive number eps. */
subdia(n-1) = a(n,n-1);
L14: eps = 0;
do k = 1 to n;
indic(k) = 0;
if k ^= n then eps = eps + subdia(k) * subdia(k);
do i = k to n;
h(k,i) = a(k,i);
eps = eps + a(k,i)*a(k,i);
end;
end;
eps = ex*sqrt(eps);
/* The QR iterative process. The upper-Hessenberg matrix H is
reduced to the upper-modified triangular form. */
/* Determination of the shift of origin for the first step of
the QR interative process. */
shift = a(n,n-1);
if n <= 2 then shift = 0;
if a(n,n) ^= 0 then shift = 0;
if a(n-1,n) ^= 0 then shift = 0;
if a(n-1,n-1) ^= 0 then shift = 0;
m = n; ns = 0; maxst = n*10;
/* Testing if the upper half of the matrix is equal to zero.
If it is equal to zero, the QR process is not necessary. */
do i = 2 to n;
do k = i to n;
if a(i-1,k) ^= 0 then go to L18;
end;
end;
do i = 1 to n;
indic(i) = 1;
evr(i) = a(i,i);
evi(i) = 0;
end;
return;
/* Start the main loop of the QR process. */
L18: k = m-1;
m1 = k;
i = k;
/* Find any decomposition of the matrix. */
/* Jump to L34 if the last sub-matrix of the decomposition is
of the order one. */
/* Jump to L35 if the last sub-matrix of the decomposition is
of the order two. */
if k < 0 then return;
if k = 0 then go to L34;
if abs(a(m,k)) <= eps then go to L34;
if m = 2 then go to L35;
L20: i = i - 1;
if abs(a(k,i)) <= eps then go to L21;
k = i;
if k > 1 then go to L20;
L21: if k = m1 then go to L35;
/* Transformation of the matrix of the order greater than two. */
s = a(m,m) + a(m1,m1) + shift;
sr = a(m,m) * a(m1,m1) - a(m,m1)*a(m1,m) + shift*shift/4;
if k < 1 then do; put skip list ('HESQR: k out of range'); stop; end;
a(k+2,k) = 0;
/* Calculate x1, y1, z1 for the sub-matrix obtained by the decomposition. */
x = a(k,k)*(a(k,k)-s) + a(k,k+1)*a(k+1,k) + sr;
y = a(k+1,k)*(a(k,k)+a(k+1,k+1)-s);
r = abs(x) + abs(y);
if r = 0 then do; shift = a(m,m-1); go to L21; end;
z = a(k+2,k+1) * a(k+1,k);
shift = 0;
ns = ns + 1;
/* The loop for one step of the QR process. */
do i = k to m1;
/* Calculate XR, YR, ZR. */
if i ^= k then
do;
x = a(i,i-1); y = a(i+1,i-1);
z = 0;
if i+2 <= m then z = a(i+2,i-1);
end;
sr2 = abs(x) + abs(y) + abs(z);
if sr2 ^= 0 then
do; x = x/sr2; y = y/sr2; z = z/sr2; end;
s = sqrt(x*x + y*y + z*z);
if x >= 0 then s = -s;
if i ^= k then a(i,i-1) = s*sr2;
if sr2 ^= 0 then go to L26;
if i+3 > m then go to L33;
go to L32;
L26: sr = 1 - x/s;
s = x - s;
x = y/s;
y = z/s;
/* Pre-multiplication by the matrix PR. */
do j = i to m;
s = a(i,j) + a(i+1,j)*x;
if i+2 <= m then s = s + a(i+2,j)*y;
s = s * sr;
a(i,j) = a(i,j) - s;
a(i+1,j) = a(i+1,j) - s*x;
if i+2 <= m then a(i+2,j) = a(i+2,j) - s*y;
end;
/* Post-multiplication by the matrix PR. */
L = i + 2;
if i >= m1 then L = m;
do j = k to L;
s = a(j,i) + a(j,i+1)*x;
if i+2 <= m then s = s + a(j,i+2)*y;
s = s * sr;
a(j,i) = a(j,i) - s;
a(j,i+1) = a(j,i+1) - s*x;
if i+2 <= m then a(j,i+2) = a(j,i+2) - s*y;
end;
if i+3 > m then go to L33;
s = -A(i+3,i+2)*y*sr;
L32: a(i+3,i) = s;
a(i+3,i+1) = s*x;
a(i+3,i+2) = s*y + a(i+3,i+2);
L33:
end; /* of loop for one step of the QR process. */
if ns > maxst then return;
go to L18;
/* Compute the last eigenvalue. */
L34: evr(m) = a(m,m);
evi(m) = 0;
indic(m) = 1;
m = k;
go to L18;
/* Compute the eigenvalues of the last 2 x 2 matrix obtained by
the decomposition. */
L35: r = (a(k,k) + a(m,m))/2;
s = (a(m,m) - a(k,k))/2;
s = s*s + a(k,m)*a(m,k);
indic(k) = 1;
indic(m) = 1;
if s >= 0 then
do;
t = sqrt(s);
evr(k) = r-t; evi(k) = 0;
evr(m) = r+t; evi(m) = 0;
m = m - 2;
go to L18;
end;
L36: t = sqrt(-s);
evr(k) = r; evi(k) = t;
evr(m) = r; evi(m) = -t;
m = m - 2;
go to L18;
end hesqr;
/* This subroutine finds the real eigenvector of the real
upper-Hessenberg matrix in the array A, corresponding to
the real eigenvalue stored in EVR(IVEC). The inverse
iteration method is used.
Note the matrix in A is destroyed by the subroutine.
N is the order of the upper-Hessenberg matrix.
M is the order of the sub-matrix obtained by a suitable
decomposition of the upper-Hessenberg matrix if some
sub-diagonal elements are equal to zero. The value of M is
chosen so that the last N-M components of the eigenvector
are zero.
IVEC gives the position of the eigenvalue in the array EVR
for which the corresponding eigenvector is computed.
The array EVI would contain the imaginary parts of the N
eigenvalues if they existed.
The M components of the computed real eigenvector will be
found in the first N places of the column IVEC of the two-
dimensional array VECR.
IWORK and WORK are the working stores used during the
Gaussian elimination and backward substitution process.
The array INDIC indicates the success of the routine as
follows:
Value of INDIC(I) Eigenvector I
1 not found
2 found
EPS is a small positive number that numerically represents
zero in the program. EPS = (Euclidian norm of A)*EX, where
EX = 2**(-T). T is the number of binary digits in the
mantissa of a floating-point number.
*/
REAL_VECTOR:
procedure (n, m, ivec, a, vecr, evr, evi, iwork, work, indic, eps, ex);
declare (n, m, ivec, iwork(20), indic(20)) fixed binary,
(a(8,8), vecr(8,8), evr(20), evi(20), work(20), eps, ex) float;
declare (s, sr) float;
declare (bound, evalue, previs, r, r1, t) float;
declare (i, j, k, L, iter, ns) fixed binary;
vecr(1,ivec) = 1;
if m = 1 then go to L24;
/* Small perturbation of equal eigenvectors to obtain a full
set of eigenvectors. */
evalue = evr(ivec);
if ivec = m then go to L2;
k = ivec + 1;
r = 0;
do i = k to m;
if evalue ^= evr(i) then go to L1;
if evi(i) = 0 then r = r + 3;
L1: end;
evalue = evalue + r*ex;
L2: do k = 1 to m;
a(k,k) = a(k,k) - evalue;
end;
/* Gaussian elimination of the upper-Hessenberg matrix A. All
row interchanges are indicated in the array IWORK. All the
multipliers are stored as the sub-diagonal elements of A. */
k = m - 1;
do i = 1 to k;
L = i + 1;
iwork(i) = 0;
if a(i+1,i) ^= 0 then go to L4;
if a(i,i) ^= 0 then go to L8;
a(i,i) = eps;
go to L8;
L4: if abs(a(i,i)) >= abs(a(i+1,i)) then go to L6;
iwork(i) = 1;
do j = i to m;
r = a(i,j);
a(i,j) = a(i+1,j);
a(i+1,j) = r;
end;
L6: r = -a(i+1,i) / a(i,i);
a(i+1,i) = r;
do j = L to m;
a(i+1,j) = a(i+1,j) + r*a(i,j);
end;
L8: end;
if a(m,m) = 0 then a(m,m) = eps;
/* The vector (1, 1, ... , 1) is stored in the place of the right-
hand side column vector. */
do i = 1 to n;
if i <= m then
work(i) = 1;
else
work(i) = 0;
end;
/* The inverse iteration is performed on the matrix until the
infinite norm of the right-hand side vector is greater
than the bound defined as 0.01/(n*ex). */
bound = 1.00000000000000000e-2 / (ex*n);
ns = 0;
iter = 1;
/* The backward substitution. */
L12: r = 0;
do j = m to 1 by -1;
s = work(j);
if j ^= m then
do;
do k = j+1 to m;
sr = work(k);
s = s - sr*a(j,k);
end;
end;
work(j) = s / a(j,j);
t = abs(work(j));
if r < t then r = t;
end;
/* The computation of the right-hand side vector for the new
iteration step. */
do i = 1 to m;
work(i) = work(i) / r;
end;
/* The computation of the residuals and comparison of the
residuals of the two successive steps of the inverse
iteration. If the infinite norm of the residual vector is
greater than the infinite norm of the previous residual
vector, the computed eigenvector of the previous step is
taken as the final eigenvector. */
r1 = 0;
do i = 1 to m;
t = 0;
do j = i to m;
t = t + a(i,j) * work(j);
end;
t = abs(t);
if r1 < t then r1 = t;
end;
if iter = 1 then go to L19;
if previs <= r1 then go to L24;
L19: do i = 1 to m;
vecr(i,ivec) = work(i);
end;
previs = r1;
if ns = 1 then go to L24;
if iter > 6 then go to L25;
iter = iter + 1;
if r >= bound then ns = 1;
/* Gaussian elimination of the right-hand side vector. */
do i = 1 to m-1;
r = work(i+1);
if iwork(i) ^= 0 then
do;
work(i+1) = work(i) + work(i+1)*a(i+1,i);
work(i) = r;
end;
else
work(i+1) = work(i+1) + work(i)*a(i+1,i);
end;
go to L12;
L24: indic(ivec) = 2;
L25: do i = m+1 to n;
vecr(i,ivec) = 0;
end;
end REAL_VECTOR;
/* This subroutine finds the complex eigenvector of the real
upper-Hessenberg matrix of order N corresponding to the
complex eigenvalue with the real part in EVR(IVEC) and the
corresponding imaginary part in EVI(IVEC). The inverse
iteration method is used modified to avoid the use of
complex arithmetic.
The matrix on which the inverse iteration is performed is
built up in the array A by using the upper-Hessenberg
matrix preserved in the upper half of the array H and in
the array SUBDIA.
M is the order of the sub-matrix obtained by a suitable
decomposition of the upper Hessenberg matrix if some
sub-diagonal elements are equal to zero. The value of M is
chosen so that the last N-M components of the complex
eigenvector are zero.
The real parts of the first M components of the computed
complex eigenvector will be found in the first M places of
the column whose top element is VECR(1,IVEC) and the
corresponding imaginary parts of the first M components of
the complex eigenvector will be found in the first M
places of the column whose top element is VECR(1,IVEC-1).
The array INDIC indicates the success of the routine as
follows:
Value of INDIC(I) Eigenvalue I
1 Not found
2 Found
The arrays IWORK, WORK1, WORK2, and WORK are the working
stores used during the inverse iteration process.
EPS is a small positive number that numerically represents
zero in the program. EPS = (Euclidean norm of H)*EX, where
EX = 2**(-T). T is the number of binary digits in the
mantissa of a floating-point number.
*/
COMPLEX_VECTOR:
procedure (n, m, ivec, a, vecr, h, evr, evi, indic, iwork, subdia,
work1, work2, work, eps, ex);
declare (n, m, ivec, indic(20), iwork(20)) fixed binary,
(a(8,8), vecr(8,8), h(8,8), evr(20), evi(20), subdia(20),
work1(20), work2(20), work(20), eps, ex) float;
declare (d, d1) float;
declare (b, bound, eta, fksi, previs, r, s, u, v) float;
declare (i, i1, i2, iter, j, k, L, ns) fixed binary;
fksi = evr(ivec);
eta = evi(ivec);
/* The modification of the eigenvalue (fksi + i*eta) if more
eigenvalues are equal. */
if ivec ^= m then
do;
r = 0;
do i = ivec+1 to m;
if fksi ^= evr(i) then go to L1;
if abs(eta) ^= abs(evi(i)) then go to L1;
r = r + 3;
L1: end;
r = r * ex;
fksi = fksi + r;
eta = eta + r;
end;
/* The matrix ((H-FKSI*I)*(H-FKSI*I) + (ETA*ETA)*I) is
stored into the array A. */
r = fksi*fksi + eta*eta;
s = 2 * fksi;
L = m - 1; /* This appears to be able to be moved to next loop */
/* & combined with DO parameter. */
do i = 1 to m;
do j = i to m;
d = 0;
a(j,i) = 0;
do k = i to j;
d = d + h(i,k) * h(k,j);
end;
a(i,j) = d - s*h(i,j);
end;
a(i,i) = a(i,i) + r;
end;
do i = 1 to L;
r = subdia(i);
a(i+1,i) = -s*r;
do j = 1 to i+1;
a(j,i) = a(j,i) + r*h(j,i+1);
end;
if i ^= 1 then a(i+1,i-1) = r*subdia(i-1);
do j = i to m;
a(i+1,j) = a(i+1,j) + r*h(i,j);
end;
end;
/* The Gaussian elimination of the matrix
((H-FKSI*I)*(H-FKSI*I) + (ETA*ETA)*I) in the array A. The
row interchanges that occur are indicated in the array
IWORK. All the multipliers are stored in the first and in
the second sub-diagonal of the array A. */
k = m - 1;
do i = 1 to k;
i1 = i + 1;
i2 = i + 2;
iwork(i) = 0;
if i = k then go to L10;
if a(i+2,i) ^= 0 then go to L11;
L10: if a(i+1,i) ^= 0 then go to L11;
if a(i,i) ^= 0 then go to L18;
a(i,i) = eps;
go to L18;
L11: if i = k then go to L12;
if abs(a(i+1,i)) >= abs(a(i+2,i)) then go to L12;
if abs(a(i,i)) >= abs(a(i+2,i)) then go to L16;
L = i + 2;
iwork(i) = 2;
go to L13;
L12: if abs(a(i,i)) >= abs(a(i+1,i)) then go to L15;
L = i + 1;
iwork(i) = 1;
L13: do j = i to m;
r = a(i,j);
a(i,j) = a(L,j);
a(L,j) = r;
end;
L15: if i = k then i2 = i1;
L16: do L = i1 to i2;
r = -a(L,i) / a(i,i);
a(L,i) = r;
do j = i1 to m;
a(L,j) = a(L,j) + r*a(i,j);
end;
end;
L18: end;
if a(m,m) = 0 then a(m,m) = eps;
/* The vector (1, 1, ... , 1) is stored into the right-hand side
vectors vecr( ,ivec) and vecr( ,ivec-1) representing the
complex right-hand side vector. */
do i = 1 to n;
if i <= m then
do; vecr(i,ivec) = 1; vecr(i,ivec-1) = 1; end;
else
do; vecr(i,ivec) = 0; vecr(i,ivec-1) = 0; end;
end;
/* The inverse iteration is performed on the matrix until the
infinite norm of the right-hand side vector is greater
than the bound defined as 0.01 / (n*ex). */
bound = 1.00000000000000000e-2 / (ex*n);
ns = 0;
iter = 1;
do i = 1 to m;
work(i) = h(i,i) - fksi;
end;
/* The sequence of the complex vectors Z(S) = P(S) + I*Q(S) and
W(S+1) = U(S+1) + I*V(S+1) is given by the relations
(A - (FKST - I*ETA)*I)*W(S+1) = Z(S) and
Z(S+1) = W(S+1)/MAX(W(S+1)).
The final W(S) is taken as the computed eigenvector. */
/* The computation of the right-hand side vector
(A-FKST*I)*PS-ETA*Q(S). A is an upper-Hessenberg matrix. */
L23: do i = 1 to m;
d = work(i) * vecr(i,ivec);
if i ^= 1 then d = d + subdia(i-1) * vecr(i-1,ivec);
do k = i+1 to m;
d = d + h(i,k) * vecr(k, ivec);
end;
vecr(i,ivec-1) = d - eta*vecr(i,ivec-1);
end;
/* Gaussian elimination of the right-hand side vector. */
k = m-1;
do i = 1 to k;
L = i + iwork(i);
r = vecr(L,ivec-1);
vecr(L, ivec-1) = vecr(i,ivec-1);
vecr(i,ivec-1) = r;
vecr(i+1,ivec-1) = vecr(i+1,ivec-1) + a(i+1,i)*r;
if i ^= k then
vecr(i+2,ivec-1) = vecr(i+2,ivec-1) + a(i+2,i)*r;
end;
/* The computation of the real part U(S+1) of the complex */
/* vector W(S+1). The vector U(S+1) is obtained after the */
/* backward substitution. */
do i = 1 to m;
j = m-i+1;
d = vecr(j,ivec-1);
if j ^= m then
do;
do k = j+1 to m;
d1 = a(j,k);
d = d - d1*vecr(k,ivec-1);
end;
end;
vecr(j,ivec-1) = d / a(j,j);
end;
/* The computation of the imaginary part V(S+1) of the vector
W(S+1), where V(S+1) = (P(S) - (A-FKSI*I)*U(S+1)) / ETA. */
do i = 1 to m;
d = work(i) * vecr(i,ivec-1);
if i ^= 1 then d = d + subdia(i-1) * vecr(i-1,ivec-1);
do k = i+1 to m;
d = d + h(i,k) * vecr(k,ivec-1);
end;
vecr(i,ivec) = (vecr(i,ivec) - d) / eta;
end;
/* The computation of (infinite norm of W(S+1))**2 */
L = 1; s = 0;
do i = 1 to m;
r = vecr(i,ivec)*vecr(i,ivec) + vecr(i,ivec-1)*vecr(i,ivec-1);
if r > s then do; s = r; L = i; end;
end;
/* The computation of the vector Z(S+1), where
Z(S+1) = W(S+1) / (component of W(S+1) with the largest absolute
value). */
u = vecr(L,ivec-1);
v = vecr(L,ivec);
do i = 1 to m;
b = vecr(i,ivec);
r = vecr(i,ivec-1);
vecr(i,ivec) = (r*u + b*v) / s;
vecr(i,ivec-1) = (b*u - r*v) / s;
end;
/* The computation of the residuals and comparison of the
residuals of the two successive steps of the inverse
iteration. If the infinite norm of the residual vector is
greater than the infinite norm of the previous residual
vector, the computed vector of the previous step is taken
as the computed approximation to the eigenvector.
*/
b = 0;
do i = 1 to m;
r = work(i)*vecr(i,ivec-1) - eta*vecr(i,ivec);
u = work(i)*vecr(i,ivec) + eta*vecr(i,ivec-1);
if i ^= 1 then
do;
r = r + subdia(i-1) * vecr(i-1,ivec-1);
u = u + subdia(i-1) * vecr(i-1,ivec);
end;
do j = i+1 to m;
r = r + h(i,j) * vecr(j,ivec-1);
u = u + h(i,j) * vecr(j,ivec);
end;
u = r*r + u*u;
if b < u then b = u;
end;
if iter = 1 then go to L42;
if previs <= b then go to L44;
L42: do i = 1 to n;
work1(i) = vecr(i,ivec);
work2(i) = vecr(i,ivec-1);
end;
previs = b;
if ns = 1 then go to L46;
if iter > 6 then go to L47;
iter = iter + 1;
if bound > sqrt(s) then go to L23;
ns = 1;
go to L23;
L44: do i = 1 to n;
vecr(i,ivec) = work1(i);
vecr(i,ivec-1) = work2(i);
end;
L46: indic(ivec-1) = 2;
indic(ivec) = 2;
L47:
end COMPLEX_VECTOR;
end eigenp;