/***************************************************************************/
/* */
/* ACM Algorithm 365 - COMPLEX ROOT FINDING */
/* */
/* H. Bach, Laboratory of Electromagnetic Theory, Technical University of */
/* Denmark, Lyngby, Denmark. */
/* */
/***************************************************************************/
/* Transcribed to PL/I by R. A. Vowels, 7 March 2006. */
/* UNTESTED AS YET. REQUIRES COMPLEX ARITHMETIC. */
/* MORE TESTING IS REQUIRED, INCLUDING CHECKING AGAINST THE ORIGINAL */
/* FORTRAN PROGRAM. */
/* This subroutine determines, within a certain region, */
/* a root of a complex transcendental equation f(x) = 0, */
/* on which the only restriction is that the function w = f(x) must be */
/* analytic in the region considered. The iterative method used, */
/* the downhill method, was originally described in H. Bach, */
/* "On the downhhill method", in Communications of the ACM 13, */
/* Dec. 1969, pp. 675-677 and is discussed and modified in J. A. Ward, */
/* "The downhill method of solving f(x) = 0", in Journal of the ACM 4, */
/* March 1967, pp. 148-150. */
/* The program uses a complex function FUNC(Z) for the */
/* computation of f(x). From a given complex starting-point ZS, */
/* the iteration is performed in steps of initial length HS. The */
/* iterations stop at the root approximation ZE when either the */
/* function value DE at the end point is less than the prescribed */
/* minimum deviation DM or when the step length HE has become */
/* less than the prescribed minimum step length HM. For reference, */
/* the subroutine also returns DS, the function value at the starting */
/* point ZS, and N, the number of iterations used. There are thus */
/* four input parameters, namely the starting point ZS, the initial */
/* step length HS, the minimum step length HM, and the minimum */
/* deviation DM. */
/* Translator's comments:
/* 1. The conditional statement at label L6 has two paths */
/* chosen by GO TO statements. Therefore the pair of */
/* statements on the following line can never be executed. */
/* 2. A subscript error is obtained with I = 0 at label L16. */
/* A patch has been put in temporarily. */
(SUBRG, FOFL, SIZE):
CRF: PROCEDURE (ZS, HS, HM, DM, DS, ZE, HE, DE, N);
DECLARE (ZS, ZE) FLOAT COMPLEX,
(HS, HM, DM, DS, HE, DE) FLOAT,
N FIXED BINARY;
/* INCOMING: ZS = starting point (complex), */
/* HS = initial step length, */
/* HM = minimum step length, */
/* DM = Minimum deviation. */
/* OUTGOING: DS = Deviation at start */
/* ZE = the root approximation (complex), */
/* DE = the deviation at the end point, */
/* HE = the step length at the end point, */
/* N = number of iterations used. */
DECLARE (Z0, ZD, ZZ, Z(3), CW, A, V,
U(7) STATIC INITIAL (
1,
.8660254+0.5I,
.0+1I,
.9659258+.2588190I,
.7071068+.7071068I,
.2588190+.9659258I,
-.2588190+.9659258I )) FLOAT COMPLEX;
DECLARE W(3) FLOAT;
DECLARE (H, W0) FLOAT;
DECLARE (I, K, NR) FIXED BINARY;
ON SUBSCRIPTRANGE SNAP BEGIN;
PUT SKIP DATA (I);
PUT SKIP DATA (NR);
END;
H = HS;
Z0 = ZS;
N = 0;
/* Calculation of DS. */
CW = FUNC(Z0);
W0 = ABS(REAL(CW)) + ABS(IMAG(CW));
DS = W0;
IF W0 <= DM THEN GO TO L18;
K = 1; I = 0;
L2: V = -1;
/* Equilateral triangular walk pattern. */
L3: A = -0.5+.866I;
/* Calculation of deviations W in the new test points. */
L4: Z(1) = Z0 + H*V*A; /* >>>>>>>>>> was Z(I) <<<<<<<<<<<<<< */
CW = FUNC(Z(1));
W(1) = ABS(REAL(CW)) + ABS(IMAG(CW));
Z(2) = Z0 + H*V;
CW = FUNC(Z(2));
W(2) = ABS(REAL(CW)) + ABS(IMAG(CW));
Z(3) = Z0 + H*CONJG(A)*V;
CW = FUNC(Z(3));
W(3) = ABS(REAL(CW)) + ABS(IMAG(CW));
N = N + 1;
PUT SKIP DATA (N);
PUT SKIP DATA (Z(1));
PUT SKIP DATA (Z(2));
PUT SKIP DATA (Z(3));
/* Determination of W(NR), the smallest of W(I). */
IF W(1) > W(3) THEN GO TO L6;
IF W(1) >= W(2) THEN GO TO L8;
L6: IF W(2) <= W(3) THEN GO TO L8; ELSE GO TO L9;
NR = 1; GO TO L10;
L8: NR = 2; GO TO L10;
L9: NR = 3;
L10: IF W0 < W(NR) THEN
DO;
IF K = 1 THEN /* Reduction of step length */
DO; K = 2; IF H < HM THEN GO TO L16; H = H/4; GO TO L3; END;
ELSE IF K = 2 THEN /* Restoration of step length. */
DO; K = 3; H = 4*H; GO TO L2; END;
ELSE /* Rotation of walk pattern. */
DO; I = I + 1; IF I > 7 THEN GO TO L17; END;
END;
K = 1; I = 0;
/* Forward directed walk pattern. */
A = .707+.707I;
V = (Z(NR) - Z0)/H;
W0 = W(NR);
Z0 = Z(NR);
IF W0 > DM THEN GO TO L4; ELSE GO TO L18;
L16: put skip list ('At label L16. I = ', i);
if I = 0 then I = 1; /* A patch. <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
V = U(I); GO TO L3;
/* Reduction of step length. */
L17: IF H >= HM THEN
DO;
H = H/4;
I = 0; GO TO L2;
END;
L18: ZE = Z0;
HE = H; DE = W0;
RETURN;
END CRF;