Automatic optimization of symbolic algol programs. I. General principles

JOURNAL

OF COMPUTATIONAL

Automatic

PHYSICS

10,

503-533 (1972)

Optimization of Symbolic Algol Programs. I. General Principles M. PETRAVIC AND G. KUO-PETRAVIC

Department

of Engineering

Science, Oxford

University,

Parks Road, Oxford,

England

AND

K. V. ROBERTS Culham Laboratory,

Abingdon, Berkshire,

England

Received June 18, 1971

The symbolic style of programming referred to as Symbolic Algol I [l] appears to have a number of advantages when applied to the solution ,of sets of nonlinear partialdifferential equations. Programs written in that style are clear, elegant, and concise and their modular structure enables large parts of the programs to be used over and over again for many different problems. Such programs, however, tend to be slow because they involve a large number of nested procedure calls at execution time. Finite-difference methods in several dimensions require in general that a relatively small number of equations be solved a large number of times and much is gained if these nested procedure calls are executed only once. This is achieved by a generator or translator program, written in Algol, which processes input written in a related style named Symbolic Algol II. Usually only finite-difference equations in very compact symbolic form are input, while output is completely explicit and can be in a number of computer languages. Of greatest interest are Assembler code modules automatically produced in this way. They are competitive in speed with fully hand-optimized Fortran versions and are produced effortlessly and error free, so that complex sets of equations can readily be programmed or alterations made. For production runs these modules can be incorporated,into a Fortran control program.

1. INTRODUCTION This paper describes how Algol 60 can be used as a powerful macro-processor which enables the symbolic expressions of classical vector analysis to generate efficient target code automatically by the use of controlled side effects. The target lanuages produced so far have been IBM 360 Assembler code, Fortran, Algol and ICL KDF9 Usercode, but it appears that any language might be generated in 503 Copyright All rights

0 1972 by Academic Press, Inc. of reproduction in any form reserved.

504

PETRAVIC,

KUO-PETRAVIC,

AND ROBERTS

a similar way. The target code can be optimized by physical symmetry declarations; for example if I’, = 0 (no rotation) and af/laz = 0 for all functions f (no z dependence), then appropriate declarations can be used to suppress terms in which such quantities occur as products. The method can be used for generalized orthogonal curvilinear coordinates and an example will be given. Finally it would seemthat the method might be extended without difficulty to other kinds of symbolic formalism. A previous paper [l] showed how Algol could be used for the symbolic solution of problems in computational physics, especially those in which sets of partialdifferential equations are solved by finite-difference methods using a discrete mesh. A vector equation such as aB/at = curl(V x B) + qV2B

(1)

can be programmed symbolically in the closely similar form AB[Cl, Q] : = CURL(CROSS(V, B)) + ETA x DELSQ(B);

(2)

a style of programming which has been termed [I] Symbolic Algol I. Here the left side of (2) is an array element representing the magnetic field B in the Cldirection, (Cl = I,2 or 3) at the mesh point Q, while most of the identifiers on the right side are symbolic operators or functions which are closely analogous to their counterparts in Eq. (1) and are representedby real procedures. Details of the choice of coordinate system, the number of dimensions, the boundary conditions, and the difference scheme are excluded from Eq. (2) and are dealt with at a lower level just as in the familiar symbolic notations of mathematical physics. Symbolic Algol I(SA/I) enables complex problems to be coded in a concise form which is virtually system independent and should be readily intelligible to physicists becauseit is close to the mathematical language which they normally use. By way of example, Table I shows the partial-differential equations which are used in the 3D magnetohydrodynamic TRINITY code [l], while Table II TABLE 3D MHD Continuity

equation

Momentum Magnetic

equation equation

Temperature

equation

Equations

I

Used in the TRINITY

Code

ap/at = -v .p~ a(puiyat = - ajax,(pij) + vv2p~i aBlat = v x (V x B) + +ZB aTjar = -V

. (TV) + (2 - y) TV ’ v + KPT

+(Y - lW/p + (Y - 1)4P x 3” + P . $7 Pressure

Pi9 = pT&, +

Current

j=VxB

pvivj + (B2/2)S,, - B,B,

OPTIMIZATION OF SYMBOLIC ALGOL

505

TABLE II 3D MHD EquationsProgrammed in SymbolicAlgol I procedureINVOKE DIFFERENCEEQUATIONS; hegin CONTINUITY EQUATION: DT: = 2 x DELTA T; Cl: = C2: = 1; Q: = 1 + I + 1 + (J + 1) x PI + (K + 1) x PI x PJ; NEW RHO: = RHO - DT x DIV(RH0 x V); MOMENTUM EQUATION: DT: = 2 x DELTA T/(1 + NU/EPS); for Cl: = 1, 2, 3 do AV[CI, Q]: = (RHO x V + DT x (-DIV 2(P) + NU x DELSQ(RH0 x V)))/NEW RHO; ARHO[Q]: = NEW RHO; MAGNETIC EQUATION: DT: = 2 x DELTA T/(1 + ETA/EPS); for Cl: = 1, 2, 3 do AB[CI,Q]: = B + DT x (CURL(CROSS(V, B)) + ETA x DELSQ(B)); TEMPERATUREEQUATION: DT: = 2 x DELTA T/(1 + KAPPA/EPS);Cl: = 1; ATEM[QJ: = TEM + DT x (-DIV(TEM x V) + KAPPA x DELSQ(TEM) + (2 - GAMMA) x SAV(TEM)x DIV(V) + (GAMMA - 1) x (ETA x SQM(CURL(B))/SAV(RHO) + NU x (SQM(CURL(V))+ DIV(V) t 2))); end; real procedureP; P:=ifCI=C2then(RHOx(TEM+VxV)+0.5xDOT(B,B)-BxB)eIse (RHO x V x V2 - B x B2); Note. Thedifferences in the treatmentofp andj betweenTablesI, II, andIII arenot essential. shows the same equations programmed in SA/I. The differences are fairly minor and are partly due to the restricted classof symbols currently available on computer input devices such as the teletype or card punch. The advantages of symbolic notation are clear enough. The Algol procedure operators are neat and concise and have the same formal properties as their mathematical counterparts, so that the manipulation of statements and the construction of new expressionsare quick, intelligible, and easy to check for errors. A typical example is the operator CURL, represented in a Cartesian coordinate system by the short procedure real procedure CURL(A); real A; CURL := RP(DEL(RP(A))) - RM(DEL(RM(A))); 5wo/3-9

(3)

506

PETRAVIC, KUO-PETRAVIC, AND ROBERTS

Here RP and RM are rotation operators which rotate the 1,2, 3-components of vectors or tensors in either the positive (RP) or negative (RM) directions (Fig. l), while DEL is a finite-difference operator. These rotation operators are reciprocal to one another so that RP(RM(A)) = RM(RP(A)) = A,

(4)

RP(RP(A)) = RM(A)

(5)

while the property

is also often used in 3 dimensions. The use of vector and tensor operators ensures that statements are independent of the coordinate system (covariant), the components being hidden and appearing only at execution time.

FIG. 1. Positive and negative rotation operators. The operators RP, RM rotate vector components cyclically in the positive and negative directions, respectively, and satisfy the symbolic relations R+8 = RwS = R+R_ = R-R, = 1, from which R+2 = R- , etc. (Here R, = RP, RI RM.)

All these properties combined with modularity and portability of programs [2] make SA/I a powerful tool for the quick and error-free development of large and complex physics or engineering programs. However, SA/I executes quite slowly because of the great number of nested procedure calls and is therefore not too useful for 2D and 3D production runs although this depends on how well the Algol 60 compiler has been written. Its main application to date lies in the testing of prototype programs on a coarse mesh over a few timesteps. In this way standard test results are obtained for comparison with future better-optimized and faster versions of the same program, written for example in ordinary Algol, Fortran, or Assembler code [l]. The aim of the present paper is to carry the theory of Symbolic Algol one stage further. We shall show that by a further slight transformation of a vector expression such as Eq. (2) it can be made to generate optimized code automatically. In this

OPTIMIZATION

OF SYMBOLIC

507

ALGOL

new style of programming, which is termed Symbolic Algol II (SA/II), Eq. (2) in fact becomes EQUATE(B, SUM(B, MULT(DT, SUM(CURL(CROSS(V, B), MULT(ETA, DELSQ(B))))));

(6)

and the statements of Table II are replaced by those of Table III. The reason for this transformation is to replace the arithmetic operators +, -, x , /, : = which occupy a privileged position in high-level languages by their generalized counterparts SUM, DIFF, MULT, QUOT which are real procedures, and EQUATE TABLE 3D MHD CONTINUITY EQUATE

in Symbolic

Algol II

DIFF(RH0,

MULT(DT(l),

DIV(MULT(RH0,

V)))));

EQUATION:

= 1,2,3do

EQUATE(V,

QUOT(SUM(MULT(RH0,

(DIFF(TEN(B,

B), MULT(RHO(TEN(V,

GRAD(SUM(MULT(RH0, MAGNETIC forC1:

TEM),

V), MULT(DT(2),

DIFF(SUM(DIV

V)))),

DELSQ(MULT(RH0,

MULT(NU,

MULT(RNUM(O.$

DOT(B,

B))))))),

V)))),

NEWRHO));

EQUATION: = 1,2,3do

EQUATE(B,

SUM(B,

TEMPERATURE EQUATE(TEM, (RNUM(2.0), (TEM,

Programmed

EQUATION’:

(NEWRHO,

MOMENTUM forC1:

Equations

III

MULT(DT(3),

SUM(CURL(CROSS(V,

SUM(TEM, GAMMA),

MULT(DT(4),

SAV(TEM)),

SAV(RHO)),

(SQUARE (CWW%

SUM(DIFF(SUM(MULT(MULT(DIFF

DIV(V)),MULT(KAPPA,DELSQ(TEM))),DIV(MULT

V))), SUM(QUOT(MULT(DIFF(GAMMA,

(CURL(B)))),

B)), MULT(ETA,DELSQ(B))))));

EQUATION:

RNUM

MULT@IFF(GAMMA, EXP (DWV),

RNUM

(LO)), MULT(ETA, (LO)), MULT

(NU,

SQUARE SUM

INUMt2))))))))));

Note. Equations (1.5) and (1.6) have been incorporated into the main equations although they can be defmed separately. The 4 calls of the procedure DT take into account the Dufort-Frankel factors [l].

which is a procedure. Once this has been done these procedures can, of course, be given any interpretation that we choose, and they can in particular be made to generate optimized code in any desired programming language by means of side effects as explained in Section 2. The languages generated so far have been IBM 360 Assembler code, optimized

508

PETRAVIC,

KUO-PETRAVIC,

AND

ROBERTS

Algol and Fortran, and ICL KDF9 Usercode which is similar to Reverse Polish and therefore has a theoretical as well as a practical interest. The transformations from Table I and Table II to Table III obey prescribed rules and we have in fact carried them out automatically by two separate methods. One method uses the STAGE 2 Macro-Processor [3,4], and the other usesAlgol60 as a string processor and syntax analyzer [5]. In practice thesetransformations are not difficult to perform by hand. An SA/II generator program looks very like the corresponding SA/I calculational program except that some of the auxiliary statements must also be changed from form (2) to form (6), so that, for example, CURL becomes real procedure CURL(X); real (X); CURL := DIFF(RP(DEL(RP(X))), RM(DEL(RM(X))));

(7)

The purpose of (7) is however rather different from that of (3) becauseinstead of actually calculating numerical values directly, the program now works out which programming instructions are needed to calculate these values and then generates the instructions, either printing them or punching them out on cards or placing them in a file for subsequent execution. Clearly there is a close correspondence between the two processesof calculation and code generation-for example, in both SA/I and SA/II the storage locations are evaluated by manipulating symbolic operators-and so CURL, DEL and other operators appear virtually unaltered. The automatically generated modules can be incorporated into control programs which are written in Fortran, Algal, or any other convenient language, as discussed in Section 6. An SA/II generator program is not a compiler, and it is shorter, more flexible, and easier to write than most compilers. We shall endeavor to explain the relation in Section 6. Since it is coded in Algol and has at least as many nested procedure calls as SA/I it may also not be as fast as some compilers. However, the point is that each statement in Table III is executed only once, instead of many millions of times as in a normal run. Code generation from Table III actually occupies about 10 set on the IBM 360/91, which is comparable with ordinary job overheads and much less than the duration of a typical production run which may last for several minutes or even hours. The only real requirement is that the generated Assembler code should be efficient, and it turns out in practice to be slightly more efficient than the corresponding Fortran written by a good programmer and compiled with the IBM Fortran H Option 2 Compiler. The basic idea of SA/I is quite straightforward and will be explained in Section 2, namely, that side effects of typed procedures can be used to generate code. Some complications arise in a practical program because it is desirable to make the generated code as efficient as possible. As we have already mentioned, it should be

OPTIMIZATION OF SYMBOLIC ALGOL

509

possible in a problem in which V, = 0 or af/az = 0 to declare these symmetry conditions to the generator program and to have it automatically eliminate all terms which are then known to be identically zero. This can be achieved by making use of the fact that the real procedures in a statement such as (6) pass numerical values to one another. The simplest version of SA/II outlined in Section 2 treats these values as dummies, but they can be employed for a variety of purposes including optimization. A zero value is assigned to basic terms or products which are known to be identically zero in the physical problem, and this value is used to control the way in which the code is generated. Unnecessarybrackets and signs are eliminated in a similar fashion. Code generation is carried out by a real procedure TRIPOP (triple operator), and operators such as SUM, DIFF, MULT QUOT simply call TRIPOP to generate the appropriate output. Although TRIPOP is quite short its working is fairly complex, and to avoid burdening the present discussion it will be described in detail elsewhere [6]. A brief account is, however, given in Section 3. Since Algol is likely to be useful for other types of symbolic manipulation it may be appropriate, after some examples have been given, to list those features of the language that we have found to be important, and this is done in Section 8. It will, however, become evident almost from the beginning that recursion and call-by-name play a significant part, together with a number of other facilities which are not available in a language such as Fortran. A more practical point is that although Algol allows identifiers to have arbitrary length, somecompilers (including that of the IBM 360) distinguish only the first six characters. To make the discussion of this paper clearer we have in some casesused long English-language identifiers, but in the published version of the program itself [7] the identifiers are restricted to six characters or less to avoid possible clashes.

2. CODE GENERATION

The generation of code by an expression such as SUM(X, MULT(Y, Z))

63)

is simple to explain. We tist need a procedure PRINT((string)) which will output an arbitrary string of characters, taking into account any special requirements to leave gaps at the beginning and end of the line. In terms of this we define the algebraic symbols procedure PLUS; PRINT(‘+‘); procedure STAR; PRINT(‘*‘);

(9)

510

PETRAVIC,

KUO-PETRAVIC,

ANP

ROBERTS

and the symbols to be used as identifiers in the generated code real procedure X; begin PRINT(‘X’); X : = 1; end; real procedure Y; begin PRINT(‘Y’); Y := 1; end; real procedure Z; begin PRINT(‘Z’); Z := 1; end;

(10)

Finally, there are two arithmetic operators real procedure MULT(A, B); real (A, B); begin real CALL; CALL : = A; STAR; CALL := B; MULT := 1; end; (11) and real procedure SUM(A, B); real(A, B); begin reai CALL; CALL : = A; PLUS; CALL : = B; SUM : = 1; end;

(12)

Explanation Exactly what happens is illustrated by the tree shown in Fig. 2. When the statement ‘CALL := A’ is encountered, the real procedure X is invoked by name from SUM, causing a symbol ‘X’ to be generated. The local variable CALL in

PRINT

i ‘Y’

f ‘*’

i ”2

FIG. 2. Generation of code by the expression SUM(X, MULT(Y,

a.

SUM is set equal to 1 but this value is disregarded. A ‘ + ’ symbol is next produced. Finally, MULT is entered by ‘CALL := B’ which generatesthe symbols ‘Y’, ‘*‘, ‘Z’ in the same way so that we get x+y*z as required. Several points may be noticed: (a) Typed procedures are used to generate the code.

(13)

OPTIMIZATION OF SYMBOLIC ALGOL

511

(b) The actual generation takes place by means of side effects; the statement CALL := A; causes the string of symbols associated with the formal parameter A to be constructed, however complex an expression A mqy represent. (c) A precise order is forced on these side effects, independently of the order in which the Algol compiler writer may choose to evaluate the terms in an arithmetic expression such as P + Q or P x Q. (d) It is immaterial whether real or integer procedures are used, but we have chosen real procedures for compatibility with SA/I. (e) At this stage the values associated with the real procedures have no significance; only the side effects are important. In Section 3 we shall however indicate how a suitably chosen assignment of values can be used to remove expressions which are known to be identically zero. Generation of Indices

The variable indices defining a particular component of a vector variable, and the mesh point at which it is to be evaluated, are specified by global variables whose values are controlled by the cyclic rotation operators RP(X), RM(X) and the displacement operators EP(X), EM(X). The action of these operators relies on the Algol “call by name” facility [S, 91.In the example real procedure RP(X); real X; begin Cl := CP[Cl]; RP := X; Cl := CM[Cl]; end;

the procedure X is made to use the value of Cl defined within the procedure RP, becauseit is evaluated at the time the statement RP := X; is executed and not at the time of the call to RP, as would be the case if the “call by value” alternative were used. Generation of Reverse Polish

The reader might be excused for wondering what exactly has been achieved so far: starting from the algebraic expression (13) we have converted it automatically or by hand into theI SA/II form (88), and then proved that this is capable of reproducing the original expression. The real advantagesare first, that the generated code can represent a considerable expansion of the original which may be in symbolic vector or tensor form, and second, that by making small changes in the basic procedures such as PLUS, STAR, X, Y, Z, MULT, SUM we can generate the code in a lower-level and therefore more-efficient language.

512

PETRAVIC,

KUO-PETRAVIC,

AND

ROBERTS

To illustrate this second point we show how to convert (8) into the Reverse Polish form x, y, z, *, +,

(14)

which is in l-l correspondence with the Assembler Language (or “Usercode”) of the ICL KDF9 on which Symbolic Algol was first developed. To do this we simply modify the five procedures (9) and (10) in order to add an extra comma, e.g., procedure PLUS; PRINT((‘+,‘);

(15)

and then reverse the order of the second and third statements in the arithmetic operators (11) and (12), e.g., CALL := A; CALL := B; STAR; MULT := 1;

(16)

which then enables (8) to generate (14). Brackets

The question of brackets has not so far been taken into account in generating algebraic expressions such as (13). Thus at present x x (Y + Z)

(17)

MULT(X, SUM(Y, Z))

(18)

would be converted by into X*Y+Z,

(19

which is wrong, although the corresponding Reverse Polish string is still correct: x y, z, +, *,

(20)

The simplest way of correcting this is to enclose the output from each arithmetic operator in a pair of brackets by means of two procedures OPENB, CLOSEB, so that, for example, (12) becomes real procedure SUM(A, B); real(A, B); begin real CALL; OPENB; CALL := A; PLUS; CALL := B; SUM := 1; CLOSEB; end;

W)

Then expressions (8) and (18) generate the correct output m + gc * ZN

(22)

OPTIMIZATION OF SYMBOLIC ALGOL

513

(X * (Y + m

(23)

and respectively. A complicated expression will now produce a large number of unnecessarypairs of brackets but although these make the expression difficult for a human to understand they should not trouble a Fortran or Algol compiler; in fact they might well shorten the compilation time since questions of operator precedenceneed no longer be resolved. However, in order to make the generated code more elegant and intelligible we have implemented a method for removing the unnecessarybracket pairs and this will be described in Ref. [6]. Another difficulty which occurs with the KDF9 is that the nesting store or arithmetic stack has finite depth (in practice about 13), so that a sufficiently long string of identifiers unrelieved by operators could causeit to overflow. To monitor this situation (which has not yet occurred in the problems that we have treated) one can introduce a level counter and print out a warning to the user when overflow is detected. In most caseshe can then permute the order in the SA/II expression to produce a better pattern. For example SUM(MULT(Y, Z), X)

(24)

y, z, *, x, +,

(25)

instead of (8) leads to which requires a smaller stack. Overflow is more likely to occur with the IBM 360 which has only 4 floating-point registers, and here we perform automatic dumping and restoring of registers when necessaryas explained in Ref. [6]. 3. PROCESSING THE EQUATIONS The central procedure of the SA/II generator program is EQUATE, which for Algol output is procedure EQUATE(X, Y): real X, Y; begin real CALL; NEXTLINE; PRINT := false; if X = 0 then go to EXIT; PRINT := true; SIGN := 0; (26) CALL := X; ASSIGN; if Y = 0 then TEXT (1, ‘0’); SEMICOLON; EXIT: end;

514

PETRAVIC,

KUO-PETRAVIC,

AND

RbBERTS

A similar version is used for any other target code but the explanation is simpler for this particular example and may give some idea of the power of the symbolic method. ASSIGN generatesthe assignment symbol, ‘= .’ for the IBM 360. Suppose that (26) is called by the statement EQUATE(B, SUM(B, MULT(DT, CURL(CROSS(V, B)))));

(27)

as in Table III but without the resistive diffusion term. Because this is a vector equation it will be called 3 times with Cl = 1,2, 3 and we shall suppose that B, has been declared identically zero. The tree of real procedure calls is indicated in Fig. 3. EQUATE

FIG. 3. Proceduretree for the magnetic equation (27). Eacharithmeticoperatorindicated by * calls the TRIPOP operator recursively to control the optimization process. Branches broken off with // are similar to the full branch shown on the diagram. Terminal symbols V, B, DT, R2DS generate the code, while RP, RM manipulate the global component variable Cl, and EP, EM manipulate the mesh location.

The auxiliary call to NEXTLINE does shy editing that is required, e.g., ruling a line across the output page, and printing is then’ switched off. The left side X of the expression, in this caseB, , B, or B, , is next tested to seeif it is identically zero, which as explained in Section 1 will be represented by a zero value. If so,

OPTIMIZATION

OF SYMBOLIC

ALGOL

515

there is no point in generating anything and the right side is skipped. (Cl = 3). Note that since X is a real procedure the statement if X = 0 then ... involves further operations at a lower 1evel.l Assuming that the component is not identically zero we switch printing on again, make sure that no + sign will be printed, and referenceX once more, subsequently calling the procedure ASSIGN. This seriesof actions will generate some expanded form such as

BlUQ/).=

(28)

depending on the hardware representation and on the way in which array components are being referenced.Note that the numerical value in (28) is that of the formal parameter V. The whole of the testing and generation of the right side is now contained in the deceptively simple statement if Y = 0 then TEXT( 1, ‘0’);

(29)

Some possible casesfor the right side Y are a. Constant. Since printing is now switched on the numerical representation of the constant is generated. b. Nonsubscripted variable. The character string corresponding to this variable is generated. c. Array variable. If this variable has been declared to be identically zero it will not be printed and the value 0 will be returned, causing the second part of the if statement to print ‘0’. If it is not identically zero it will print its own character representation, including any necessaryvector or tensor components. d. Arithmetic expression. The outermost arithmetic operator calls TRIPOP which switches off the printing and tests the whole of the expression, TRIPOP being called again recursively by each of the internal arithmetic operators, including those in CURL and CROSS. If it is finally determined to be identically zero then ‘0’ is printed. Otherwise this outermost TRIPOP initiates a second scan which prints all those subexpressions which were found to be not identically zero the first time they were tested. As an example, Table IV shows Algol 60 target code generated for the IBM 360 for Maxwell’s two equations aE/at = curl H,

aH/Ct = -curl E,

1 Since X is real it is perhaps not good practice to test whether or not it is exactly zero but we have found no difficulty on the ICL KDF9 or IBM 360. The procedure values are always set to integers or sums or differences of small integers and it is hard to see why trouble should arise with any computer. If it does, one can simply test for X < E where 0 < E < 1.

516

PETRAVIC,

IWO-PETRAVIC,

AND

ROBERTS

TABLE 1V Maxwell’s Equations in 3D Orthogonal Curvilinear Coordinates Generated in IBM 360 Algal” ‘COMMENT’----

____

---------

________________;

HFILDl(/Q/). = HFILDl(/Q/) - DCT7 * ((EFILD3(/Q + 141) * H3(/ + 2/) - EFILD3(/Q-1 41) * H3(/ - 20) * R2DQ2 - (EFILDZ(/Q + 140/) * H2(/ + 31) - EFILDZ(/Q - 140/) * H2(/ - 3/)) * R2DQ3) * RHPHMl; ‘COMMENT’---------------------------------; HFILD2(/Q/). = HFILD2(/Q/) - DCT7 * ((EFILDl(/Q + 140/) * Hl(/ + 3/)-EFILDl(/Q-14 Oi) * Hl(/ - 31)) * R2DQ3 - (EFILD3(/Q + l/) * H3(/ + l/) - EFILD3(/Q - I/) * H2(/ - l/)) * R2DQl) * RHPHMZ; ‘COMMENT---------------------------------; HFILD3(/Q/). = HFILD3(/Q/) - DCT7 * ((EFILD2(/Q + l/) * H2(/ + l/) - EFILDZ(/Q li) * H2(/ - l/)) * R2DQl - (EFILDl(/Q + 14/) * Hl(/ + 2/) - EFILDl(/Q - 14j) * Hl(/ - 2/)) * R2DQ2) * RHPHM3; ‘COMMENT---------------------------------; EFILDl(/Q/). = EFILDl(/Q/) + DCTS * ((HFILD3(/Q + 141)* H3(/ + 2/) - HFILD3(/Q-1 4/) * H3(/ - 20) * R2DQ2 - (HFILD2(/Q + 1401)* H2(/ + 3/) - HFILD2(/Q - 140)/ * H2(/ - 3/)) * R2DQ3) * RHPHMI: ‘COMMENT’---------------------------------; EFILD2(/Q/). = EFILDZ(/Q/) + DCT8 * ((HFILDl(/Q + 1401) * Hl(/ + 3/) - HFILDl(/Q-14 O/) * Hl(/ - 31)) * R2DQ3 - (HFILD3(/Q + l/) * H3(/ + l/) - HFILD3(/Q - l/) * H3(/ - l/)) * RZDQl) * RHPHMZ; ‘COMMENT’---------------------------------; EFILD3(/Q/). = EFILD3(/Q/) + DCT8 * ((HFILDZ(/Q + l/) * H2(/ + 10 - HFILD2(/Q l/) * H2(/ - l/)) * R2DQl - (HFILDl(/Q + 141) * Hl(/ + 2/) - HFILDl(/Q - 14/) * Hl(/ - 2/)) * R2DQ2) * RHPHM3; a In this example, PI = 14 and PJ = 10, hence the displacements in the 1, 2, 3 directions are (1, -1), (14, -14) and (14 x 10, -14 x lo), respectively.

in orthogonal curvilinear coordinates using a mesh of size (14 x 10 x 42). The source statements were for Cl = 1, 2, 3 do EQUATE(EFIELD, SUM(EFIELD, MULT(DCT(8), CURL(HFIELD)))); for Cl = 1,2,3 do EQUATE(HFIELD, DIFF(HFIELD, MULT(DCT(7), CURL(EFIELD)))); (30)

OPTIMIZATION OF SYMBOLIC ALGOL

517

In this caseno symmetry conditions were imposed. DCT7 and DCT8 are time-step factors which will usually be the same.

4. PHYSICAL CONSTANTS AND VARIABLES Provision is made for handling nonsubscripted variables which may be constants or functions only of the time, and scalar, vector, and tensor functions. A typical “declaration” specified by the user is, for example, real procedure B; B : = VECTOR(5, 1, ‘B’, 1, 1,O);

(31)

This calls the real procedure VECTOR which in the Algol60 output version reads real procedure VECTOR(ORDINAL NUMBER, LENGTH, S, Vl, V2, V3); integer ORDINAL NUMBER, LENGTH, VI, V2, V3; string S; begin VECTOR := if Cl = 1 then Vl else if Cl = 2 then V2 else V3; if not PRINT then go to EXIT; SIGN IT; TEXT(LENGTH, S); COMPONENT; SHIFT; EXIT: end; (32) The parameter ORDINAL NUMBER is not being used here; it controls the actual storage region used by B which is important in some Assembler Language versions and is retained for consistency. LENGTH gives the length of the identifier string ‘B’ which will be used in target statements, instructions, and comments, in this case 1. The last 3 parameters VI, V2, V3 define whether or not the x, y, z components are identically zero. In this case the user has specified V3 = 0 so that B, is taken to be identically zero. Because VECTOR returns the value 0 whenever Cl = 3, the z component of the magnetic equation will be suppressedaltogether (Section 8), and all terms in which B, occurs as a product will be suppressedon the right side of any other equation. Otherwise when PRINT = true the procedure VECTOR will proceed as follows: SIGN IT TEXT(LENGTH,

Examine the value of a global variable SIGN and output either a preliminary ‘+‘, ‘-‘, or no sign. S) Output the string S of length LENGTH, example ‘B’.

in our

518

PETRAVIC,

KUO:PETRAVIC,

AND

ROBERTS

COMPONENT

Output the current value of the component Cl, e.g., ‘2’.

SHIFT

Output ‘[’ or ‘(/’ followed by the signed numerical value of the current displacement from the local mesh origin, and finally a closing bracket, ‘I’ or ‘/)‘.

Typical output for Algol on the IBM 360 is

+WQ

+ 730

(33)

and similarly for Fortran except for the absence of the slashes. Several other versions have been developed; for example,, the displacements can be handled symbolically so that the code does not have to be regenerated when the mesh size is changed or one might generate ‘BY’ instead of ‘B2’ for clarity. The operators RP, RM act on the global variable Cl as determined by algebraic and analytic operators such as DOT, CROSS, DIV, CURL and therefore enable the correct component label to be calculated. Similarly, vector translation operators EP, EM which are called whenever a spacederivative occurs are used to calculate the position on the lattice relative to the central point Q of the local mesh ‘molecule’. In Eq. (33) this appears as a numerical displacement of 73 words within the region of core store occupied by the array B2. A scalar is handled in a precisely similar way by the real procedure SCALAR which requires only a single value 0 or 1 and does not need to call COMPONENT, and correspondingly for constants on the one hand, and tensors on the other. Adressing of Variables

We have already mentioned that the generated code can be in any chosen computer language. Whatever the language, however, we obviously have to be able to print some form of address at which the value of a particular variable defined on a specific mesh point is stored. The addressing can be completely symbolic, e.g., in Algol via an array B[CI, Q + DX - DY],

(34)

it can be coded numerically as in KDF9 Usercode V34P6M 15

(35)

(location 34 of variable block 6 modified by register 15), or it can be partly numerical and partly symbolic as in IBM 360 assemblycode. 0052(,DISPQ)

BY(I,J- 1,K)

0312(PLUSDZ,DISPQ)

BZ(I,J,K+ 1)’

(36)

OPTIMIZATION

OF SYMBOLIC

ALGOL

519

Here the displacement in bytes defines both the variable and the position in the xy plane relative to the centre of the molecule, DISPQ is a general register which contains the current center, while PLUSDZ is a register which shifts one mesh unit in the z direction. The comments on the right side give the Fortran notation. A fuller example of IBM 360Assembly code produced in this way is shown in Table V. Several other choices are possible. In any case the address will depend on the variable in question, on the coordinate direction on to which a vector variable is being projected, and on the meshpoint at which the expression should be evaluated. These 3 quantities are defined in all the Symbolic Algol programs [1] by the variable name and by the global integer variables Cl and Q. Finite-DifSerence Mesh Used for TRINITY

The mesh is assumed to be Cartesian, of up to 3 dimensions, and equispaced. The numbers of mesh points in the x, y, and z directions are denoted by PI, PJ and and PK, respectively, so that the total number of mesh points is SIZE = PI x PJ x PK. This includes six guard planes introduced to enable the samedifference expressions to be used on the physical boundaries of the volume as inside it. The mesh points are counted in the x direction, starting along the intersection of the first y- and the first z-plane numbered by integers Z, J, and K. The numerical relation is Q = I + (J - 1) x PI + (K - 1) x PI x PJ.

(37)

For output in Algol or Fortran we choose separate arrays for each vector component in which the values are packed with increasing number Q. All the arrays are then of the samesizeequal to PI x PJ x PK and are functions only of one index Q. This way of mapping the 3-dimensional mesh onto l-dimensional arrays implies that two mesh points which are adjacent in the y direction correspond to array elements PI storage locations apart, and that if the points are adjacent in the z direction the corresponding elements will be PI x PJ locations apart. The single mesh index changes by DELTA = PI x PJ if the shift is in the z direction. It is the vector-displacement operators EP and EM that cause,on a single application, a shift by one mesh point in the direction determined by the current value of Cl. The number of shifts in the x, y, and z directions is denoted in the program by the global integers Kl, K2, and K3 from which, knowing PI and PJ, the corresponding change in Q can be calculated. It is sometimes desirable to interleave the variable arrays so that the variable values at one mesh point Q are stored sequentially in the memory, e.g., in the order RHO, Vl, V2, V3, Bl, B2, B3, TEM

(38)

520

PETRAVIC,

KUO-PETRAVIC, TABLE

AND ROBERTS

V

IBM 360 Assembly Code Generated for the Second Component of the Magnetic Equation The code is efficient, with no extra address-manipulation instructions. Out of the 48 instructions, 4 are concerned with stack overflow. Note that comments are also generated automatically.

Operation code

Location

Variable

*

Level 00

LE

0,0308(, DISPQ)

LE

2, CONSTANT

LE

4,0296(PLUSDZ,

B2+0+0

001

01

002

DT

02

003

DISPQ)

v2+0+0

03

004

+ 008

ME

4,0312(PLUSDZ,

DISPQ)

B3+0+0

03

005

LE

6,03C@@‘LUSDZ,

DISPQ)

v3+0+0

04

006

ME

6,0308(PLUSDZ,

DISPQ)

B2+0+0

04

007

SER

4,6

03

008

LE

6,0296(MINSDZ,

DISPQ)

v2+0+0

04

009

ME

6,03 12(MINSDZ,

DISPQ)

B3+0+0

04

010

SER

496

03

011

LE

6,03OO(MINSDZ,

DISPQ)

v3+0+0

04

012

ME

6,0308(MINSDZ,

DISPQ)

B2+0+0

04

013

03

014

AER

496

ME

4, CONSTANT

RECDS2

03

015

LE

6,0324(, DISPQ)

v1+1+0

04

016

ME

6,0340(, DISPQ)

B2+1+0

04

017

STE

0, STORAGE

05

018

LE

0,032s

v2+1+0

05

019

ME

0,0336 (, DISPQ)

Bl+l+O

05

020

SER

690

04

021

LE

0,0260(, DISPQ)

Vl-1+0

05

022

ME

0,0276(, DISPQ)

B2-110

05

023

04

024

05

025

+ 016

+ 000

LE

60 0,0264(, DISPQ)

v2-1+0

ME

0,0272(, DISPQ)

Bl-l+O

SER

AER

60

ME

6, CONSTANT

SER

4,6

+ 016

RECDS2

05

026

04

027

04

028

03

029

OPTIMIZATION

OF SYMBOLIC

521

ALGOL

TABLE V (continued) Operation code

Variable

Location

ETA

04

030

DISPQ)

B2+0+1

05

031

0,0052(, DISPQ)

B2+0-1

05

032

DISPQ)

B2+0+0

05

033

DISPQ)

B2+0+0

05

034

B2+1+0

05

035

B2-110

05

036

06

037

06

038

06

039

032

05

040

690

04

041

04

042

03

043

02

044

02

045

01

046

01

047

01

048

LE

6, CONSTANT

LE

0,0X4(,

AE AE

0,0308(PLUSDZ,

AE

0,0308(MINSDZ,

AE

0,0340(, DISPQ)

AE

0,0276(, DISPQ)

STE

2, STORAGE

LE

2, = E’ + 6.0000’ + 00

ME

2,0308(, DISPQ)

SER MER ME

6, CONSTANT

AER

436

LE

2, STORAGE

MER

294 0, STORAGE

LE AER STE

Level

+ 004

+ 004 B2+0+0

+ 020

REDSSQ

+ 004 + 000

092 0,0308(, DISPQ)

B2+0+0

followed by the samesequencefor the next meshpoint (Q + 1) and so on. Adjacent values of the samevariable are then (say) 32 bytes apart. This enablesphysical data planes to be transferred readily to and from the backing store as a single block, and can be done either in Assembler code or, in Fortran, by the use of EQUIVALENCE statements.In all casesthe SA/II program can readily be adapted to calculate the correct word or byte position in the store. A similar situation exists when the core store is too small to hold the complete set of physical data. The largest size of mesh on which TRINITY has been run has 60 x 80 x 48 points and requires two IBM 2301 drums to accommodate the 7 Mbytes of data. All the calculation then takes place within three sectors of a rotating quadruple buffer, the fourth sector being used to transfer data to and from the core store in parallel with the calculation. The SHIFT procedure has been adapted so that it always refers to the correct areas of buffer storage. 581/10/3-10

522

PETRAVIC, KUO-PETRAVIC, AND ROBERTS 5. ORTHOGONAL CURVILINEAR COORDINATES

Symbolic vector algebra and analysis on a Cartesian lattice in SA/I have been described elsewhere [l] and few changes are required for SA/II. It may, however, be of interest to demonstrate how the method has been extended to generalized orthogonal curvilinear coordinates, with spherical polar coordinates as a special case. The generalized definitions for divergence and curl can conveniently be written as div A = &

c +, (hi+&Ai), z a

ai cur1 A = c -jq&= z

a(hi-Ai-) &Ii+

3(hi+Ai+) aqi’

where (a,, a,, a3) are unit vectors, h = (h, , h, , h3) are scale factors, and +, - denote positive and negative cyclic rotation, respectively. These are translated into SA/II as real procedure DIV(A); real A;

DIV := MULT(DOT(DEL(MULT(A,

HPHM)), R2DQ), RHlH2H3);

(41)

real procedure CURL(A); real A;

CURL := MULT(DIFF(RP(MULT(DEL(RP(MULT(A, H))), R2DQ)), RM(MULT(DEL(RM(MULT(A, H))), R2DQ))), RHPHM); (42) where real procedure DEL(X); real X; DEL := DIFF(EP(X), EM(X));

(43)

This leaves to be defined the real procedures H RHl H2H3 RHPHM R2DQ

-+ (h, , h, , h3) (vector function) --f l/(h,h,h,) (scalar function) (vector function) + l/(h+h-) + l/(2 * DQ) (vector).

(44)

The mnemonic ‘R’ means reciprocal and signifies that division is avoided in the interests of efficiency, and for similar reasonsh,h,h, and hfh- are defined separately instead of being constructed from h. When using the leapfrog scheme [l] we find it useful to store the scale factors and quantities that are constructed from them on a subsidiary mesh centered on

OPTIMIZATION OF SYMBOLIC ALGOL

523

the point Q, which may either contain 7 points if only one displacement &Aqi occurs at a time (Fig. 4) or 27 if they occur together. Then, for example, RHPHM becomesa real procedure which generatesthe code RHPHMl [I], RHPHM%[I], RHPHM3[1]

(45)

with I in the range (-3,3) as in Table IV. The generated target module is therefore still independent of the coordinate system although it does depend on the symmetry. To run the target program we must reload the variables on the subsidiary mesh whenever they alter. In spherical polars the scale factors are h,=

FIG.

1,

ho = r,

h, = r sin 8,

(46)

4. Subsidiary mesh used to store the scale factors.

so that we can minimize the amount of recalculation by assigning the variables in order of q1 = q.-+q2 = 8, q3 = r with the innermost scan over q1 , although this is not necessarily the best choice on other grounds. Some further improvements can be made if the coordinate system is specified at generation time; for example h, might be automatically suppressedsince it is known to be unity.

6. THE PROBLEM PROGRAM AND THE GENERATOR PROGRAM

The problem program may be in any desired combination of languages, and by no means all of it need be constructed automatically. One convenient approach is to use Fortran for the more straightforward control sections of the program, and for setting up core storage, and to use automatically generated assembly-code modules to handle the solution of the partial-differential equations themselves as

524

PETRAVIC,

KUO-PETRAVIC,

AND

ROBERTS

well as any complexities due to boundary conditions. Changes to the physics can then be made reliably and quickly, simply by generating a new module and linking it to the remainder of the program. The automatically generated modules operate on variables in the COMMON storage area, and communicate with the rest of the program either via COMMON or through argument lists. To switch to a new type of computer system, it is only necessaryto alter the generator program so that it outputs the appropriate assembly code. An alternative approach is to use a combination of Algol 60 control program and assembly-code modules. This approach has advantages if the local Algol system is sufficiently powerful, but in most cases Fortran is at present to be preferred. One therefore starts by developing a strategy for the problem which takes into account its modular structure, the layout of core and backing storage, the communication between program modules, the choice of language and programming style for each module, and so on. Those parts of the program whose execution count is low should be optimized according to considerations other than those of CPU efficiency; for example, they should be made intelligible or easy to write and to debug. Portability should also be planned right from the start, so that low-level modules designed for one computer system can be rapidly replaced by those written for another. The Generator Program

An Assembler code replacement for an SA/I module of the prototype program is constructed by incorporating the corresponding SAjII module into the generator program which is then run (Fig. 5). We visualize a series of a related physical problems, run either on one computer at a single laboratory or on a range of different computers at several laboratories. If the generator program has been constructed properly only a small part of the work has to be repeated for each new situation. This saves time and makes it easier for one person to understand a range of programs, since they share a family resemblancelike members of related biological species. The generator program ought to be highly modular since there are a number of requirements which may well need to be changed independently. Examples are a. Computer System on which the code is generated. i. Algol character representation. ii. System output procedures. b. Computer System for which the code is generated. i. Character representation.

OPTIMIZATION OF SYMBOLIC ALGOL

ii. iii. c. d. e. f.

525

Language. System facilities.

Number of dimensions and coordinate system. Mesh and storage layout Symbols used for variables and constants. Numerical methods and difference schemes.

xx PHYSICAL PROBLEM

/

-

I

Express algerithm symboblly Plan progmm structure Work out storage

FIG. 5. Programming strategy. The method of solution is first tried out by writing a prototype program which is used to produce sets of test results. Modules of the prototype are next reconstructed either by hand or with the SA/II generator program to make them more efficient. Test results from the production program are then compared with those from the original tests to make sure that the updating has been carried out correctly.

526

PETRAVIC, KUO-PETRAVIC, AND ROBERTS

g. Physical equations. h. Boundary conditions. i. Physical coefficients (e.g., thermal conductivity). Figure 6 shows the relation between the modules which have been developed so far while Table V provides a list. In Section 7 we shall go through this list briefly, indicating what the various modules do. A detailed description together with a program listing and test runs will be published elsewhere [6]. Relation to Compilers and Macro-Generators

In mathematics or theoretical physics it is always possible for an author to devise a new notation or extension to the accepted“language” and having defined it for the reader, to use it throughout a paper or a course of research.This flexibility has been largely unavailable in computing science,where it has been customary to use rather standardized and limited languages such as Fortran, Algol or PL/I which are constructed either by manufacturers or by international committees, and translated by compilers which are expensive and time-consuming to write and to maintain. The only degree of freedom left to the individual has been the ability to define sets of library subroutines or procedures which in effect become additions to the standard notation. There is no limit on the extensions which can be achieved in this way but programs tend to run slowly if they make considerable use of procedure calls, as in SA/I. Macro-processors such as STAGE2 [4] enable any string of symbols to be given a meaning. The user is free to define sets of macros which convert any string into any other string and eventually, into the code of some high- or low-level language whose efficiency depends solely on his own ingenuity. Thus the full flexibility of mathematics is achieved provided that the character set is wide enough. SA/II appears to lie somewhere in between. The formal structure of the input string is constrained by the syntax of Algol so that there is usually an excessive number of brackets as in Table III, but the manipulations that can be carried out are quite general and it is remarkably easy for the individual user to make alterations or additions to the “language” by changing the basic procedures of the generator program. An SAjII generator program is also quite short; usually only a few hundred Algol statements.Thus we have effectively at our disposal an ultra high-level “language” compiler which is difference-schemeand problem dependent, but is also easily changed by any user. 7. STRUCTURE OF THE GENERATOR PROGRAM

The current version divides logically into 6 main modules and 15 submodules as shown in Table VI.

527

OPTIMIZATION OF SYMBOLIC ALGOL

VI.2 CONTROL STATEMENTS

)a V6

INITlALl2AlUl

FIG. 6. Relation between the modules of the generator program. An arrow indicates that one module makes use of procedures belonging to the other. Note that all output is channelled through a single short module BASIC OUTPUT so that only this module has to be changed to runlthe program on another computer.

I. BASIC OUTPUT

In transferring

the generator

program

to a new computer

system

the iirst task is

528

PETRAVIC,

KUO-PETRAVIC, TABLE

List of Generator

AND

ROBERTS

VI

Program

Modules

and Submodules

I. BASIC OUTPUT II. CHARACTERS III.

MATHEMATICS III.

A

1. ARITHMETIC

111.2. ALGEBRA IV, MATHEMATICS IV.1.

VECTOR

B ANALYSIS

IV.2. SPACE AND V. OUTPUT

TIME

SCALES

ORGANIZATION

V. 1. TRANSLATION

LOGIC

V.2. SIGNS, SPACES, NUMBERS V.3. VARIABLE

AND

FUNCTIONS

CLASSES

V.4. COMPONENTS

AND

VS. COMMANDS

AND

DERIVATIVES REGISTER

CHECKS

V.6. INITIALIZATION VI. PROBLEM

DEFINITION

VI.1. PHYSICAL

CONSTANTS

VI.2. CONTROL

STATEMENTS

VI.3. SOURCE

AND

VARIABLES

STATEMENTS

to rewrite this module, which forms a link to the standard Algol output procedures of the computer on which the generator is being run. (This need not of course be the same as the computer for which the optimized code is being produced.) The ICL KDF9 version occupies about 25 cards and it can usually be i-written for another system in a few hours. The module defines the output channel and sets up standard formats, and contains procedures which enable the output to be manipulated in a straightforward system-independent way; e.g., Output N blank spaces BLANKS(N) Start a new line LINE OUTNUM(L, F, X) Output the value of an arithmetic expression X in format F, length L Output a string S of length L? TEXT(L, S) The complete set is enough to generate code and comments in any language. e The simpler procedure PRINT discussed in Section 2 did not contain the parameter helps the output procedures to organize the layout of the line.

L, which

OPTIMIZATION

OF SYMBOLIC

ALGOL

529

II. CHARACTERS This also contains no submodules. It is made up of a dozen or so simple statement procedures which enable one to refer to symbols like (, + ; = etc., by symbolic names: OPENB, COMMA, PLUS, SEMICOLON, EQUALS. Although this somewhat slows down code generation it makes the program much easier to read and quicker to write, particularly in view of the awkward way in which string quotes are often represented in Algol. Examples are procedure EQUALS; TEXT( 1, ‘=‘); procedure CLOSEB ; TEXT( 1, ‘)‘); procedure OPENSB; TEXT(1, ‘[‘) procedure MINUS; TEXT(1, ‘-‘);

(47)

The first argument gives the lenth of the string, in this case1. III. MATHEMATICS A The tist mathematics module comprises the submodules ARITHMETIC ALGEBRA. The former of these contains a set of procedures SUM(X, Y) DIFF(X, Y) MULT(X, Y)

MULT21(X, Y)

QUO’WL Y) SUM3(X, Y, Z)

which are all dealt with by a single fairly complex procedure TRIPOP (triple operator) to be explained elsewhere [5], e.g., real procedure SUM3(X, Y, Z); real X, Y, Z; SUM3 := TRIPOP(5, X, Y, Z, 1); (48) The first argument of TRIPOP is an operation code, the second to fourth are the operands, and the last two specify the first or secondindexes (in the caseof a secondrank tensor). Procedure SUM3 is convenient when handling 3D scalar products and divergences while MULT21 is used for handling tensor contractions. The others deal with the arithmetic operations +, -, x , /. The ALGEBRA submodule contains the rotation operators RI? and RM and the vector-algebraic operators DOT and CROSS.

530

PETRAVIC,

IV. MATHEMATICS

KUO-PETRAVIC,

AND

ROBERTS

B

In the simplest version the second mathematics module contains the two submodules 1. CARTESIAN ANALYSIS LEAPFROG 2. SPACE AND TIME SCALES. As its name implies, the first of these submodules depends on the mesh geometry and on the difference schemesused (although not on the physical problem or on the computer system). The procedures EP, EM, DEL, GRAD, DIV, CURL, SAV, DELSQ are fairly direct translations of the SA/l versions already published [I]. The other submodule deals with the constants dt, 2dS and (&)z and can be also readily extended. V. OUTPUT

ORGANIZATION

This module dependson the output language. It may contain up to 6 submodules of which number 5 is omitted in the Algol target version: Contains most of the logic needed to generate V.l. TRANSLATION LOGIC. the output code and to eliminate expressions which are identically zero, as well as unnecessary brackets. It therefore deservesa more detailed description which is given in Ref. [6]. A number of standard V.2. SIGNS, SPACES NUMBERS AND FUNCTIONS. utilities are provided here, some of which depend on the output language or format. For example a line overflow in Fortran requires that a continuation symbol should be punched in column 6; overflow in assembly language is handled in a different way while at the end of an Algol statement a semicolon is required. Signs, integer and real numbers and elementary mathematical functions are also provided for. The most important procedure is EQUATE which has been discussed in Section 3. V.3. VARIABLE CLASSES. Contains procedures which deal with constants, scalars, vectors and tensors of which (32) is an example.

Contains the procedures COMV.4. COMPONENTS AND DERIVATIVES. PONENT and SHIFT which generate the code for referencing storage locations, including any subsidiary calculations that are needed. V.5. COMMANDS AND REGISTER CHECKS (IBM 360 Assembler code only). Contains procedures REGISTER REGISTER, REGISTER STORAGE which issue IBM 360 instruction mnemonics, and STORE IF OVERFLOW

OPTIMIZATION OF SYMBOLIC ALGOL

531

which determines whether or not the required register is already in use, if so copying it into a reserved location in core store. V.6. INITIALIZATION. Contains a procedure START which initializes the variables of the generator program. VI. PROBLEM

DEFINITION

This module is provided by the user and consists of three submodules which have been kept as simple and as close to the physics as possible: 1. PHYSICAL CONSTANTS AND VARIABLES 2. CONTROL STATEMENTS 3. SOURCE STATEMENTS. Table III gives the source statements for TRINITY while (31) is an example of a variable declaration. Typical initialization statements are NDIM := 3; PI := PJ := PK := 8;

(49)

(use an 8 x 8 x 8 mesh in 3 dimensions). 8. CONCLUDING REMARKS

A satisfactory solution to the slownessof programs written in Symbolic Algol I has been found. Using a style known as Symbolic Algol II it has been possible to translate finite-difference equations automatically into fully explicit codes in a number of target languages. The best of these codes are fully competitive in speed with hand-optimized Fortran and are fast enough to make the solution of timedependent magnetohydrodynamic equations in three space dimensions a feasible proposition. In a recent exercise,a 3D plasma code in rotating spherical coordinates was designed and written in about four days using the SA/II method. The same translator program can be used for other systemsof fluid equations and also for problems in two dimensions. Although the advantages of the method are smaller when applied to simpler problems, repetition of effort can still be avoided. More importantly, the use of well-tested procedures reduces programming errors, in particular those of a numerical nature which are often impossible to detect except experimentally through a comparison with another calculation. Though at first sight trivial this may prove to be one of the important attractions of the method. The underlying principle is that instead of writing a problem program by hand, one constructs a generator problem which writes it automatically. Because this generator program is built up from prefabricated modules and is also highly symbolic it can be developed and altered very quickly.

532

PETRAVIC, KUO-PETRAVIC, AND ROBERTS

In essencewe are using Algol as a powerful macro-generator which is capable of substituting one expression into another as well as performing many subsidiary calculations. Becausea value is associatedwith each substitution, extra information can be carried along which allows some optimization to be done. A further extension might be to relate this value to a generalized variable type (e.g., logical, integer, real, complex, quaternion, matrix or whatever), so that any necessary conversions can be carried out and the basic operators SUM, DIFF, MULT, etc., can be interpreted in the appropriate way in each case.This is close to the procedure which is followed in mathematics which allows operators such as +, -, x to be freely generalized to new classesof object. Other interesting possibilities are to apply the SA/II technique to other kinds of program such as operating systems and compilers, and to other types of computer such as the CDC STAR. The features of Algol 60 that appear to be necessaryor useful for this kind of work are [8,9] (a) Call by name. Needed for the symbolic substitution of one expression into another. (b) Parameterless typedprocedures. Just as in mathematics, a function need have no explicitly indicated arguments. (Fortran does not allow this.) (c) English-language ident$ers of any length, with blanks ignored. Can be used to make programs more intelligible and to avoid bulky comments. (d) Elimination of the unnecessary word ‘CALL’. (e) Ability to have several statements on one line. Both (d) and (e) make programs more concise and attractive. (f) No overhead on procedure declarations. Often these declarations are only one card long, and one line in the compiler listing, instead of several pages as in Fortran. (g) Block structure for variable scopes. Global variables can be passedinto a procedure implicitly without the need for a bulky COMMON deck or argument lists. (h) Recursion. Typed procedures can be substituted into one another without restriction as required by the mathematical physics. (i) Side efsects. Available also in other languages, but mentioned here as being crucial to the whole method. ACKNOWLEDGMENTS The early planning of SA/II was carried out in collaboration with Dr. J. P. Boris. We should like to thank Dr. F. Hertweck and his colleagues at the Institut fiir Plasmaphysik, Garching,

OPTIMIZATION OF SYMBOLIC ALGOL

533

Federal Republic of Germany, for making available to us the excellent facilities of the IBM 360/91 Computing Centre at the Institute. We should also like to thank Mr. R. S. Peckover for many discussions on Symbolic Algol, and Dr. N. K. Winsor for providing an improved compiler.

&WRENCE?S

1. K. V. ROBERTSAND J. P. BORIS,The solution of partial differential equations using a symbolic style of Algol, J. Comput. Phys. 8 (1971), 83. 2. K. V. ROBERTS AND R. S. PECKOVER, Symbolic programming for plasma physicists, in “Proceedings of the Fourth Conference on Numerical Simulation of Plasmas,” U.S. Naval Research Laboratory, Washington D.C., November 1970, p. 165. Office of Naval Research, Department of the Navy, July 1971. 3. G. KUO-PETRAVIC,M. PETRAVIC,AND K. V. ROBERTS,The translation of Symbolic Algol I to Symbolic Algol II by the STAGE 2 Macro-Processor, Culham Laboratory preprint CLM-P 275. 4. W. M. WAITE, The mobile programming system STAGE 2, Comm. ACM 13 (1970), 415. 5. G. KUO-PETRAVIC,M. PETRAVIC,AND K. V. ROBERTS,The translation of Symbolic Algol I to Symbolic Algol II using Algol60, to be submitted for publication in Comput. Phys. Commun. 6. M. PETRAVIC,G. KU~PETRAVIC, AND K. V. ROBERTS,Automatic optimization of Symbolic Algol programs. II. Code generation, to be published. 7. M. PETRAVIC, G. KUO-PETRAVIC,AND K. V. ROBERTS,The Symbolic Algol II generator program, to be submitted for publication in Comput. Phys. Commun. 8. B. HIGMAN, “A Comparative Study of Programming Languages,” McDonald/Elsevier Computer Monographs (1967). 9. R. BAUMANN, M. FELICIANO, F. L. BAUER, AND K. SAMIJELSON,“Introduction to Algol,” Prentice-Hall, Englewood Cliffs, NJ., 1964.