AUTOMATIC GENERATION OF SOLUTION PROCEDURES
FOR INDEXED EQUATION SETS USING GENIE
By
GARY LOUIS ALLEN
A DISSERTATION PRES:NTKD TO THE GRA JA'\E CO^i(CII. O
THE UiIIVERSily OiF FLORIDA
S ,ARTiA! FULFIL!E' :" OF THE RP,UIR ,M F ?:TS FO
DEGREE OF DOCTOR CF PiILOSOPHY
UNIVERSITY OF FLORIDA
1974
To Jan
I
ACKNOWLEDGEMENTS
The author wishes to thank the chairman of his supervisory
committee, Dr. A. W. Westerberg, Associate Professor of Chemical
Engineering, for suggesting the research topic and for his guidance
through the course of the research. Thanks are also extended to
Dr. S. W. Director, Associate Professor of Electrical Engineering, and
Dr. F. D. Vickers, Associate Professor of Computer and Information
Sciences, for serving on the supervisory committee.
The financial support of the Graduate School and of the College
of Engineering is gratefully acknowledged. Thanks are also extended
to the National Science Foundation for the financial support afforded
by grants GK18633 and GK41606. In addition to assistantship
support, these grants provided computer funds for the development of
the programs in GENIE. The author accepted these funds with mixed
feelings, as he believes that in a free society research should be
supported by private concerns, rather than by the goi'ernet.
The a.ithor extends special thanks to nis wife, Jan, iho not only
provijdd encouragement when it was ,, os:t neded1, but also served as
copy editor and typist for this dissertation.
TABLE OF CONTENTS
ACKNOWLEDGEMENTS....................................
ABSTRACT..............................................
CHAPTERS:
1 INTRODUCTION..................................
2 INDEXED EQUATIONS.............................
2.1 Conventions Used in Describing Indexed
Equations ...............................
2.2 FunctionVariable Incidence Matrices....
2.3 Index Display Matrices..................
2.4 Index Imbedding.........................
2.5 Output Set Assignments .................
2.6 Acceleration of Variables..............
3 SOLUTION PROCEDURES FOR INDEX EQUATION SETS...
3.1 Definitions.............................
3.2 An Example Problem......................
3.3 Output Sets and Decision Variables......
3.4 Index Imbedding ........................
3.5 Decoupling .............................
3.6 Blocking Factors...............................
4 ALGORITHMS FOR DERIVING SOLUTION PROCEDURES FOR SETS
OF INDEXED EQUATIONS ...............................
4.1 Decoupling..................................
4.2 Index Output Set Assignments ..................
iii
vi
4.3 FunctionVariable Output Set Assignment.........
4.4 Minimum Weighted Tearing..... .................
5 COMPUTER IMPLEMENTATION ...........................
5.1 Data Structures.............................. .
5.2 A Versatile Memory System .....................
5.3 Flow Diagram and Subroutine Descriptions........
6 CONVERGENCE PROPERTIES ..................................
6.1 Review.........................................
6.2 Modification of Solution Procedures............
6.3 Convergence Properties of Combined
GaussSeidel and NewtonRaphson................
7 EXAMPLE PROBLEMS .............................. ......
7.1 Distillation Model .............................
7.2 Examples ........................... ......... .
7.3 Discussion................................... .
8 CONCLUSIONS AND RECOMMENDATIONS......................
APPENDICES:
A THE SIM DATA STRUCTURE..............
B SUBROUTINE DESCRIPTIONS.............
BIBL IOG APHY..................................
BIOGRAPHICAL SKETCH ...........................
Page
75
77
80
80
88
96
102
102
105
107
114
114
1 !
123
131
135
139
158
160
. .. .
. w D t
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
AUTOMATIC GENERATION OF SOLUTION PROCEDURES
FOR INDEXED EQUATION SETS USING GENIE
By
Gary Louis Allen
June, 1974
Chairman: Arthur W. Westerberg
Major Department: Chemical Engineering
Chemical engineering design and simulation problems generally
involve solving simultaneously a large number of equations. Many
algorithms are available for structurally analyzing a set of equations
with the goal of developing a solution procedure for the equations.
Some of these algorithms have been implemented in a package of
computer programs called "GENDER," an acronym for General Engineering
Design Routines.
Many chemical engineering processes are ma.iherr,atica'iiy mindeled
with indexed equations. The index ranges are often largo, resulting
in very large numbers of equations to be solved simultaneously.
Conventional algorithms, such as those in GENDER, require that each
equation to be analyzed by explicitly represented. This in itself is
a difficult chore for indexed equations with large index ranges. A
second shortcoming of existing algorithms is the time required to
analyze such a large set of equations.
The GENIE system, which is an acronym for GENDER with Indexed
Equations, provides a more efficient means of analyzing sets of indexed
equations. A compact means for representing the structure of indexed
equation sets is presented. This representation allows both existing
and newly developed algorithms to be applied in order to generate a
solution procedure for the indexed equations. An evaluation of the
convergence characteristics of the proposed solution procedure is made,
using order of magnitude estimates of variable values. Should a
proposed solution procedure be found to be nonconvergent, means are
provided for modifying it to enhance the likelihood of convergence.
Formulas are presented which allow the estimation of the convergence
properties of the new solution procedure.
The solution procedures generated are independent of the ranges
on the indices. Conventional methods derive solution procedures for
a particular set of equations. The solution procedures are generated
in the form of F0RTRAN "doloops," so the index ranges become the
doloop parameters and can change from application to application.
Solution procedures have been derived for a typical chemical
engineering problem, a mathematical model of a disti llation column.
Two different problems were solved by the solution procedure. GEN;I
successfully predicted convergence difficulties for or= of the problems
and modified the solution procedure. The modified solution procedure
did solve the problem.
CHAPTER 1
INTRODUCTION
The design or simulation of a chemical process requires finding
the solution to a large set of equations (many of them nonlinear).
Oftentimes these equations represent either a physical recycling of
a process stream from one unit to another, or an interaction among
variables within a unit; both conditions force a simultaneous solution
of the equations. Since the equations are, in general, nonlinear, an
iterative solution is necessary. Frequently convergence of the itera
tive solution depends on the manner in which the solution is carried
out. Unless the engineer has special familiarity with the particular
set of equations he is trying to solve, he essentially must choose
blindly from among the possible methods apparent to him. An alternative
to this is to develop algorithms suitable for implementation on a
digital computer capable of analyzing the structural and numerical
properties of the set of equations with the aim oF finding a means of
solving the equations which is not only efficient, but also likely to
converge to the solution. The method chosen to solve a set of equations,
whether chosen by the engineer or analytical algorithms, is called the
"'solution procedure."
Two commonly employed types or solution procedures are tearing
procedures (such as GaussSeidel) and "iNewton" type procedures (such as
NewtonRaphson). Because of the cost (usually in computer time) of
solving large sets of equations simultaneously, it is desirable to avoid
NewtonRaphson type solution procedures whenever possible. The alterna
tive, tearing procedures, requires that certain variable values be
guessed or torn (these variables become the iterative variables to be
converged), all variable values are calculated algebraically including
the torn variables, and the new values of the torn variables are used
for the next iteration. The method continues until the criteria for
convergence are satisfied.
A tearing solution procedure requires an output set assignment
for the algebraic solutions. An output set assignment is the assign
ment to each equation of one of the variables occurring in that
equation is such a way that each variable is assigned to only one
equation. Not all variables need be assigned to equations; the
unassigned variables are "decision" variables. Early studies of the
nature of output set assignments were concerned with structural proper
ties (Steward, 1962,1965). Later studies addressed the problem of
assigning an output set with the aim of enhancing the likelihood of
convergence (Edie, 1970) and developed effective algorithms for assign
ing outputs when some criterion for optimality has been defined (Gupta,
1972; Gupta et al., 1974).
An iterative solution procedure is made iterative by the presence
of tear variables which must be converged. iMany algorithms have been
developed for choosing tear variables in order to satisfy any of several
criteria, such as minimum number of torn variables or minimum sum of
weights of torn variables, where each variable is given a weight
(Sargent and Westerberg, 1964; Christensen and Rudd, 1969; Barkley and
Motard, 1972; Ramirez and Vestal, 1972; Pho and Lapidus, 1973).
The choice of decision variables has received comparatively less
attention although the problems of simplifying the resulting solution
procedure (Lee, Christensen, and Rudd, 1966) and avoiding singular
systems (Edie, 1970) have been treated.
The problem of deciding in what order to solve the equations has
been solved (Sargent and Westerberg, 1964) to the extent that acyclic
(nonrecycle) systems can be discovered and equations involved in
various recycle calculations can be identified.
All of the algorithms developed thus far require an explicit
representation of each equation and variable before the analysis can
begin. For units described by a large number of equations (for example,
a 50 plate distillation column with three components, which results in
350 equations) the total number of equations becomes quite large, and
application of the above algorithms becomes time consuming and requires
large amounts of computer storage space. In chemical engineering design
and simulation the proliferation of equations and variables is quite
often dlia to staged processes, which are modeled with indexed equations
and variables, which are the same on every stage except for index value.
Great savings of time and computer storage are realized by analyzing
the set of indexed equations, each equation writtt.n enly once, rather
than having to expand the set of equations and write separate equations
for different index values. Unfortunately., existing algorithms are not
directly applicable to these indexed eqaLtioris. Another shortcoming of
existing algorithms is that the solution procedures derived for indexed
equations (written in expanded form) would not be valid for an arbitrary
range on index values.
Algorithms for generating solution procedures for nonindexed
equation sets have been programmed for use on digital computers by
Cunningham (1972, 1973). The package of computer programs is called
GENDER, for General Engineering Design Routines. The GENDER system is
designed to provide an integrated computer program package for auto
matically generating solution procedures.
The goal of this dissertation is to provide to the GENDER system
a means of generating solution procedures for sets of indexed equations.
The computer programs developed to do this are called the GENIE system,
GENIE being an acronym for GENDER with Indexed Equations. Chapters 2
to 4 detail modifications both to existing algorithms and to the means
of representing the set of equations which allow an automatic derivation
of a solution procedure.
The way that the data are represented within the computer storage,
i.e., the data structure, has a profound effect on both the total
storage required and on the time required to execute an algorithm.
Chapter 5, along with Appendix A, discuss the various data structures
used.
Although the algorithms developed in Chapters 2 to 4 allow a
solution procedure based on structural considerations to be derived,
there is, unfortunately, no guarantee that the solution procedure will
converge unless the system of equations fully precedence orders, i.e.,
is completely acyclic. It is desirable to base a solution procedure
not only on structural considerations, but also on numerical consider
ations. Chapter 6 discusses some of the relevant numerical consider
ations and gives some insight into the prediction of convergence
properties.
5
Chapter 7 gives some worked example problems, including distil
lation problems, which illustrate the importance of numerical consider
ations. Conclusions and recommendations are given in Chapter 8.
CHAPTER 2
INDEXED EQUATIONS
In order to discuss indexed equations a terminology must first be
developed. This is especially true of the indices themselves since no
existing terminology seems well suited to them. This is done in the
next section of this chapter. Also, because it is desirable to have
algorithms which efficiently derive solution procedures it was necessary
to restrict both the types of equations allowed and the types of
solution procedures which can result. Since it has been the aim to
treat equations describing physical systems, those equations are among
those allowed. The remainder of the chapter is devoted to further means
of representing the equations and to a discussion of the restrictions
imposed on the problems.
2.1 Conventions Used in Describing Indexed Equations
The set of material balance equations described dhove might be
written as:
M 1 +v .1 .v .f. = 0 21
3i I1 ij i,] i+l 1,j ij
where i = stage number
i = component number
M = name of function type, i.e., Material balance
., = liquid flow rate of component j leaving stage i
v. = vapor flow rate of component j leaving stage i
f.. = feed flow rate of component j entering stage i.
1i7
Thus, M.. represents the steady state material balances for all
components on all stages.
The name of an indexed function is used to distinguish functions
of one function type from functions of another function type. Therefore,
material balance functions, designated "M," and energy balance functions,
designated "E," would belong to different function types. Similarly,
variables occurring in indexed functions are classed by "variable type."
For example, component liquid flow rate is of a different variable type
from component vapor flow rate. It should be noted that in modeling a
multiunit process, material balance equations in a given unit contain
different function types from material balance equations in another
unit, even though the units are similar, e.g., distillation columns.
It is convenient to make a distinction between the indices which
subscript a function type and those which subscript a variable type.
Consider equation 21. The indices "i" and "j" which subscript the
function type M are said to be function indices. Similarly, the indices
which subscript the variables 1, v, and f are said to be variable
indices.
In order to structurally analyze a set of equations it is not
necessary to know the exact form of the equations. It is sufficient to
kno;: wh.ih variables occur in which equations. All existing algorithms
which perform structural analysis use only this information. This being
the cas, ai further simplification of the representation of indexed
euuations is made, with the ain of facilitating i2e structural analysis
of the equation set.
The convention is adopted that function indices will be designated
by the letter "i." Thus, M. would become M,. ., the subscript on the
SJ 1112
L
subscript being used to distinguish between the two different function
indices. Variable indices are designated by the letter "j." Since in
the original representation of the function,the indices on the variables
were functions of the indices on the functions, the variable indices
must be expressed as functions of the function indices.
First, some restrictions will be made on the possible values which
the function indices can assume. Function indices will always have
their possible values defined by a "range" which consists of a lower
limit, an upper limit, and an increment. An index range is very similar
to a FORTRAN doloop range and consequently makes the eventual writing
of the solution procedure in the form of doloops convenient. Examples
of function index ranges are:
a) ii: L = 1 b) i2: L = 2 c) i3: L = 1
U = 10 U = 29 U = 13
A= 1 A = 1 A = 2
In case a) the function index would have possible index values of all
integers between 1 and 10 inclusive, in case b) i2 would have possible
index values of all integers between 2 and 29 inclusive, and in cazs c)
13 would have possible index values of all odd integer values between
1 and 13 inclusive. Although restricting the function indices to this
form reduces the possible values that function indices in a function
index set can have, it does not exclude any of the index sets an
engineer is likely to create in modeling a chemical process. The stages
in a staged process are numbered consecutively; a distillation column
which has a third and fifth stage has a fourth stage. Certainly the
index sets which arise from staged processes and from discretization of
models are describable in terms of index ranges.
More flexibility is desired in describing variable indices than
is available in describing function indices. Variable index sets have
their possible values defined by a range similar to the function index
range or by a "list." The variable index range differs from the
function index range in that its upper and lower limits will be
permitted to be functions of either the function index value itself or
the function index range limits. This restriction allows for many
different arrangements of indexed variables. The equation
f. = 0 = x. +x x i=1,2,...10 22
1 ii" +1
represented by a tridiagonal incidence matrix, would have its index
sets described as:
i: L = 1 j: L = ii
U = 10 U = i+l
A= 1 A=
where both the lower and upper limits of the variable index range are
defined as offsets from the function index value. Other examples of
variable index ranges are:
ji: L = L. j2: L ij3: L = L.
U = U U. U U.
i 1
A =1 A = 1 A 1
These; three examples illustrate the other three basic forms of the
variable index range. The first, jil, defines a lower triangular
incidence matrix, the second, j,, defines an upper triangular incidence
matrix, and th third, J3, defines a full incidence mritrix.
A second means of defining variable indices, the variable index
list, allows even greater flexibility. Suppose we have an indexed
equation of the following form:
f. = x. +x.x 23
S1i2 1 1+1
The variable indices for this function cannot be expressed as an index
range depending on the function index. In order to describe the index
relationships in an equation of this type, we use a list. The variable
index list is a list of offsets from the function index value; it has
arbitrary length. The variable index list for equation 23 would have
length 3 and would be:
j: ii2
il
There is some overlap between index lists and ranges. For example,
either a list or a range can be used to describe the tridiagonal
relationship in equation 22. The index list and range provide for a
structural description of a wide range of indexed functions, and in
particular seem to describe any which are likely to occur in modeling
any chemical process.
2.2 FunctionVariable Incidence Matrices
Since the incidence matrix is widely used to display structural
information about a set of equations, it would be desirable to be able
to represent a set of indexed equations in an incidence matrix.
Unfortunately, for a typical problem such as a 20 stage, 4 component
distillation column the number of rows would be 180 and the number of
columns would be even greater. Incidence matrices of this size make
analysis by hand time consuming and error prone. Even computer analysis
becomes expensive, especially when much larger problems are encountered.
A convenient, compact representation of the structure of a set of
indexed equations is afforded by the "FunctionVariable Incidence
Matrix" (FVIM). The FVIM is not a true incidence matrix, and because
of that has some unusual properties. A FunctionVariable Incidence
Matrix is a matrix whose rows each correspond to a function type and
whose columns each correspond to a variable type. An element a.. of
an FVIM exists if the variable type corresponding to j occurs in the
function type corresponding to i, otherwise the element does not exist.
An element which exists does not assume a value in the conventional
sense. Instead, the element assumes the names of the variable indices
which correspond to that variable type's incidence in the appropriate
function type. As an example consider the following set of equations:
f = 0 = x +x y. +y 24
1112 i11,i2 1ii2 112 t12+1
g = 0 = x. +y y
gil1i 112 i li 2 +1,1i2
The FunctionVariable Incidence Matrix would be:
x y
g(il,i2) 3?j4 1 324
it
__ I
where I iL1 J2: il h3" O
iL i]+!
J4: i2 Js: i2
i +2
In the position corresponding to function type f and variable type x
is the element jij4. From the definitions of the variable indices we
can deduce that in function type f we would find variables x. il 2 or
more compactly noted, variables x..
S13 L
I
It should be noted that variable index names cannot be assigned
casually. An example will serve to point out a possible pitfall. The
problem occurs when there is more than one variable index. Consider
the following set of equations:
f. = 0 = x +x. y +y 25
1112 111,i2 1li2 ,1i2 11,i2+1
One might be tempted to assign four variable index names as follows:
iJl i1 j2:' il J3: i2 j4: i2
i ii+1 12+1
This, however, does not accurately reflect the structure of the
equations 25 for f. does not contain the variables y. Variables
11i2 j2j4
y j2j include yi. i+y. y ,i2. ,y.i+,i2"1, the second and
third of which do not occur in equation f. This error can be
1112
avoided by several methods. First note that the problem arose from
attempting to define a single set of variable indices for the variable
y. Since this is not possible in this case, this equation could not
be represented as it is in a FunctionVariable Incidence Matrix. By
rewriting the variables y as y and z, this shortcoming can be overcome.
It should be noted that this type of problem is unlikely to occur, and
does not occur in any of the several chemical process models examined
in the course of this work.
The FunctionVariable Incidence Matrix then provides the means
for compactly displaying the structure of a wide variety of sets of
indexed equations and variables. Further, with slight modifications,
virtually all sets of indexed equations and variables can be included.
The FVIM will subsequently be shown to be valuable in effecting a speedy
structural analysis of large sets of indexed equations.
2.3 Index Display Matrices
The method used until now to describe variable indices, i.e.,
lists and ranges, is not well suited for a structural analysis of the
relationship between function indices and variable indices. Since many
existing algorithms are either designed to treat incidence matrices, or
can be easily modified to do so, it is desirable to somehow convert the
lists and ranges to incidence matrix representation. To do this the
rows of the incidence matrix are made to correspond to the values of the
function index and the columns are made to correspond to the values of
the variable index. An element a.. is zero unless the variable index
137
can have the value corresponding to the column j when the function index
is equal to the value corresponding to the row i. The resulting matrix
is called an "Index Display Matrix" (IDM).
The number of rows in an Index Display Matrix is determined by the
range for the function index upon which it depends. Thus,.for the
function index ii and the variable index j, which depends on it defined
as follows:
ii: L = 1 ji(ii): iu1
U = 20 ii
A = 1 il+i
the Index Display Matrix (with zero elements blank and nonzero elements
denoted by 'x') would be that of Fig. 21. This representation of the
IDM occupies considerable space and is certainly larger than is neces
sary to convey the structural pattern either to a person studying it or
to an algorithm analyzing it. In addition, as the range of the function
index changes, so does the Index Display Matrix. In order to reduce the
size of Index Display Matrices to something small enough to allow
I I 1 I 1 I I i I 2 2
0 1 2345 6 7 890 12 345678 9 0
I X X X
4 XXX
2 XXX
3 XXX
5 XXX
5I ~XXX
6 XXX
7 XXX
8 XXX
9i XXX
10 X XX
V \e V
3 XXX
14 X XX
15 XXX
16 X XX"
17 X X X
18 X, XX
2O1 ^ ^
FIGURE 21
INDEX DISPLAY MATRIX
efficient structural analysis and to free the Index Display Matrix from
its dependence upon the function index range the concept of the
"blocking factor" is introduced. The blocking factor for a function
index is the number of rows which will be included in the IDM for
variable indices which depend upon that function index. If the function
index ii above were assigned a blocking factor of 3 then the IDM for j1
would be:
0 1 2 3 4
1 x xx
2 X X X
2 x x x
3 x xx
This Index Display Matrix conveys the notion of a tridiagonal matrix
while occupying considerably less space both on paper and in any data
structure used in computer implementation. Since the blocking factor is
arbitrary and is defined by the person defining a problem it should be
pointed out that the blocking factor must be large enough to avoid any
ambiguity in describing the IDM. In particular a blocking factor of 1
should never be used.
2.4 Tndex Imbedding
The value of FunctionVariable Incidence Matrices and Index
Display Matrices is that while they are considerably smaller than the
actual incidence matrix which they represent, they contain all of the
inorrmation contained in the incidence matrix. This can best be illus
trated by the fact that the incidence matrix can be constructed from the
Function Variable Incidence Matrix and its Index Display Matrices.
First it is necessary to introduce the concept of index ordering.
In naming the rows of an incidence matrix for a set of indexed equations
it is desirable to have an orderly way to "step through" the various
index values. The reasons for this are first, to insure that each
function is represented and that no functions are included more than
once,and second, to provide a logical means of stepping through the
various functions which is easily adaptable to an iterative solution
procedure. The method adopted is similar to the nesting of FORTRAN
"doloops." There is one important deviation, however, in that the
function type is also treated as an index. The reasons for this modifi
cation become clear upon examination of the problem in Fig. 22. If we
consider the three indices (function type, ii, and i2) as doloop
indices, obviously there are six ways in which these indices can be
nested. This gives rise to six different structures for the indidence
matrix. Here we have assumed that the ordering of the functions which
would be dictated by the doloop would be the order of the functions
occurring in sequential rows. Also, the variables would be ordered
similarly, i.e., each variable index would have the same nested position
as the function index upon which it depends and the variable type would
be assumed to depend upon function type. While it is obvious that the
index I precedes the index 2 in a normal ordering, there is no such
normal ordering when considering the function or variable type. It is
not known whether f precedes g or whether x precedes y and at this point
it does not really matter and the functions and variables will be put
into an arbitrary order. Later the actual order will be dictated by
precedence ordering considerations, similar to those applied to non
indexed equations and variables. The description of the index orderings
FunctionVariable Incidence Matrix
x y
a) f(il,i2) jlj2 J4j3
g(il,i2) jlj3 jlJ3
Index Definitions
i,: L = 1
U=3
A = 1
Jl: L = 1
U = ii
A=
J4: L = 1
U = U.
ii
i: L = l
j2: i21
FIGURE 22
AN EXA:''LE PROBLEM
U=2
A=
i3: i2
i2+1
Index Display Matrices
J2 J4
123 123
b) 1 x c) 1 x xx
i 2 x x i1 2 x x x
3 x x 3 x x x
J2 J2
012 123
0 1 2 1 2 3
d) 1 x x e) 1 x x
i2 2 x i2 2 xx
FIGURE 22 (cont.)
(i.e., the order of nesting) in terms of doloops makes the transition
to computer implemented solution procedures in FORTRAN using doloops
straightforward. The advantages of using the doloop in a solution
procedure are primarily the compactness of code and the speed of
execution. Indeed, most solution procedures implemented in order
to solve indexed equations use doloops or code similar to doloops.
In order to fully realize the consequences of the nesting of the
indices and the effect of the FunctionVariable Incidence Matrix and
Index Display Matrices on the actual incidence matrix, examine Fig. 23.
This is the expanded incidence matrix for the problem in Fig. 22.
First look at Fig. 23 a), the entire incidence matrix. The matrix is
partitioned into four submatrices, the divisions being chosen along the
boundaries between function types and between variable types. These are
the outermost nested "indices" and hence the ones examined first. The
structure of the partitioned matrix is that of a 2x2, nonzero elements
in all positions. That is the same as the structure of the FVIM, which
has its rows and columns labeled by function and variable type. The
next matrix, actually the (1,1) partition of Fig. 23 a), is itself
partitioned along boundaries between the indices nested next innermost
(corresponding to i2). Figure 22 d) exhibits the s;me structure as
Fig. 23 b) (in a blockwise nature). This is what might be expected
since Fig. 22 d) is the IDM for j2, which is a subscript of x in f.
The pattern now becomes clear and as would be expected, Fig. 23 c)
exhibits: the same structure as Fig. 22 b). Figures 24 a)e) are
representations of 23 a) achieved by varying the index orderings among
the five other possible orderings. The type of behavior described above
is evident in all of these incidence matrices also. Knowledge of this
NM1rot0c K iy NN0 r
to N r0  0 H0 Z 1
x x >. >, ,., ,, >.: K. x ,.
b)
fl1
f21
f12
f22
fs2
O000 011e0
0 N 8 0 N t 0oN N
x t >s >c x: x x x x f e x
X XX
X X iX X I
v X
x x
XX Xxx
S I d
oo00o
21 IX X
f3! XXX
FIGURE 23
EXPANDED INCIDENCE MATRIX
00
f X :X
N ck
c3^":
fli
f21
fS1
f12
f22
f52
eil
g21
312
g22
g52
XX XX XX XXXXXX
XX XX X XX XXX
x %If r %
XXXXXX XXX X
S X X X X
XX XX XX XX XX
X XX X X XXX X
,_,__ X XX Xx___x x x x _x
_ _ii__ i_ _ i i i X . .
o C~ to oM r 0 m w ro cj to Nr K)
0 e J rOr C0 V K o) r o  \ cr N O ro to ro
fi x x xx X Xxx
f12 X X XX XX X X
f21 XX Xx xX XX X
f22 X X XX XX XX XX
f. IX X X X X X X X X X X X
f 2 XX X X X X X XX
gfj XX xx
919 x x x X
gl2 xx xx
c21 xx xx xx xx
g22 XX XX XX XX
g1 XX X XX XX XX XX
if, ix i2
0 N,. > o0 o3I rO 0 to 2
 c Q cq cQ3 c" cj o j to r no ro r ')
fi x xx x x
2 X x X X x X X
gxs x x x xx
f2 XX Y, Xx
91,/ A AX1
t2 XX XX XX XX XX
f221 XX XX XX x x
T2 x x xx xx xx x x
fa xx xx xx x xX
f32 XX XX XX XX XX XX
'32. X X X X X X X X X,
T?3 XX v/ v "^ B, / ? / '* \
g ~ r! ?x I ____ A x^ A ,, A ^ __^. ^ AT.^.
il, if, 2
FIGURE 24
ALTERNATE INCIDENCE MATRICES
o0  c~ MC r0  WM w "z 0 J rl
SX X X XXX
g11 xxx x
Xa xxx x w x x x
S xxx xxxx
9 x x x x
z21 X XX X X X X X
gz x x x X X X X X X
f, xxx x Xxx x X X 
fi xx x x x x x x x x x
X^. xx Xx X
531 X X X X X X X X X X
fx2 X XX X X XX X X XX X
9Q32 XXXX XXXX XXXX
'i l2,i'f
fi
f21
f31
g l
921
g31
f!2
f22
f32
g 12
g22
g 32
i ,if i
FIGURE 24 (cont.)
o_ 0 c\ _r ro
X XXX XXX
x x xxI x xxx
XX XX XXX XXX
X X X X X X X X XX
X X X X
X_ XX XXXXX XXX
x x X X X X xX
XX XX XXX XXX
x x x
XX XX XX >X
xx xx x xx
X X X X X X X X X x
d)
e)
i2viI lif
FIGURE 24 (cont.)
o 0  cJ 0 C C C N NN NN N tO"M tont
 c 3 N M o CJ o to o C< c tor
X X >c >4,X > x >% x >% x >, x >% x
A xx x x X x
xx xx
XXXX XXX x
XXX X XXXX X X X
X X XXXX X X X X
I Xxxx X X X X
XXXX XXXX
X X X XXXXXX X X X
.. x. x x XX x x x x x
fi
f21
g21
f3!
t3i
S}2
f?22
f32
932
behavior, called "Index Imbedding," allows the structure of the full
incidence matrix to be predicted, and in fact visualized from knowledge
of the FunctionVariable Incidence Matrix and the Index Display Matrix.
2.5 Output Set Assignments
A notion central to iterative solution of a set of equations is
that of the output set assignment. Briefly an output set assignment is
the assignment of one variable to each equation such that no variable
is assigned to more than one equation. Variables which are not assigned
are decision variables and have their values defined prior to the
initiation of a solution procedure. There are, however, many output set
assignments for a set of indexed equations which could be derived by
existing algorithms operating on the full incidence matrix which would
be extremely difficult, if not impossible, to solve using doloops. As
an example consider Fig. 25, where the circled elements are the out
puts. There is no reasonable or logical way to write a solution pro
cedure for that output set assignment using doloops.
If, however, we assign outputs to the FunctionVariable Incidence
Matrix and to the Index Display Matrices, a logical, concise solution
procedure for a set of indexed equations can be implemented using do
loops. This means that, whenever a function of a given type is being
solved, the output variable will always be of the type assigned to that
function in the FVIM. Similarly, if a function being solved has index
ii with a given value, the variable which is the output will have its
index value corresponding to ii known. The IDM's which need to have
outputs assigned (called "index outputs") are those which occur in
output set assigned incidences of the FVIM. To illustrate the restric
cj C']  (\J (V (j
7'"* ri. II
xxCi xx.
.I'.
s A / '^ 1 '*" i/;
A 4> (I) A_
IR 25
> 0 a X A
;FI GRE I2
/.* ^ ..' >. A t
~ v \ *' ,^ *
^ ^~\ \'\" \ "f
E 1 .,*.} A A A ; *Vn
FIGUJRE 25
ARBITRARY OUTPUT SET ASSIGNMENT
Fl
21 2
Gii
Gi2
.i ;'
tions on output set assignments consider the system of equations from
Fig. 25, defined structurally by
f(ili2)
g(i ,i2)
x y
JlJ4 J2J3
j274 J193
ji(il): L = L.
U = i
A= 1
j3(i2): i2
j2(il): L = L.
U = U
A 11
j4(i2): i2
U
A
with IDM shown
index outputs
the following
fl1
f12
f21
f22
911
912
921
922
=2
=l
i2+1
in Fig. 26. Assigning x to f
indicated by the circles on the
incidence matrix:
X2
xll X12 X21 X22 Y11 Y12 Y13 Y91
and y
IDM's
( x x x x
( x X X x
x x ) x x
x ( x x x x
x x x (
x x x (
x x X X
Y22 Y23
to g and making the
in Fig. 26, we have
with
ii: L
U
A
=1
= 2
=1
i2: L = 1
I 2
v
J3
L 2 3
x o
12
I1 2
I x x
X
x
J4
I 2
FIGURE 26
INDEX DISPLAY MATRICES
A\
CX
The variables marked with the "D" are decision variables. This is
evident because they have index j3 equal to 1 and, by examining the
index outputs, we see that for j3 the value 1 must be a decision. This
output set assignment easily converts to a doloop type iterative
solution procedure such as:
C READ DECISIONS AND TEARS
READ (5,5000) Y
5000 FORMAT (10F8.0)
10 Do 20 IF=l,2
DO 20 11=1,2
DO 20 12=1,2
IF (IF .EQ. 1) X(I1,I2) = F(X,Y,I1,I2)
IF (IF .EQ. 2) Y(I1,I2+1) = G(X,Y,I1,I2)
20 CONTINUE
C CHECK CONVERGENCE
If y not converged, go to 10
The convergence check would probably check the change in the torn
variables y from iteration to iteration, convergence being defined as
an arbitrarily small change in y. The functions "F" and "G" would be
whatever f and g happen to be.
Although these restrictions on the allowable output set assign
ments have been necessary to facilitate the writing of iterative
solution procedures using doloops, they do not weaken to any extent the
solution procedure deriving power of the algorithms which operate under
them. In fact these restrictions point the way, as it were, to the
types of output set assignments engineers traditionally have employed.
For example, in modeling a distillation column the material balance
equations would either be solved for component liquid or vapor flow
rates, not some liquid and some vapor flow rates. Similarly the flow
rates solved for would always be those leaving the tray over which the
material balance is being made. These are exactly the type of restric
tions which are being made on the FVIM and index outputs. In addition,
because the outputs are made with small matrices (FVIM's and IDM's)
existing algorithms can be used for the most part, which avoids the
necessity of including in a computer program package, such as GENDER,
separate routines for output set assigning indexed and nonindexed
equations. Finally, by assigning the outputs to the small submatrices
much effort is saved when compared to output set assigning the full
incidence matrix.
2.6 Acceleration of Variables
A possible consequence of defining variable indices with a range
is a phenomenon called "acceleration of variables." Acceleration of
variables is said to occur when the number of variable indices in a
problem increases faster than the number of function indices as the
range on a function index is increased. Acceleration of variables
occurs at the lowest level of decomposition, i.e., on function indices
and their associated variable indices. It is, of course, possible to
have the number of variables in an entire problem increase faster than
the number of equations and have no acceleration of variables. Accel
eration of variables is actually a possible result of the manner in
which the indices are defined. To illustrate, consider the following
function index and variable index:
ii: L = 1
U = 2n+l
A = 2
j1: L = L.
U = U.
A=
For n=3 the IDM is
1
3
ii
5
7
1 2 3 4 5 6 7
x
X X x
x x x x x
x x x x x x x
with three more variables than equations. If n is increased to 4,
i.e., if one equation is added, the IDM becomes
j1
1 2 3 4 5 6
1
3
ii 5
7
9
7 8 9
which has four more variables than equations. Acceleration of
variables produces the condition where a solution procedure must be
derived for particular sets of function index ranges, and thus precludes
any universal applicability of solution procedures. For this reason it
is desirable to avoid acceleration of variables. The following
theorem indicates how that can be insured.
x
X x X
x x x x x
x x x x x x x
x x x x x x x x x
X X X X XX X
31
Define the function index range as:
L = M
U = A.(N)(A.M)
A = A.
1
where N and A. are positive constant integers, M a constant integer.
Define the variable index range (it is only with the range that the
problem can occur) as:
L = aL.+(la)i+ki aE{O,l}
U. = bU.+(1b)i+k." bc{O,1}
3 2 3
A = A.
3 2
where k. and k .u are integers and A. is a positive integer, again all
3 3 3
constants.
Theorem: If A. = A. there is no acceleration of variables.
3 1
Proof: Assume A. = A.. The number of functions is:
j 1
A.(N)(A .M)M
Hf A.
1
= N
The number of variables
U. L.
n= max min +1
v A.
U.
U3max =bU,+(1b)imax+k.u
i 3
but imax = U. so,
U = bU +(lb) U+k .
gmax 2 2 g
= U .+k, .
similarly,
L = L .+k .'
rmin J 7
U .+k .(L .+k )
so n +1
A.
U.L. k.uk .
S2 + 3 +1
A. A.
1 2
k .k .
=n f+ 3 3
but since k.u, k and A. are all constants the difference between
number of functions and the number of variables is always the same and
hence there is no acceleration of variables.
There can be no acceleration of variables for a variable index
defined by a list since the difference between the number of functions
(or function index values) and variables (or variable index values) is
always the number of list elements minus one.
CHAPTER 3
SOLUTION PROCEDURES FOR INDEXED EQUATION SETS
The last chapter detailed some restrictions which are placed on
the solution procedure developed. Other restrictions, along with their
motivations, will be presented in this chapter. In order that efficient
algorithms be written which can derive a solution procedure for indexed
equations, it is desirable that theorems be developed which can be used
to eliminate many of the possible choices open to the algorithms. This
chapter presents many such useful theorems, along with comments on their
application to problems. Most of the properties are related to a
single example problem presented in the beginning of the chapter.
3.1 Definitions
Before proceeding with the example it is helpful to define some
of the terms used throughout the chapter.
3.1.1 FVIM Outputs
As discussed in the last chapter, the FVIM must be output set
assigned. The variable types chosen as outputs are called the "FVIM
outputs."
3.1.2 FVIM Decisions
In general, the number of variable types in an FVIM will be
greater than the number of function types. When this is the case, some
of the variable types cannot be assigned as outputs. These variable
types must then be decision variables. The variable types which are
decisions are called the "FVIM decisions."
3.1.3 Index Outputs
The concept of index outputs was also introduced in Chapter 2.
When the IDM is treated as an incidence matrix and output set assigned,
the assignments made are called the "index outputs."
3.1.4 Index Decisions
Just as the FVIM may have both outputs and decisions, so may the
IDM. If any columns in an IDM remain unassigned after the index outputs
are chosen those columns are the "index decisions."
3.1.5 Imbedded Loops
The concept of index imbedding was introduced in the last chapter.
The natural way to achieve index imbedding in a FORTRAN program is to
define each index in a "doloop" and nest the doloops. Figure 31
illustrates a typical nesting of doloops. The innermost loop, loop 3,
is executed seven times for each pass through the next innermost loop,
loop 2. Similarly, loop 2 is executed five times for each pass through
loop 1. The index imbedding represented by this nesting of doloops
would be 13 imbedded in 12 which is in turn imbedded in Il. In terms of
incidence matrices, each 12 incidence would actually represent an 13
IDM. The nested doloops are termed "imbedded loops." The indices are
said to be "nested inside" or "nested outside" the other indices.
3.1.6 Decomposed Problem and Expanded Problem
As was seen in the last chapter, the FVIM and IDM's contain all of
the structural information necessary to describe a set of indexed
equations. The problem representation in terms of the FVIM and IDM's
is called the "decomposed problem." When represented in terms of the
1 /
30
20
10
Q. I
D
D00(
0
AT= 1210, 1
12= I,5,
I3= I, 7,v
0 CONTINUE
20 CONTINUE
30 CNNTINUEU
loop I
loop 2
loop 3
end loop 5
end loop 2
end loop I
FIGURE 31
NESTED DOLOOPS
full incidence matrix it is called the "expanded problem."
3.1.7 General Solution Procedure
Existing algorithms generate a solution procedure for a particular
set of equations. It is the aim of the algorithms here, not only to
treat indexed equation sets, but to generate solution procedures which
are independent of the index limits. A solution procedure which
satisfies this requirement is called a "general solution procedure."
3.1.8 Function Ordering
There is a natural ordering of function indices. They are either
incremented or decremented in steps of equal size. There is no such
natural ordering for the function types. In fact, the order of the
functions is usually determined from precedence ordering considerations.
The order in which the function types appear in the FVIM is called the
"function ordering." The function ordering can change as the solution
procedure is generated.
3.2 An Example Problem
Suppose that a solution procedure is desired for this set of
equations:
f(il,i2) = 0 = f[x(j1,j3),y(jl,ji),z(jl)]
31
g(il,i2) = 0 = g[x(j2,j4),y(j2,j4),z(j2)]
The decomposed representation of these equations is shown in Fig. 3Z,
which also defines the indices shown in the equations. To solve these
equations, which may be assumed to be nonlinear, a solution procedure
must be specified.
Suppose the following solution procedure is chosen. Assign x to f
and y to g. Assign the index outputs as shown in the IDM's. Both
FV M
f ( i,i, )
Index
1, L=
Definitions
jr: j1 l
II
A= I
i2: L=I
U=2
j : L= Li
U= Ui2
A=I
SD s
'2
4
i27,
," w C, !
"i"
; v/ 1!
FIGURE 32
DECOMPOSED PROBLEM
J2: ii
i,+i
J4: i2
j3
2 x x
indices will be incremented and the function ordering will be as it is
in the FVII. The index nesting will be il,if(function type),i2. The
expanded incidence matrix resulting from these choices is shown in
Fig. 33. It so happens that the equations fully precedence order for
the solution procedure chosen.
This example problem illustrates many things described in Chapter
2. Assigning x to f and y to g is an example of an FVIM output set
assignment. Since z was not assigned it was an FVIM decision variable.
Similarly the IDM's contain index outputs and index decisions. All
variables with index ji were decisions for ji equal zero, and those with
index j2 were decisions for j2 equal one. The problem decoupled in the
variable type y for index i1. This example will be referred to by other
sections in this chapter because it illustrates so many valuable
concepts.
3.3 Output Sets and Decision Variables
The choice of an output set and decision variable set almost
entirely defines a solution procedure. This section discusses these
choices as applied both to function and variable types and to indices.
3.3.1 FunctionVariable Output Sets
In the example problem a variable type was assigned to a function
type. When this is the case, implementing the solution procedure is.
simplified by the fact that each function type only has to be solved for
one variable type. This is particularly helpful in computer implemented
solution procedures, such as those to be generated by GENIE. For
purposes of deriving solution procedures for sets of index equations,
the requirement is made that outputs be assigned such that all outputs
DDD DD
I I, I
D DD
 CIA" cQ CM Cki
0 0 0 C0 0 NC% C4 (C C\J C I fO rC O tO' 0 t~ C S
x x x xlx) x x
x xxxx xxx
xx xxx x @ x
Xx x x X X
x __x x x _____
xxxx @xYx
X X XX XXXX
_x xx xx x
X VX ," V ^7
A ^f\,.f) /A /
,, V XX v /^? v
I. L~ Y
1,
D = DECISION
VARIABLE
FIGURE 33
EXPANDED INCIDENCE MATRIX
f l
f12
gTI
g2
f22
9
g22
T32
f32
gs32
of a particular function type are of the same variable type. This
eliminates the need for more than one representation of a function type
in a computer program. It also avoids having to determine during
solution which output variable type would be required for a particular
set of function index values.
This restriction is a reasonable one to make, and actually is
widely used in accepted solution procedures for typical modeling
problems in chemical engineering. One shortcoming might appear to be
the inability to assign different outputs to the equations describing
different sections of a process, such as a distillation column. The
rectifying and stripping sections are often treated separately and this
restriction would seem to prevent that. This is not the case, however,
since it is possible to describe the rectifying and stripping sections
as separate problems with the appropriate connection equations. This
has the added advantage that different solution procedures can be
derived for the two sections.
It should be noted that each output variable in the example
problem had the same number of indices as the function for which it was
the output. Clearly, z could not have been chosen as the output for
either f or g since there are six equations of each type, but only five
variables of type z. If the number of indices on a function type
exceeds the number of indices on a variable type which occurs in that
function, the total number of functions can be made to exceed the total
number of variables by choosing appropriate index limits. For this
reason, a variable type cannot be assigned as the output for a function
type which has fewer indices than the function type.
If the number of function indices were less than the number of
variable indices for an output selection there could be many more
variables than functions. While an output set assignment would be
possible, determining the output variable index values for a particular
set of function index values would be difficult or impossible. For this
reason it will be required that the output variable type has the same
number of indices as the function type for which it is the output.
Suppose next that each of the variable indices for an output
variable type does not depend upon a different function index. Then
there must be a function index which does not appear in the output
variable type. The range on this index could be increased until there
are more functions than variables. This would mean that the variable
type could not be the output. Therefore, each of the variable indices
must depend upon a different function index.
The following restrictions, then, are imposed on the FVIM output
set assignment:
1. A separate variable type must be the output for each function type.
2. The number of indices for an output variable type must be the same as
for the function type for which it is the output.
3. For an output variable type, each variable index must depend upon a
different function index of the function type for which the variable
type is the output.
An interesting observation about the restrictions imposed on the
output choice is that it is possible to have a set of equations which in
expanded form has an output set assignment, and yet not be able to
assign outputs to the FVIM. As an example, consider the constraint
equation for the sum of the mole fractions in the liquid phase on a
stage in a distillation column.
NC
LM. 0 = x..l 32
j=1
where: LM = Function Type
i = Stage number
j = Component number
NC = Number of components
x.. = Liquid mole fraction of component j on stage i
The function type, LM, has one index. The only variable type is x,
which has two indices. By the restrictions stated above, there is no
possible output set assignment. There are two possible solutions to
this problem. The first is to rearrange algebraically the equations in
an attempt to eliminate the problem. The second is to assign, from the
other variable types in the problem, a different variable type, with the
correct number of indices, as an "implicit output." This variable would
have to be connected to the function through the other function types
(which means that through algebraic substitutions the implicit output
could be made to appear in the function, although this is not done).
A technique such as NewtonRaphson is then employed to converge the
implicit output.
3.3.2 Variable Type Decisions
In the example problem the variable type z was a decision variable
type. The only restriction on the choice of decision variable types is
that they be made so as not to reduce the number of variables below the
number of equations. If the FVIM outputs are chosen before the
decisions there can be no problem. If, however, any decision variable
types are to be chosen before the FVIM outputs are assigned care must be
taken to insure that the FVIM remains output set assignable. For the
example problem, had x or y been chosen as the decision variable type,
an output set assignment could not have been made. FVIM decisions may
not be chosen so that the restrictions on FVIM outputs must be violated
in assigning outputs.
Consider the FVIM in Fig. 34 a). There are three function
types and four variable types. Choosing x as the decision variable
type results in the FVIM in Fig. 34 b) and choosing z results in
Fig. 34 c). The FVIM in b) is not output set assignable according to
the rules whereas the one in c) is. The decisions must be chosen so
that there are the same number of variable type as function types
with any given number of indices.
3.3.3 Index Outputs
Index outputs are the outputs assigned in the Index Display
Matrices. They indicate which of the variables of the output variable
type will be the output for particular function index values. In
extending index outputs from the IDM to a general solution procedure
knowledge of the following properties is necessary.
Theorem 31: All Index Display Matrices, when treated as incidence
matrices are output set assignable provided they are not null.
Proof: The proof is by induction and is divided into two parts, one for
a list and one for a range.
List: There must be at least one list element, which is an offset from
the function index value. This is known from the definition of variable
index lists. Therefore, for one function index value there is an output
set assignment. Now suppose that there are output set assignments for
an IDM with k rows. Row k+l introduces a new function index value,
1
a)
f(il ,i2) J, J4 J2 J4 j2 i3
g (1,12) Ja J5 j I 4 J
h(i ) j j5 ij j2
= .
b)
f(igi2)
h(i ,)
c)
.'(i
g(i,
,i2)
,'2)
h(i, )
i J4 i2 I
j3 j5 ij
3 1 2
w x y
JI J4 J2 J4 J2
j3 j5 iJ 1 4
2 i5 j3
FIGURE 34
FVIM DECISION CHOICES
larger than any of the other function index values. The largest offset
in the list, when added to the new index value, must produce a new
variable index value which can be assigned as an output. The first part
is proved.
Range: Define an index range as follows:
i: L = L. j: L = a(i)+(la)L.+kz a{O0,1}
U = U. U = b(i)+(1b)U.+k bk{0,1}
2 1 U
A = A. A = A.
2 1
where k and k are offsets. Note that L. cannot be greater than U. as
this could result in a null IDM (for k=kU for example).
Assume L. = U.
i 2
then L. = k
and U, = k
In order that the IDM be nonnull k must be greater than or equal to
U
k For one row the IDM is output set assignable since it is nonnull.
Now assume that the IDM is output set assignable for k rows. The number
of rows can be expressed as:
U.L
r i
k A.
and the number of columns as
b(imax)+(lb)U.+k a(imin)(la)L k~
uc A + 1 33
but imadx = U,
and imin = L.
so eqn. 33 becomes
Ui L +k ki
k u + 1
a.
k k.
k A.
Now suppose that a new row is added. The number of columns becomes:
U.+A.L.+k k
Ck+1 + 1 34
1
= k +
so there is a new column which can be assigned as the output for the new
row. Theorem 31 is proved.
Corollary: It is always possible to specify the index outputs as
offsets from the index values in an IDM.
Proof: For the list choose any list element as the index output offset.
For the range note that for L.=U. the output can be expressed as an
7 1
offset from the function index value, say k Now assume that for
L.
3 1 U
U. is increased by A.. A new function index, the new U.,
1 1
has been introduced. At the same time a new variable index value,
U. +k not previously assigned as an output is introduced, which is
offset from U. by the same factor as all previous outputs. The
corollary is proved.
The theorem and its corollary not only guarantee that index
outputs are possible for arbitrary limits, but also that there is an
index output set assignment which makes the determination of output
variable index values easily obtainable from the function index values.
3.3.4 Index Decisions
In an IDM the columns not assigned as index outputs are termed
index decisions. Index decisions result from the lower limit offset
for a variable index range being less than the upper limit offset, or
from there being more than one list element for the variable index list.
Theorem 32: For a general solution procedure, any index decisions must
be declared as offsets from the function index upper and lower limits.
Proof: Since for a general solution procedure the index range is arbi
trary, the only function index values which can be guaranteed to occur
are L. and U.. Thus the only variable index values which can be
7 1
guaranteed to exist are offsets from L. and U.. (Remember that an
1 1
offset can be equal to zero.) The index decisions must therefore be
chosen as offsets from L. and U.. Theorem 32 is proved.
1 2
To see the effect of index decisions on a problem recall the
example problem at the beginning of this chapter. For j, the value zero
(actually L. 1) was an index decision. All variables with index j,
11
equal to zero were decision variables. Similarly, all variables with
index j2 equal to one were decision variables.
3.4 Index Imbedding
As was seen in Chapter 2, the order of nesting chosen for the
indices has a profound effect on the appearance of the expanded
incidence matrix. Also, as in the case of the example problem at the
beginning of this chapter, it will affect the efficiency of the solution
procedure.
First, index ordering (the order of nesting) can affect the
number of function evaluations necessary to reach convergence. To see
this, consider the incidence matrices c) and d) in Fig. 35, which are
generated from the Index Display Matrices a) and b) with different index
orderings. Suppose that these incidence matrices represent a set of
equations to be solved for the outputs indicated by the circles. To
FIGURE 35
DIFFERENT INDEX ORDERINGS FOR INDEX IMBEDDING
6ii
Iv
Ex xxB
L\ {4
('A ^l* 'A*ut* ^^'w<^^^
t: S.? V
6 (1S/ A 5
bs) ,
i~XXX
,.
c~caum0a
3.
0,
V\`
'A, A
X
)r '< i
A' A,
'\ J
S',: ".
,..
X
X XX
IkX XX
7 ,
Ki
'A. 'A I 'A 'A 'A
.', 5. ', s
s s
'A A
isi
\f 'V v' '1
' 1
' X( j
/* S 59 
 s.
x'x x I. x
i
6 f '
I i '^ ."t Ii t r
S'V 'P v V .
* J \.
a '
49
Cm w. wVZXIllW YU rm
r_^ ~ Y
t
d)
solve for the (1,1) variable in matrix c) the (1,5) and (1,9) variables
must be torn and converged. To do this involves evaluating all of the
functions, since the order of the rows in the incidence matrix is the
order of solving the equations. To solve for the same variable in the
second matrix involves evaluating only three functions until conver
gence. Since the equations in both incidence matrices are the same,
the extra function evaluations in the first method would be wasted. The
same number of iterations would be required to converge the first
variable for both methods, but the number of function evaluations would
not be the same.
Theorem 33: Let il be a function index whose IDM does not fully
precedence order and i2 be a function index whose IDM does fully
precedence order. If the two indices are adjacent in index ordering,
the number of function evaluations necessary for convergence of the set
of functions they describe for the ordering i2il will be less than for
the ordering ili2.
Proof: Convergence is necessary for blocks of variables with i2
constant, iI ranging over its allowable values. Call the set of
variables in a block with a particular value of i2 an is block. Let
the number of iterations required to converge each i2 block be k, .
'2
For the ordering i2il the number of function evaluations to converge
the entire set of variables would be
M2 M2
NI = (k i2. ) 1= M1 k. 35
i2=1 i= =1
Mi and M2 are the number of rows (and columns) in the IDM's for i1 and
i2 respectively. Equation 35 states that the total number of function
evaluations is the sum of the number of iterations for each i2 block
times the number of functions in each i, block, over all such blocks.
Now suppose that the order of the indices is ili2. The solution
procedure would converge the first variable in the first i, block
simultaneously with the first variable in each ii block. This is
equivalent to converging all variables in the first i2 block. The
number of function evaluations would be
N2 = kiM1M2 36
since all functions must be evaluated. Equation 36 represents only the
number of function evaluations required to converge the first variable
in each iI block. The total number of function evaluations required
would be
M2
N = MM, I k. 37
T i 2
i2=1
N is greater than N1 by a factor of M2. Since the IDM for i2 does not
fully precedence order, M2 must be greater than one. The theorem is
proved.
This theorem indicates the desirability of nesting indices whose
IDM's fully precedence order outside those whose IDM's do not. A
problem arises when a function index is to be nested outside the index
for function type.
Consider the incidence matrices in Fig. 36 a) and b). Any of the
four partitions in the incidence matrix a) themselves have legal output
set assignments. In the incidence matrix b), however, the (1,0) and
(2,1) blocks do not have legal output set assignments. This condition
is caused by ii being nested outside i The function g contains only
variable index j2, and hence no variables with index value equal to
zero. The result is that if the index output is to be assigned inde
pendent of function type it must be chosen from a special IDM. This IDM
Nt?
0 * C jO O cj fl)
1 / /
SV v"
 *"* ''' di
) } : t, n
I si
p 1
o c)  C~ Cxl iro to
A. s*** .^,; >' ^ >'
t  ) :,*.;, ixiE.
mea r eaev *mm "I e.e
f. '' 4
I I. I
H fi i>
*j :* 'f
Fj I
h F .. Li V
i, t
02~
tt~~' '"~
FIGURE 36
INDEX DISPLAY MATRICES
is the result of performing a logical "AND" operation on all IDM's which
occur in incidences of the FVIM which are assigned as outputs. For
example, refer to Fig. 36. If y is assigned to f and x to g, either
ii or i1+i can be the index output. If, however, x is assigned to f
and y to g, only ii can be the index output. The logical "AND" then
provides a means for discovering which index outputs can be chosen inde
pendent of function type. When only structural considerations are being
made, this can save time. This is because the logical "AND" can be
performed faster than an output set can be assigned and because only one
output set assignment is required in this case.
An important result is that index ordering does not affect the
number of tears in the problem. Consider the two different index
orderings in Fig. 37. Both incidence matrices can be solved with four
tears.
Theorem 34: If two function indices are adjacent in ordering, the
number of tears for the expanded incidence matrix is independent of the
ordering of the indices.
Proof: Let the IDM's for the two indices be called a and b. Let the
number of tears for a be NT and for b be N ,. Let Nca and Ncb be the
number of columns (and rows) in the IDM's a and b respectively. Assume
b is nested inside a. Then the total number of tears is
N = NTa x No. of columns for each occurence of a
+ NTb x No. of occurences of b not already torn
= Ta x Nb + NTb(ca NTa
= NTa x N b + N Tb x Nca NTbNTa 3
which is symmetric in NTa, N Tb NTc and Ncb. Thus, the number of tears
is independent of the index ordering. The theorem is proved.
a) i
T>* t f ^
X Xi
V
S v B
r,
"I
'.^a J
T T T T
4 o
XX / "
x x x x
xxx x x xY
X X X X NX
X X X X v !
Si,
GX XEX X3
I
x x":x x:
j ^_^ x x'
FIGURE 37
TWO DIFFERENT INDEX ORDERINGS
b) tJ
L [ /u
',7"'7:'!~
As shown earlier, a function index may have several different
choices for the index output. It has also been shown that an index
which fully precedence orders should be nested outside those which do
not. If a function index has only one possible index output, it must
fully precedence order. Since all IDM's are output set assignable, for
a function index to fully precedence order, with only one index output,
all of its IDM's must have the same index output and that index output
must be the only possible one for each of the IDM's. This is equivalent
to saying that the logical "OR" of all IDM's for a function index
produces an IDM which fully precedence orders. This, then, provides
an easy way to discover function indices which will give a full
precedence ordering of blocks.
It has been stated that function indices will either be incre
mented or decremented. Some variable indices are best treated by
incrementing and some by decrementing. For example, a lower triangular
IDM would fully precedence order for an incremented index and an upper
triangular IDi would fully precedence order for a decremented index.
Therefore, for some functions it may be desirable to increment a
particular index and for others to decrement it. If the index in
question is nested inside function type a simple test of function type
can be made to determine which of the two methods is to be used. If,
however, the function index is nested outside function type, no such
test is possible and the index must be either incremented or decre
mented.
3.5 Decoupling
Freedom in the choice of index decisions often affords the
opportunity to choose a solution procedure much simpler than would be
expected from an examination of the decomposed problem. Consider the
decomposed problem in Fig. 38. At first glance it would appear that in
order to solve for either x or y in f, the variable type not chosen
as the output must be torn. This is not the case. By assigning y to f
and x to g, and by making the index output assignments ii and i2+1, the
incidence matrix in Fig. 39 results. For the outputs shown, the first
nine columns represent decision variables. The first block of variables
of output type y can be calculated without tearing any variables of type
x. In fact, the entire set of equations can be solved without tearing
any variables of type x. This phenomenon is called decouplingg."
Decoupling is made possible by a judicious choice of index outputs.
This problem is said to have decoupled in the variable x. Decoupling
of a function type f in its output variable type x occurs when the
function ordering would appear to require that the variable x be torn to
solve functions not of type f, but the variable type x actually does not
need to be torn for that reason.
Theorem 35: A variable type x will decouple a function type f, for
which it is the output, if and only if the following conditions exist:
1) A variable index for x must have index output offsets strictly
greater than (or less than) the possible index output offsets for all
other incidences of that variable in any other function type. 2) The
index in 1) must be nested outside the index for function type.
Proof: The proof will be presented for the "greater than" case. The
proof for the "less than" case is analogous.
If part The proof will be by induction. Assume 1) and 2) hold
for some set of indexed equations. Let the index in question be called
x y
f(igi2) i j2 j4 ij
g(il'i2 )Jl J3 J, J3
i i I
x
%x x x
j2
X X
2 r x x
j4
Tx x xB
IXXX
%x X x~
13
x x
2x xX
FIGURE 38
A DECOMPOSED PROBLEM
fl
911
921
931
g,
f12
f32
912
922
932
EXAMPLE OF DECOUPLING
0 0 0  J J c N NCJ CUN N re o tO M
 cr cn N N I
x< X x X X X X X >< x
x x x x
xx x A xx
x x xxx !xx
Sx xx x x x x x xx
Sx x R x
XX X XXX
Y x x x xx x x o xx
xX X X X X x (X
[ xx x xxxx
y v w ,,,'"
x x x x Ix x x x( x x xI
F E 3if 1
FIGURE 39
ii. For il=L. all incidences of x in functions other than f which have
11
an index depending on ii, have index values less than the lowest valued
index output (by condition 1), and thus are decisions. Since 1i is
nested outside function type (by condition 2), the first pass for il
can be completed without tearing x for functions other than f. Now
assume that k passes through the i1 loop have been completed without
tearing any x's for functions other than f. The index outputs for il
in x, for iteration k, were greater than any of the other index offsets
for ii in x. Thus these index outputs of x are greater than or equal to
any of the index values for x in functions other than f for iteration
k+l. The functions other than f, then, do not require any tear of x
for iteration k+l. The if part is proved.
Only if part The proof will be by contradiction. Assume 1)
does not hold. The index output offset for x in f is less than or
equal to some other index value for il in an incidence of x in a
function, other than f, which must be solved before
f. The variable x in this function must be torn for iteration k since
there are index values of x present which will be outputs for iteration
k. This contradicts the definition of decoupling, so condition 1)
is required for decoupling. Assume 2) does not hold. Then a function
not f, which contains the variable x must be solved for all values of
ii before any function f is solved. This would require tearing x which
would contradict the definition of decoupling, thus condition 2) is
required. The theorem is proved.
This theorem provides the basis for algorithms which discover
conditions which allow decoupling.
3.6 Blocking Factors
The blocking factor defines the size of the Index Display Matrices
for the various function indices. It is desired to perform analyses on
the IDM's and have the results of those analyses be applicable to an
expanded incidence matrix with arbitrary index limits.
It is possible, by employing different blocking factors, to
analyze an IDM in the same manner twice and reach different results.
For example, the blocking factor can affect the ratio between the number
of tears and the number of rows in a problem. (See Fig. 310.) In this
figure, a) and b) are Index Display Matrices representing a tridiagonal
matrix with first and last columns declared as decisions and omitted
from the figure. For a blocking factor of 2, either of the two possible
output set assignments results in one tear for the two rows. When the
blocking factor is increased from 2 to 3, two of three possible output
set assignments only require one tear for the three rows. The require
ment that the order of solution be dictated by the IDM (i.e., be top to
bottom or bottom to top) is relaxed for the analysis within the blocks.
When the problem is expanded, the IDM for an index is made up of many
blocks. These blocks must be solved in either ascending or descending
order.
It is not always the case that the solution procedure derived
for a block will hold for the expanded IDM. Consider the two expanded
IDM's in Fig. 311, each representing one of the onetear, blocking
factor of 3, solution procedures from Fig. 310. The first IDM reflects
a solution procedure which still allows for solution of the equations by
solving the (1,1) block then the (2,2) block. The second, however, does
not allow for solving either block before the other. When a solution
b)
X X
X X
T T
vR
l x
LX@X
xx
x8
T = TEAR VARIABLE
FIGURE 310
ALTERNATE BLOCKING FACTORS
a)
X X
XX
T
xi
IxI
T
x x
xx
X 6)
A A l i
!
xx
xx
xl x
x x
x
FIGURE 311
EXPANDED INDEX DISPLAY MATRICES
I
i
x
procedure for a blocking factor remains effective as the IDM expands,
the solution procedure is said to "extend."
Theorem 36: Call a block with blocking factor equal to k a kgroup.
For the solution of kgroups in ascending order the following condition
is necessary and sufficient for extension of a solution procedure. The
decision variables offset from U. must all be k greater than tear
variables within the last kgroup. An analogous condition exists for
solution of kgroups in descending order.
Proof: The proof for ascending order only will be presented. Analogous
arguments hold for the case of descending order.
Sufficient: Suppose that all decision variables offset from U. are
offset from the kgroup tear variables by k. When the next kgroup is
added, the decisions offset from U. are no longer decision variables.
1
In order to solve the next to last kgroup, values for these variables
must be available. The values are available since these variables have
become tear variables. Thus if a solution procedure holds for n
kgroups it holds for n+l kgroups. It must hold for one kgroup,
otherwise it would not be a solution procedure. The condition is
sufficient by induction.
Necessary: Assume that some decision variable offset from U. is not
k greater than a tear variable in the last kgroup. When the next
kgroup is added, this variable's value must be available to solve the
next to last kgroup. It is no longer a decision and is not a tear, but
is an output of a function in the last kgroup and hence its value is
not yet calculated. For its value to be known it must be a decision or
tear variable, but this is a contradiction. Thus the solution procedure
does not extend and the necessary condition holds by contradiction. The
theorem is proved.
This theorem provides an easy method for determining which
solution procedures to IDM's are worth analyzing and which are not. If
there are no decision variables (decision indices actually) involved,
the analysis is even simpler. For a variable index defined by a list
with no decisions, letting L.=U. proves there is only one list element.
Thus for the index there is a full precedence ordering. For a variable
index defined by a range, four cases must be considered. First consider
L. and U. offset from L. and U. by the same factor (to result in no
3 3 1 1
decisions). In this case a solution procedure does not extend since,
for any blocking factor k, k1 tears are made in the first block and k
tears must be made in the second block. Next consider L. offset from L.
7
and U. offset from i, again by the same factor. A full precedence
3
ordering exists as this results in a lower triangular IDM. Similarly L
offset from i and U. offset from U. results in full precedence order.
7 .7
With both limits offset from i by the same amount the IDM is diagonal
and fully precedence orders.
Knowledge of the behavior of solution procedures from index
definitions is very useful in choosing among the various solution
procedures. Many possible solution procedures can be rejected without
extensive analysis because they are known to be inefficient for the
problem at hand.
CHAPTER 4
ALGORITHMS FOR DERIVING SOLUTION PROCEDURES
FOR SETS OF INDEXED EQUATIONS
The algorithms implemented in the GENDER system for deriving
solution procedures are capable of treating nonindexed equations only.
In deriving a solution procedure GENDER performs the following analyses.
First the acyclic algorithm is applied to determine which parts of the
problem are acyclic and which contain recycle calculations. An output
set assignment is then made. After the outputs are assigned the
precedence ordering is found and the groups of functions and variables
requiring simultaneous solution are identified. The last step is to
choose the tear variables for each group.
The treatment of indexed equations by GENIE is similar in approach
to that of GENDER. First, the FVIM is analyzed by the acyclic algorithm
to discover the acyclic nature of the set of equations. The decoupling
algorithm is then applied to the set of function and variable types
involved in recycle calculations in an attempt to identify further
acyclic structure. Outputs are next assigned to both function types
and function indices. A precedence ordering of function types is
determined and groups are identified. Finally, the tear variables
within each group are determined by a minimum weighted tearing algo
rithm.
The algorithms presented in this chapter perform the particular
analyses necessary for deriving a solution procedure for indexed
equations. The manner in which these algorithms are used is explained
in section 5.3. Some of the algorithms are new, while some are adap
tations of algorithms previously developed. Where existing algorithms
are used, the alterations to the problems necessary to allow their use
are detailed.
4.1 Decoupling
Decoupling can result in significant simplification of the solu
tion procedure for a set of indexed equations. Since the presence of
decoupling depends upon the functionvariable and index output set
assignments made, it is desirable to choose these in such a way that
decoupling will result. In addition, the function ordering also
affects the degree to which decoupling is possible. Finally, the
direction in which the function indices are incremented (ascending or
descending), influences any decoupling. For a particular direction of
index incrementing, some IDM's may fully precedence order while others
may not. All IDM's for a function index must fully precedence order
for it to be a candidate for a decoupling index. Recall that decoupling
results in a full precedence ordering of function types. This is
equivalent to a full precedence ordering for an index, the index in this
case being function type. It has been shown that nesting a fully
precedence ordered index inside one which does not fully precedence
order results in increased function evaluations. This would be the
case if a decoupling index did not fully precedence order. What is
required, then, is an algorithm which will analyze the FVIM and IDM's
to choose the function ordering, function outputs, index outputs, and
direction of index incrementing in such a way as to decouple the
functions if possible. The following algorithm does that. (An example
which follows illustrates the algorithm.)
I. Produce the FVIM. Flag all function indices whose IDM's
do not precedence order.
II. A. Break all assignments made by this algorithm.
B. Choose a set of directions for incrementing the function
indices not previously analyzed. If there are none,
go to III.C.
C. Untag all function and variable types and function
indices.
III. A. If the set of unassigned function and variable types does
not fully precedence order, go to IV.A., otherwise,
continue.
B. A full precedence ordering has been discovered. Stop.
C. Decoupling will not produce a full precedence ordering.
Stop.
IV. A. Choose an untagged, unassigned variable type. If there
are none, go to II.A. Tag the variable type. Untag all
function indices and unassigned function types.
B. Choose an untagged, unassigned function type. If there
are none, go to IV.A. Tag the function type.
C. If the variable type chosen is not a legal output of the
function type chosen, go to IV.B.
V. A. Choose the first untagged function index. If there is
none, go to IV.A. Otherwise, go to VI.A.
B. Choose the next untagged function index. If there is
none, go to IV.B.
VI. A. Does the variable index depending upon the chosen func
tion index decouple the function type from the unassigned
functions? Tag the function indices which fail for all
functions. If no decoupling, go to V.B.
B. Assign the chosen variable type to the chosen function
type. Give the function type a position in the function
ordering immediately after all of the presently
unassigned functions.
C. Make the index output set assignments which produced
the decoupling. Flag the function index as a decoupling
index. Go to II.C.
The decoupling algorithm performs an exhaustive search over all
combinations of directions of index incrementing. Additionally, when
ever an output assignment is made, all remaining function and variable
types are untagged and the algorithm restarted on the remaining
unassigned function and variable types. Generally, an exhaustive search
of this type would be extremely long. In the case of the decoupling
algorithm, the search will be relatively short since the FVIM will
usually have only a few function and variable types and there are
usually only one or two function indices.
For some of the function index directions it may be possible to
reject certain variable types at the outset. Remember that the aim of
decoupling is to manipulate the choice of solution procedure so that
variables which would otherwise have to be torn do not have to be torn.
Suppose that the function index ii is incremented in ascending order.
Suppose, also, that the variable type x is being examined to determine
if it can be used to decouple a function type from the others. If x has
any FVIM incidence, say in function f, with a variable index defined by
a range with upper limit offset from U. decoupling will not hold for
that function index. The reason is that, for iz=Li the first pass
through the "ii loop," the function f will require values of x for all
values of ii. This being the case these values must be supplied by
tearing x, hence decoupling could not hold. An analogous condition
holds for the function index being decremented rather than incremented.
Another "shortcut" is to note that, if for some variable type all
variable indices depending upon a particular function index are the
same, decoupling cannot occur for that function index.
Figure 41 presents an example problem to which the decoupling
algorithm will be applied. The first step is to produce the FVIM, which
is at the top of the figure. No assignments need be broken as none have
been made. Choose ascending order for both indices initially. The
function and variable types are untagged. Now the algorithm begins in
earnest. The FVIM does not fully precedence order, so the algorithm
skips to step IV.A. The following steps are then carried out:
IV. A. Choose x. Tag x.
IV. B. Choose f. Tag f.
IV. C. Legal output; continue.
V. A. Choose ii.
VI. A. No decoupling (all variable indices are j1). Tag i
V. B. Choose i2.
VI. A. No decoupling (U. =U. ). Tag i .
34 12
V. B. No more function indices.
IV. B. Choose g. Tag g.
IV. C. Legal output; continue.
1
f(iI ,i2)
g(i, ,i2)
ig
L=
U=3
A=
U2
L=I
U=2
A=I
i
ii
j3
i2
J2
L= Li
L =Li,
A= I
FIGURE 41
EXAMPLE PROBLEM FOR DECOUPLING
J J4 jz J3
i J j3 j4
V. A. No untagged indices.
IV. A. Choose y. Tag y. Untag f and g. Untag i1 and i2.
IV. B. Choose f. Tag f.
IV. C. Legal output; continue.
V. A. Choose ii.
VI. A. No decoupling. (No possible index output for j2 greater
than all index values for ji.)
V. B. Choose i2.
VI. A. No decoupling (fails for all functions since U j=U ).
Tag i2.
IV. B. Choose g. Tag g.
IV. C. Legal output; continue.
V. A. Choose ii.
VI. A. Decoupling.
VI. B. Assign y to g. Solve g after f.
VI. C. Assign iii+ as index output for i1 in g.
II. C. Untag everything.
III. A. There is a full precedence ordering of unassigned function
and variable types (that is, after column For y and row
for g deleted).
III. B. Stop.
The expanded incidence matrix for this example, arranged to
reflect the decoupling, is shown in Fig. 42. The index output for ii
in f is chosen to be i1+1 and for i2 the index output is chosen to be
i2, all choices being arbitrary. There are only four tear variables in
this problem, and the convergence loops (each 2x2) are the smallest
possible.
DECISIONS
x x x xI I ,
,x X jx
x x x xx
x
r V
xx xxx x
x x x
I x x x
I xxx x@q
FIGURE 42
A DECOUPLED PROBLEM
fi
Ti2
ci
g2
C
fa
2g,
931
g32
Although full decoupling may not be possible for a set of indexed
equations, a partial decoupling may be. A partial decoupling results
when some function types are decoupled from the set of equations, but
not enough to render the FVIM acyclic. As a result a subset of the
original function and variable types must be solved simultaneously. If
full decoupling is not possible, there will, in general, be a partial
decoupling for each set of index incrementing directions. Since a
partial decoupling results in a smaller set of equations which must be
solved simultaneously than does no decoupling, it would be desirable to
choose one of the partial decoupling schemes for the solution procedure.
It is possible to characterize the "degree of decoupling" by the
expected number of functions decoupled (where the expected number for a
function is calculated from its expected index ranges, not the blocking
factor). The different partial decoupling schemes could then be
compared quantitatively and the one exhibiting the greatest degree of
decoupling chosen for the solution procedure. This is the strategy
adopted by GENIE.
4.2 Index Output Set Assignments
Index outputs are assigned by treating the IDM's for each variable
index as an incidence matrix. The IDM's are output set assigned by
conventional output set assignment methods. The resulting output set
assignment is then available for defining the index outputs, should that
be necessary.
As indicated in Chapter 3, it is sometimes desirable to assign
index outputs to the IDM resulting from a logical "OR" or a logical
"AND" being performed over all IDM's for a function index. Both of
these special IDM's are produced by the same method. An IDM for the
function index in question is chosen. The logical "AND" or logical "OR"
of each of the other IDM's (for that function index) and the chosen IDM
is then stored in place of the chosen IDM. By successively operating on
all of the IDM's required, the end result is the logical "AND" or "OR"
of all desired IDM's. The resulting IDM is then available for analysis.
In the case of the logical "AND," if it can be output set assigned, the
result is an output set assignment which can be used for all IDM's
depending on that function index. In the case of the logical "OR," if
an output set can be found to fully precedence order, then the function
index fully precedence orders. When this is the case this function
index should be nested outside those which do not fully precedence
order.
The output set assignments are made by an algorithm which allows
the incidence matrix elements to have weights and then minimizes (or
maximizes) the sum of the weights of the output elements (Gupta et al.,
1974). In order to fully utilize that capability, weights should be
assigned to the IDM elements. If this is not done an arbitrary output
set is chosen. The weights should be assigned to the elements in a
manner which indicates to the algorithm which are the preferred outputs.
One such possible weighting scheme could be to assign weights to reflect
the nonlinearity of an equation in a certain term. For example,
consider the following set of indexed equations:
f. = 0 = x +2x. 3(lnx  )+3x 14 41
11 1 11 1 x. il+1
11
For any value of ii, it would be considerably easier to solve eqn. 41
for x. or x.1 than for x This can be reflected in the IDM by
111 11+1 11
assigning the incidences corresponding to i3 a much higher weight than
the incidences corresponding to either ii or i1+i. The resulting IDM,
for a blocking factor of 6, would then appear something like the one in
Fig. 43 a) (for minimization of the sum of the output weights). The
outputs could then be as in Fig. 43 b) or c). The assignment of
weights is a subjective matter, but it is one area where engineering
intuition can be employed to guide the selection of the solution
procedure.
4.3 FunctionVariable Output Set Assignment
The assignment of variable type outputs is made by treating the
FVIM as an incidence matrix and using the same algorithm as was used
for the index outputs. The aim is to employ the maximum output product
criterion developed by Edie (1970). The maximum product criterion is
used in an attempt to enhance the convergence properties of the system
of equations. It states that the outputs should be assigned in such a
way as to maximize the product of the sensitivities of the functions to
their chosen outputs. Then sensitivities become weights for the inci
dence matrix elements. The maximum product of the output weights can
be achieved by minimizing the negative of the logarithms of the output
weights. The real problem is to assign meaningful weights to the
elements in the FVIM. What is actually desired is to assign outputs .to
the expanded problem, each element having a weight proportional to its
corresponding element in the Jacobian. If this were done the restric
tions on FVIM output set assignments would most likely be violated.
The procedure adopted for calculating the FVIM weights is the
following. The Jacobian is calculated for the expanded set of
b) __
9
9
C
It 9G
c)
I 9
FIGURE 43
91
i
S9 1
0
90
19 0
INDEX DISPLAY MATRICES WITH WEIGHTS
1 9 1
9 9
1 9 1
1 9
1 9 1
191
equations, with index limits defined by the blocking factors. For the
calculation of the Jacobian, estimated variable values are used. The
index ordering for the Jacobian has function type nested outermost,
as indicated in Fig. 44. For each functionvariable partition in the
Jacobian [Fig. 44 a)], an output product is computed by multiplying the
elements corresponding to the index outputs for that function in that
variable. If a variable type is not a legal output for a function type,
the corresponding product is set to a very large positive number. These
products are then the weights used in the FVIM output set assignment,
as is shown in Fig. 44 b).
While this procedure does not necessarily generate the true
maximum product output set assignment, it does generate the most nearly
maximum product output set assignment possible under the restrictions
for FVIM outputs.
4.4 Minimum Weighted Tearing
Tear variables for a solution procedure with assigned output
variables are chosen using the minimum weighted tearing algorithm of
Pho and Lapidus (1973). If all variable weights are equal, the tearing
scheme chosen will be one with the minimum number of torn variable
types. The tear variables, however, are chosen by analysis of the FVIM.
The columns in the FVIM do not necessarily all represent the same
number of variables. The number of variables of a variable type is
determined by the function index definition (not the variable index
definitions, which may include decisions). The expected index ranges
can be used to calculate the expected number of variables for each
variable type. This number can then be used as the weight for the
x V z
af af af
a xl a vT a zT
a xT aT a zT
a XT avT aZT
ifi ,. ,in
x y z
,4 
y( I
5 I
___ l

FIGURE 44
FUNCTIONVARIABLE OUTPUT SET ASSIGNMENT
b)
variable type. The weighted tearing algorithm can then be used to
determine the minimum number of actual torn variables directly from the
FVIM. The choice of tear variables, like the choice of output
variables, is restricted to variable type.
Another possible means of assigning weights to variable types is
to have the user assign them. If the user has any familiarity with the
problem, he may know that some variables are more difficult to guess
accurately than others. He could then assign high weights to those
difficult to guess. The GENIE system presently allows either of these
two methods of assigning weights.
An alternative method of assigning tear variables would be to
perform a sensitivity analysis on the functions and variables. Then,
the tears could be chosen so that the functions would be relatively
insensitive to bad guesses for the tear variables. This appears to
involve an expensive combinatorial analysis and is not at present
implemented.
CHAPTER 5
COMPUTER IMPLEMENTATION
Since a detailed description of the means of computer implemen
tation of GENIE would require a detailed description of GENDER, the
description of the computer implementation given here will be restricted
to a discussion of the data structures developed for indexed equations,
a discussion of the memory system developed for the GENDER system, and
a general flow diagram and description of key subroutines.
5.1 Data Structures
In designing data structures for use in the GENDER system the goal
has been to use as little space as possible while producing a structure
which is compatible with the algorithms which will use it. In order to
do this there has been an attempt to minimize the number of list pro
cessing type data structures and to use the addressing capabilities of
FORTRAN to minimize searches. In fact, there is only one list of a type
which must be searched in the data structures which contain all of the
index data for the functions and variables.
5.1.1 Index Data Structures
The first major data structure necessary to the analysis of
indexed equations is the data structure which contains the information
about the function and variable indices. This data structure is the
"INDEX" COMMON area. Two of the vectors in this data structure contain
definitions of the function indices (FIDV) and of the variable indices
1
(VIDV), where both are vectors of integer variables.
The first of these, FIDV for Function Index Definition Vector,
is set up as illustrated in Fig. 51. The length indicates how many
words of the vector are occupied by data. The pointer to VS points to
a cell within the vector "VS" which contains the name of the index. The
next three entries in an FIDV cell, L., Ui, A. are the lower and upper
limits and increment for the index. The range limits are suggested
values only since the solution procedures developed are for arbitrary
index ranges. The suggested range indicates typical values for the
indices in problems which will be solved by the generated solution
procedure. The entries L.B and U define the blocking factor and are
I 1
the range limits used to produce the Index Display Matrices. The "FPO
pointer" is a circular linked list pointer which links together all of
the function indices which fully precedence order. A scalar variable
"IFPO" serves as a pointer into the circular list so that it can be
searched. The last entry in an FIDV cell is a pointer to "IDL," the
index decision list vector. If this variable is equal to zero it does
not point to anything. The vector IDL will be explained later.
The second index data vector is VIDV, the Variable Index Defini
tion Vector, which is illustrated in Fig. 52. As in most vectors used
in the GENDER system, the first word of VIDV is a length used to
indicate how much of VIDV is occupied by data. The size of a VIDV cell
is not fixed; however, the first six words always serve the same
function. Each of the VIDV cells represents a different variable index.
The first word of a VIDV cell is a pointer to a cell within FIDV and
thus indicates upon which function index the variable index in question
depends. The second word of the cell serves as a pointer to "IDMPV,"
F I DV
Length
Pointer to VS
Lij
Lil
Ui
Ai1
FPO Pointer
Pointer to IDL
Pointer to VS
Li2
Ui2
Ai2
L 1
FPO Pointer
Pointer to IDL
  1
etc.
FIGURE 51
FUNCTION INDEX DEFINITION VECTOR
FIDV Cell for i1
FIDV Cell for i2
/
VIDV
Range
L Flag
L Offset
U Flag
U Offset
Increment
List
No. of Offsets
First Offset
Weight
Second Offset
Weight
etc.
FIGURE 52
VARIABLE INDEX DEFINITION VECTOR
Length
Pointer to FIDV
Pointer to IDMPV
Output Flag
Output Offset
Range or List Flag
Direction Flag
Z hng
the Index Display Matrix Pointer Vector; for, as mentioned earlier,
analysis on variable indices will largely be performed on their Index
Display Matrices. The third word serves as a flag which indicates
whether or not the variable index has been output set assigned. If it
has, the next word indicates what the output is if it is a simple offset
from the function index. If the output is not a simple offset the Index
Display Matrix must be examined to determine the index output. The
fifth word of the VIDV cell is a flag indicating whether or not the
variable index is defined by a range or a list. The sixth word
indicates whether the index is to be incremented in ascending or
descending order. The rest of the VIDV cell depends on whether the
index is defined by a range or a list.
If the index is defined by a range, the length of the VIDV cell is
fixed, as there are five more entries for a range defined index. The
first entry indicates whether the lower limit of the range is offset
from the lower limit of the function index or from the function index
itself. The next entry indicates what that offset is. The next two
entries serve the same purpose for the variable index range upper limit.
The last entry is the range increment. If the index is defined by a
list, the length of the VIDV cell is variable. This is because the
number of list elements is variable. Within the list definition section
of the VIDV cell, the first entry indicates how many list elements
there are. Subsequently each list element occupies two entries, the
first of these being the offset from the function index value and the
second being a weight which is used by some of the analysis routines.
To facilitate analyses performed on function and variable indices
and to prevent searches, there are two vectors of pointers, one
associated with FIDV and the other with VIDV. The form of both vectors
(VIPV and FIPV) is the same, and they will be described together. As
usual the first word is the length of the data in the vector. The
second word is a pointer to the first FIDV or VIDV cell, the third word
is a pointer to the second cell, etc. The entire data structure as
described so far fits together in the fashion illustrated in Fig. 53.
Satellite data structures to the index representation data
structures are the vectors "IDMPV," the Index Display Matrix Pointer
Vector, and "IDL," the Index Decision List. The vector IDMPV is
accessed through VIDV and merely provides the name of the Index Display
Matrix for each of the variable indices in VIDV. The name is not stored
in VIDV because that would complicate searching through the Index
Display Matrix should that be necessary. The index decision list is
the vector which contains the index decision declarations. Index
decisions are offsets either from the lower or upper limit of the
function index. Each function index has associated with it an IDL cell
which is composed of two subcells, as shown in Fig. 54. If a function
index has no index decisions, its IDL cell will have no entries. Each
IDL cell has three length specifications, the first stating the length
of the entire cell. The other specifies the length of each of the sub
cells. The first subcell, the L offset subcell, specifies the offsets
of decisions from the function index lower limit, and the second
specifies the offsets of decisions from the function index upper limit.
The index decisions can be declared in two ways; either the user can
declare them or algorithms can assign them. Declared decision variables
can never be changed, while derived decisions can be changed. In order
VIPV
etc.
FIGURE 53
INDEX DATA STRUCTURE
VIDV
ID L
Length
I Length
Length
L Offset 1
L Offset 2
etc.
Length
U Offset 1
U Offset 2
Setc.
etc.
L Offset Subcell
U Offset Subcell
FIGURE 54
INDEX DECISION LIST
IDL
Cell
to distinguish between these two types of decisions, derived decision
offsets are incremented by 10,000. Since decision offsets will never
be that large, large offsets are easily recognized as derived and are
adjusted to their true value when necessary.
5.1.2 FunctionVariable Incidence Matrix Representation
It is necessary that the FunctionVariable Incidence Matrix be
stored in a data structure which allows for easy access by the algo
rithms. Since the FVIM is a type of incidence matrix it is reasonable
to store it in the data structure used to store incidence matrices. The
incidence matrices are stored in a data structure called the SIM data
structure (Cunningham, 1973). Appendix A contains a description of the
SIM data structure. The principal difference between the normal
incidence matrix and the FVIM is the fact that an element in the FVIM is
actually a list of variable indices. To allow for that it is necessary
to modify the element entry in the normal SIM structure to accommodate an
extra word. This is easily done since the length of an entry is
variable in the SIM data structure. This extra word serves as a pointer
to the vector "VILIST." Thus, each FVIM entry has associated with it a
pointer to VILIST. The structure of VILIST is illustrated in Fig. 55.
Each FVIM element pointer points to a different VILIST cell. The VILIST
cells are merely lists of pointers to variable index definition cells.
This results in each element in an FVIM having associated with it a list
of variable indices, which is what is required.
5.2 A Versatile Memory System
A part of the design philosophy of the GENDER system has been to
construct the different data structures necessary in separate labeled
COMM0N areas. For example, a special data structure was developed for
VILIST
Length
Length
Pointer to VIDV
Pointer to VIDV
Pointer to VIDV
Length
Pointer to VIDV
Pointer to VIDV
etc.
VILIST Cell
FIGURE 55
VARIABLE INDEX LIST VECTOR
the storage of sparse incidence matrices. This data structure is built
into a set of labeled COMMON areas called "SIM COMMON." At any time
only one sparse incidence matrix (SIM) can occupt this common area.
Quite often during the analysis performed by the GENIE system it
becomes necessary to temporarily stop working with one SIM and start
working with another, without losing the data in the first. Since all
routines operate on SIM's in the SIM COMMON area, it is necessary to
move the first SIM out of common and at some later time restore it to
the common area. This type of movement is not restricted to SIM's and
is desirable for many different data structures. It became necessary to
design a general purpose memory system, easily assessible to any of a
variety of routines which need to store data in a relocatable fashion,
so that the data can be sestored to the appropriate data structure when
required.
The memory system which accomplishes this is itself a labeled
COMMON area called, appropriately enough, "MEMORY." The memory system
consists of two major storage areas, or vectors, the larger being
"MEMORY" and the smaller being "MEMDIR." MEMORY is the area of core
into which data to be saved are moved. It is a long vector (in the
first version 5,000 words) and could conceivably be much longer. MEMDIR
is a directory to the vector MEMORY which contains all of the informa
tion necessary to relocate a data item stored in MEMORY. A third major
storage area associated with the memory system is a random access mass
memory data set onto which data from MEMORY are written should MEMORY
ever not have enough space available to store a data item. It should be
noted that many different data items will in general be stored in MEMORY
at the same time, and that if enough are stored MEMORY will become full.
It is perhaps easiest to consider ME: !,' to be a buffer for the mass
memory device. The memory system then can be considered to be a mass
memory system using relatively slow speed mass memory with a higher
speed core resident beffer. Access to stored data is made through the
buffer. The memory system handles all access to the mass memroy.
Although to the user of the memory system it appears that all data
are kept in MEMORY, an understanding of the interaction between MEMORY
and mass memory helps in understanding the system. The random access
data set has associated with it a logical record length (1,000 words in
the case of the memory system). The data set is divided up into a
number of these logical records. In this particular implementation,
written for IBM 370/165, the random access data set has data written to
or read from a particular record on a disk. This makes it possible to
access a particular record. Since the records are numbered sequen
tially, the words in the first record can be thought of as being words
1 to 1,000 and the words in the second record as 1,001 to 2,000, etc.
Figure 56 illustrates this addressing scheme. If the address of a word
of data is known, then the record of the data set it is stored in can
easily be determined.
In the memory system the first record is reserved for a special
function which will be described later and is never used for data
storage. Data storage begins with the second record. Data which are
moved to MEMORY have associated with them absolute addresses, which are
the word addresses that the data would have were they stored onto the
random access data set. This can most easily be visualized by imagining
that the vector MEMORY is superimposed on the random access file in Fig.
56. The first word of MEMORY would have absolute address 1,001, the
RANDOM ACCESS
DATA SET
Word 1
Word 2
+ Word 1000
 Word 1001
+ Word 1002
Word 2000
FIGURE 56
STRUCTURE OF RANDOM ASSECC DATA SET FOR MEMORY SYSTEM
Record 1
Record 2
I
e.
etc.
second word 1,002, etc. The vector ME i RY has two states, one indi
cating that the data in MEMORY are not actually stored in the random
access locations corresponding to their absolute addresses, the other
indicating that the data in MEMORY are also in random access memory.
Until data in MEMORY have been written to the random access file, the
first word of MEMORY will correspond to word 1,001 of the random access
file. Suppose, however, MEMORY is 5,000 words long and that 4,500 words
of MEMORY are being used to store data and a request to store a data
item of length 1,000 words is made. Since there are only 500 words
available for data, some of the data in MEMORY will be written to the
random access file. Actually all of MEMORY is written to direct access.
Since data are accessed by 1,000 word records, the first 4,000 words of
MEMORY are considered as written and the next 500 words are moved to the
first 500 locations in MEMORY. The first word of MEMORY then corresponds
to absolute location 5,001. This would correspond to superimposing
MEMORY over the random access data set beginning with the sixth record.
In this manner considerably more words than can be stored in MEMORY can
be stored and addressed.
When a data item is stored into MEMORY it has a directory entry
created for it so that it can be referenced through the directory. A
directory entry consists of six words, allocated to the following
purposes:
WORD CONTENTS
1 Number of Characters in Name
2 First 3 Characters in Name
3 Second 3 Characters in Name
