/************************************************************************/
/* */
/* ACM Algorithm 405 Roots of Matrix Pencils: */
/* The Generalized Eigenvalue Problem. */
/* */
/* Alice M. Dell, Roman L. Weiss, Gerald L. Thompson.*/
/* 25 May 1970 and 12 October 1970. */
/* */
/************************************************************************/
/* Translated from Algol 60 to PL/I by R. A. Vowels, 21 January 2008. */
/* Status: This procedure gives the same results as the original for */
/* Exercises 1, 2, and 11. */
/* Correct results are obtained for Ex. 3. */
/* Notes: */
/* 1. m and n have been omitted from the parameter list of PENCIL, */
/* because A and B are of size m x n, and the bounds are retrieved */
/* using HBOUND. */
/* 2. Lambda is now a COMPLEX vector. */
/* 3. The END statement of PENCIL was altered so as to encompass */
/* procedure EIG which referenced variable Sp (the number of */
/* eigenvalues found). */
/* 4. It was necessary to create a temporary variable XP in REDUCE for */
/* parameter X, as changes to X in REDUCE must not affect the */
/* corresponding argument. */
/* PENCIL finds the generalized eigenvalues LAMBDA which solve */
/* x(A-Lambda.B) = 0 and (A-Lambda.B)y = 0 and simultaneously reduce */
/* the rank of (A - Lambda.B), where A and B are m x n matrices. PENCIL */
/* converts the m by n problem of finding the rank-reducing numbers of */
/* (A - Lambda.B) into an ordinary r by r eigenvalue problem by a */
/* sequence of elementary transformations. */
/* Example (1).
10 2 3 1 1 12 1 -1 2 1
2 12 1 2 1 1 14 1 -1 1
A = 3 1 11 1 -1 B = -1 1 16 -1 1
1 2 1 9 1 2 -1 -1 12 -1
1 1 -1 1 15 1 1 1 -1 11
The derived matrix is:
.7220 .0248 .0871 .2358 -.1426
.0160 .8663 .1882 .0636 .0015
A5 = .0781 .2339 .7967 -.0357 .2016
.3107 .0554 -.0699 .8551 -.0791
-.1792 .0261 .1447 -.0173 1.4020
A5 is effectively B**-1.A */
/* Example (2);
2 3 1 2 -1 1
A = 3 5 2 B = 2 2 2
3 4 2 2 -1 1
A and B are transformed to the 1 x 1 matrices:
A' = [.23077] B = [1.1538]
for re-entry to PENCIL, which is called recursively.
The final output from PENCIL is A1 = [.2]
These values of A', B', and A1 are obtained by this
PL/I procedure.
Note: Since parameter A is call-by-value (in the ALGOL
version) and by a dummy parameter in the PL/I version,
there can be no output of A from PENCIL unless some
PUT statements are included in that procedure.
Note also that the final value of A1 = [0.2] is
the value obtained at the termination of the recursive call. */
/* Example (3).
2 3 2 1 0 0
A = 3 5 2 B = 0 1 0
2 2 2 0 0 0
Here m-r-q = n-r-q = 0 so that there were no recursive calls of
PENCIL. On exit from PENCIL,
A2 = 0 1
1 3 */
/* Example (11).
This is the same as Example (2), except that the values of A and B
are interchanged.
3 0 1
A3 = -7 2 -1
4 -2 0
The eigenvalue is the reciprocal of the eigenvalue from exercise (2).*/
(SUBRG, FOFL, SIZE):
PENCIL: PROCEDURE (AP, BP, Lambda, Sp, Par, Tol) RECURSIVE OPTIONS (REORDER);
/* AP = input matrix of size m x n */
/* BP = input matrix of size m x n */
/* Lambda = eigenvalues found. */
/* Sp = number of eigenvalues found. */
/* Par = 0 if there are no solutions. */
/* = 1 if there are solutions. */
/* Tol = Criteria for determining accuracy. */
DECLARE (AP(*,*), BP(*,*), Lambda(*) COMPLEX ) FLOAT (18);
DECLARE (Sp, Par) FIXED BINARY (31);
DECLARE Tol FLOAT;
DECLARE (m, n) FIXED BINARY (31);
DECLARE (A, B) (HBOUND(AP,1), HBOUND(AP,2)) FLOAT (18);
DECLARE P1 (HBOUND(A,1), HBOUND(A,1)) FLOAT (18),
P2 (HBOUND(A,2), HBOUND(A,2)) FLOAT (18);
DECLARE (i, j, k, limc, limr, q, r, s, t) FIXED BINARY (31);
m = HBOUND(A, 1); n = HBOUND(A,2);
A = AP; B = BP;
Par = 1; Sp = 0;
if m < n then k = m; else k = n;
do j = 1 to k;
Lambda(j) = 0;
end;
CALL REDUCE (p1, B, p2, m, n, r, Tol);
if r = 0 then
do;
par = 0; go to Endp;
end;
CALL MATMUL (p1, A, A, m, m, n); CALL MATMUL (A, p2, A, m, n, n);
CALL MATMUL (p1, B, B, m, m, n); CALL MATMUL (B, p2, B, m, n, n);
/* Note: The last two matrix multiplications together, by definition */
/* of P1 and P2, change B to */
/* | Ir 0 | */
/* | 0 0 | */
if n = r & m = r then
do;/* Calculations for examples (1) and (11) exit here, and for */
/* example (2) cease here after one recursive call. */
CALL EIG (Lambda, A, r);
go to Endp;
end;
if n = r then limc = 1; else limc = n-r;
if m = r then limr = 1; else limr = m-r;
begin;
DECLARE (C12(r,limc), C21(limr,r), C22(limr, limc),
P1(limr,limr), P2(limc,limc) ) FLOAT (18);
if n ^= r then
do i = 1 to r;
do j = r+1 to n;
C12(i,j-r) = A(i,j);
end;
end;
if m ^= r then
do i = r+1 to m;
do j = 1 to r; C21(i-r,j) = A(i,j); end;
if n ^= r then
do j = r+1 to n; C22(i-r,j-r) = A(i,j); end;
end;
/* A has now been partitioned and the parts are referred to below */
/* as */
/* A C12 */
/* C21 C22 */
if n = r | m = r then
do;
q = 0; go to mul;
end;
CALL REDUCE (P1, C22, P2, limr, limc, q, Tol);
CALL MATMUL (C12, P2, C12, r, limc, limc);
CALL MATMUL (P1, C21, C21, limr, limr, r);
CALL MATMUL (P1, C22, C22, limr, limr, limc);
CALL MATMUL (C22, P2, C22, limr, limc, limc);
/* C22 has been "reduced" and the requisite operations have been */
/* performed on C12 and C21. That is, */
/* C22 = P1 * C12 * P2 = Iq 0 */
/* 0 0 */
/* */
/* C12 = P1 * C12, and C21 = C21 * P2. */
if q = 0 then go to mul;
begin;
DECLARE D21(r,r) FLOAT (18);
CALL MATMUL (C12, c21, d21, r, q, r);
do i = 1 to r;
do j = 1 to r;
A(i,j) = A(i,j) - D21(i,j);
end;
end;
end;
Dstep:
if (m-r-q) = 0 & (n-r-q) = 0 then
do;/* Calculations for example (3) cease here. */
CALL EIG(Lambda, A, r); go to Endp;
end;
Mul:
if m-r-q = 0 then limr = 1; else limr = m-r-q;
if n-r-q = 0 then limc = 1; else limc = n-r-q;
begin;
DECLARE (E12(r,limc), E21(limr,r) ) FLOAT (18);
if n-r-q ^= 0 then
do i = 1 to r;
do j = q+1 to n-r;
E12(i,j-q) = C12(i,j);
end;
end;
if m-r-q ^= 0 then
do i = q+1 to m-r;
do j = 1 to r;
E21(i-q,j) = C21(i,j);
end;
end;
/* The columns of C12 above Iq and the rows of C21 to the left */
/* of Iq are annihilated. The remaining submatrices are are */
/* now called E12 and E21, respectively. */
begin;
DECLARE (P1, P4) (r,r) FLOAT (18), P2(limc,limc) FLOAT (18),
P3(limr,limr) FLOAT (18);
if n-r-q ^= 0 then
CALL REDUCE (p1, E12, p2, r, limc, t, Tol);
else
t = 0;
if m-r-q ^= 0 then
CALL REDUCE (p3, E21, p4, limr, r, s, Tol);
else
s = 0;
if r = t | r = s then /* Set parameters for no solutions. */
do;
par = 0; go to Endp;
/* Calculations for examples (4) to (7) cease here. */
end;
if s+t = 0 then
do;
CALL EIG (Lambda, A, r); go to Endp;
/* Calculations for examples (8) to (10) cease here. */
end;
if n-r-q ^= 0 then
do;
CALL MATMUL (p1, A, A, r, r, r);
CALL MATMUL (p1, B, B, r, r, r);
end;
if m-r-q ^= 0 then
do;
CALL MATMUL (A, p4, A, r, r, r);
CALL MATMUL (B, p4, B, r, r, r);
/* E12 and E21 have been "reduced" and the requisite */
/* operations have been performed on B. That is */
/* E12 = P1 * E12 * P2 = It 0 */
/* 0 0 */
/* */
/* E21 = P3 * E21 * P4 = Is 0 */
/* 0 0 */
/* and */
/* Br = P1 * Ir * P4. */
end;
end;
end;
end;
/* The columns of Ar above Is and the rows Ar to the left of It are */
/* annihilated, and the remaining (r-s) x (r-t) submatrix is called */
/* G. The corresponding r-s rows and r-t columns of B are called H. */
/* The following statements build the matrices G and H. */
begin;
DECLARE (G, H)(r-t,r-s) FLOAT (18);
do i = t+1 to r;
do j = s+1 to r;
G(i-t,j-s) = A(i,j);
H(i-t,j-s) = B(i,j);
end;
end;
CALL PENCIL (G, H, Lambda, Sp, Par, Tol);
end;
Endp:
;
/* REDUCE applied to an m by n matrix X of rank dex finds an m by m matrix */
/* I1 and an n by n matrix I2 such that */
/* I1 x X x I2 = | Idex 0 | */
/* | 0 0 | */
/* The rank is found by REDUCE and returned in dex. Gaussian elimination */
/* with complete pivoting is used until the (dex + 1)st pivot element */
/* would be less than Tol, a parameter to be supplied by the user to */
/* PENCIL. This procedure is supplied to make the PENCIL routine complete.*/
REDUCE: PROCEDURE (I1, XP, I2, m, n, dex, Tol) OPTIONS (REORDER);
DECLARE (XP, I1, I2) (*,*) FLOAT (18);
DECLARE (m, n, dex) FIXED BINARY (31);
DECLARE Tol FLOAT;
DECLARE (i, j, k, l, lim, p, q) FIXED BINARY (31),
div FLOAT (18);
DECLARE X(HBOUND(XP,1), HBOUND(XP,2)) FLOAT (18);
DECLARE (CVEC, TEMP) (n,n) FLOAT (18), I3(m,m) FLOAT (18);
DECLARE (rvec(m), mvec(n)) FIXED BINARY (31);
X = XP;
if m > n then lim = n; else lim = m;
dex = 0;
do i = 1 to m;
I1(i,i) = 1; rvec(i) = i;
do j = 1 to m;
if i ^= j then I1(i,j) = 0;
end;
end;
do i = 1 to n;
mvec(i) = i;
do j = 1 to n;
if i ^= j then
I2(i,j), TEMP(i,j), CVEC(i,j) = 0;
else
I2(i,j), TEMP(i,j), CVEC(i,j) = 1; /* Missing RHS. */
end;
end;
Rowop:
do i = 1 to m;
do j = 1 to m;
if i=j then I3(i,j) = 1; else I3(i,j) = 0;
end;
end;
dex = dex + 1;
if dex <= lim then
do;
CALL SEARCH (X, dex, k, l, m, n, rvec, mvec);
/* X(k,l) is the pivot element. */
CALL ISWAP (rvec(dex), rvec(k)); CALL ISWAP (mvec(dex), mvec(l));
do i = 1 to n;
CALL SWAP (CVEC(i,dex), CVEC(i,l));
end;
if abs(X(rvec(dex), mvec(dex))) < Tol then
do; dex = 0; RETURN; end;
do i = 1 to m;
if i ^= dex then
do;
div = X(rvec(i),mvec(dex)) / X(rvec(dex), mvec(dex));
I3(i,k) = -div;
do j = dex to n;
X(rvec(i), mvec(j)) = X(rvec(i), mvec(j)) - div*X(rvec(dex), mvec(j));
end;
end;
end;
I3(dex,k) = 1 / X(rvec(dex), mvec(dex));
do j = dex to n;
if j ^= dex then
X(rvec(dex), mvec(j)) = X(rvec(dex), mvec(j)) / X(rvec(dex), mvec(dex));
end;
X(rvec(dex), mvec(dex)) = 1;
if k ^= dex then
CALL SWAP (I3(dex,dex), I3(k, dex));
CALL MATMUL (I3, I1, I1, m, m, m);
do i = dex+1 to m;
do j = dex+1 to n;
if abs(X(rvec(i), mvec(j))) > Tol then go to Rowop;
end;
end;
end;
if m < n | dex < lim then
begin;
DECLARE (p, q) FIXED BINARY (31), mul FLOAT (18);
do i = 1 to dex;
do j = dex+1 to n;
if abs(X(rvec(i), mvec(j))) > Tol then
do;
mul, TEMP(i,j) = -X(rvec(i), mvec(j));
do p = 1 to m;
X(rvec(p), mvec(j)) =
X(rvec(p), mvec(j)) + mul*X(rvec(p), mvec(i));
end;
CALL MATMUL (I2, TEMP, I2, n, n, n);
do p = 1 to n;
do q = 1 to n;
if p ^= q then TEMP(p,q) = 0; else TEMP(p,q) = 1;
end;
end;
end;
end;
end;
end;
CALL MATMUL (CVEC, I2, I2, n, n, n);
END REDUCE;
SWAP: PROCEDURE (r, s);
DECLARE (r, s) FLOAT (18);
DECLARE temp FLOAT (18);
temp = r; r = s; s = temp;
END SWAP;
ISWAP: PROCEDURE (r, s);
DECLARE (r, s) FIXED BINARY (31);
DECLARE temp FIXED BINARY (31);
temp = r; r = s; s = temp;
END ISWAP;
/* Multiplies the u x v matrix X by the v x w matrix Y, depositing the */
/* product in the u x w matrix Z. */
/* Because the result may overwrite one of the arguments corresponding */
/* to X or Y, the intermediate products are temporarily stored in TEMP. */
/* When matrix multiplication is completed, the product in TEMP is */
/* copied to Z. */
MATMUL: PROCEDURE (X, Y, Z, u, v, w) OPTIONS (REORDER);
DECLARE (X, Y, Z) (*,*) FLOAT (18), (u, v, w) FIXED BINARY (31);
DECLARE TEMP (u,w) FLOAT (18), (I, J, K) FIXED BINARY;
do i = 1 to u;
do j = 1 to w;
TEMP(i,j) = 0;
do k = 1 to v;
TEMP(i,j) = TEMP(i,j) + X(i,k)*Y(k,j);
end;
end;
end;
/* Z = TEMP; /* May need to check that dimensions are the same. */
do i = 1 to u;
do j = 1 to w;
Z(i,j) = TEMP(i,j);
end;
end;
END MATMUL;
SEARCH: PROCEDURE (Y, Lim, k, l, m, n, veci, vecj) OPTIONS (REORDER);
DECLARE Y(*,*) FLOAT (18),
(Lim, k, l, m, n) FIXED BINARY (31);
DECLARE (veci, vecj) (*) FIXED BINARY (31);
DECLARE (i, j) FIXED BINARY (31);
k, l = Lim;
do i = Lim to m;
do j = Lim to n;
if abs(Y(veci(i), vecj(j))) > abs(Y(veci(k), vecj(l))) then
do; k = i; l = j; end;
end;
end;
END SEARCH;
%INCLUDE 'ACM-343F.INC';
/* EIG calls EIGENVALUES that finds the eigenvalues of the r x r */
/* matrix X and stores them in the r-vector Lambda. */
EIG: PROCEDURE (Lambda, X, r);
DECLARE (Lambda (*) COMPLEX, X (*,*)) FLOAT (18);
DECLARE r FIXED BINARY (31);
DECLARE XS(r,r) FLOAT (18); /* The extent of X can be greater than r, */
/* hence define XS. */
DECLARE (evi, evr) (r) FLOAT (18);
DECLARE (vecr, veci) (r,r) FLOAT (18);
DECLARE t FLOAT (18);
DECLARE Indic(r) FIXED BINARY;
DECLARE (i, j) FIXED BINARY (31);
/* CALL EIGENVALUES (X, Lambda, r); */
t = places(t);
do i = 1 to r;
do j = 1 to r;
XS(i,j) = X(i,j);
end;
end;
evr = 0; evi = 0;
CALL EIGENP (XS, t, evr, evi, vecr, veci, indic);
do i = 1 to r;
Lambda(i) = COMPLEX (evr(i), evi(i));
end;
Sp = r;
END EIG;
END PENCIL;