Title Page
 Table of Contents
 Indexed equations
 Solution procedures for index equation...
 Algorithms for deriving solution...
 Computer implementation
 Convergence properties
 Example problems
 Conclusions and recommendation...
 The sim data structure
 Subroutine descriptions
 Biographical sketch

Title: Automatic generation of solution procedures for indexed equation sets using GENIE
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00084181/00001
 Material Information
Title: Automatic generation of solution procedures for indexed equation sets using GENIE
Physical Description: vii, 160 leaves. : illus. ; 28 cm.
Language: English
Creator: Allen, Gary Louis, 1948-
Publication Date: 1974
Subject: Electronic data processing -- Equations   ( lcsh )
Algorithms   ( lcsh )
GENIE (Computer program)   ( lcsh )
Chemical Engineering thesis Ph. D
Dissertations, Academic -- Chemical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis -- University of Florida.
Bibliography: Bibliography: leaves 158-159.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00084181
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 000580528
oclc - 14046231
notis - ADA8633

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
        Page vii
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
    Indexed equations
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
    Solution procedures for index equation sets
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
    Algorithms for deriving solution procedures for sets of indexed equations
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
    Computer implementation
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
    Convergence properties
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
    Example problems
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
    Conclusions and recommendations
        Page 131
        Page 132
        Page 133
        Page 134
    The sim data structure
        Page 135
        Page 136
        Page 137
        Page 138
    Subroutine descriptions
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
    Biographical sketch
        Page 160
        Page 161
Full Text







To Jan



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 GK-18633 and GK-41606. 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.





1 INTRODUCTION..................................

2 INDEXED EQUATIONS.............................

2.1 Conventions Used in Describing Indexed
Equations ...............................

2.2 Function-Variable Incidence Matrices....

2.3 Index Display Matrices..................

2.4 Index Imbedding.........................

2.5 Output Set Assignments .................

2.6 Acceleration of Variables..............


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...............................

OF INDEXED EQUATIONS ...............................

4.1 Decoupling..................................

4.2 Index Output Set Assignments ..................



4.3 Function-Variable 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
Gauss-Seidel and Newton-Raphson................

7 EXAMPLE PROBLEMS .............................. ......

7.1 Distillation Model .............................

7.2 Examples ........................... ......... .

7.3 Discussion................................... .





BIBL IOG APHY..................................

BIOGRAPHICAL SKETCH ...........................














1 !







. .. .
. 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



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 non-convergent, 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 "do-loops," so the index ranges become the

do-loop 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.



The design or simulation of a chemical process requires finding

the solution to a large set of equations (many of them non-linear).

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, non-linear, 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
"'solu-tion procedure."

Two commonly employed types or solution procedures are tearing

procedures (such as Gauss-Seidel) and "iNewton" type procedures (such as

Newton-Raphson). Because of the cost (usually in computer time) of

solving large sets of equations simultaneously, it is desirable to avoid

Newton-Raphson 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

(non-recycle) 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 non-indexed

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


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



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.



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 2-1
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.

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

multi-unit 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 2-1. 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


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


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 do-loop range and consequently makes the eventual writing

of the solution procedure in the form of do-loops 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 2-2
1 i-i" +1
represented by a tri-diagonal incidence matrix, would have its index

sets described as:

i: L = 1 j: L = i-i

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 2-3
S1i-2 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 2-3 would have

length 3 and would be:

j: ii-2


There is some overlap between index lists and ranges. For example,

either a list or a range can be used to describe the tri-diagonal

relationship in equation 2-2. 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 Function-Variable 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 "Function-Variable Incidence

Matrix" (FVIM). The FVIM is not a true incidence matrix, and because

of that has some unusual properties. A Function-Variable 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 2-4
1112 i1-1,i2 1ii2 112 t1-2+1

g = 0 = x. +y -y
gil1i 112 i li 2 +1,1i2

The Function-Variable Incidence Matrix would be:

x y

g(il,i2) 3?j4 1 324

__-- I

where I iL-1 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


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 2-5
1112 11-1,i2 1li2 ,1i2 11,i2+1

One might be tempted to assign four variable index names as follows:

iJl i-1 j2:' il J3: i2 j4: i2

i- ii+1 12+1

This, however, does not accurately reflect the structure of the

equations 2-5 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
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 Function-Variable 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 Function-Variable 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
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): iu-1

U = 20 ii

A = 1 il+i

the Index Display Matrix (with zero elements blank and non-zero elements

denoted by 'x') would be that of Fig. 2-1. 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

9i XXX

10 X XX
V \e V

14 X XX
15 XXX
16 X XX"
17 X X X
18 X, XX
2O1 ^ ^



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 tri-diagonal 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 Function-Variable 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

"do-loops." 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. 2-2. If we

consider the three indices (function type, ii, and i2) as do-loop

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 do-loop 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

Function-Variable Incidence Matrix

x y

a) f(il,i2) jlj2 J4j3

g(il,i2) jlj3 jlJ3

Index Definitions

i,: L = 1


A = 1

Jl: L = 1

U = ii


J4: L = 1

U = U.

i: L = l

j2: i2-1





i3: i2


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 2-2 (cont.)

(i.e., the order of nesting) in terms of do-loops makes the transition

to computer implemented solution procedures in FORTRAN using do-loops

straightforward. The advantages of using the do-loop 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 do-loops or code similar to do-loops.

In order to fully realize the consequences of the nesting of the

indices and the effect of the Function-Variable Incidence Matrix and

Index Display Matrices on the actual incidence matrix, examine Fig. 2-3.

This is the expanded incidence matrix for the problem in Fig. 2-2.

First look at Fig. 2-3 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, non-zero 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. 2-3 a), is itself

partitioned along boundaries between the indices nested next innermost

(corresponding to i2). Figure 2-2 d) exhibits the s;me structure as

Fig. 2-3 b) (in a blockwise nature). This is what might be expected

since Fig. 2-2 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. 2-3 c)

exhibits: the same structure as Fig. 2-2 b). Figures 2-4 a)-e) are

representations of 2-3 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 -,.




O000 011e0
0 N 8 0 N t 0oN N
x t >s >c x: x x x x f e x
X X iX X I

v X
x x
XX Xxx
S -I d


21 IX X
f3! -XXX-



f X :X

N ck-




x %If r %



,_,__ 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
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


if, ix i2

0- N,. -> o0 o3--I 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
f221 XX XX XX x x

T2 x x xx xx xx x x
fa xx xx xx x 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



o0- -- c~ MC r0 -- WM --w "z 0-- J rl

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

'i l2,i'f

g l
g 12
g 32

i ,if i

FIGURE 2-4 (cont.)

o_ 0 c\ _r- ro

x x xxI x xxx

x x X X X X xX

x x x
xx xx x xx
X X X X X X X X X x



i2viI lif

FIGURE 2-4 (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


I Xxxx X X X X

.. x. x x XX x x x x x




behavior, called "Index Imbedding," allows the structure of the full

incidence matrix to be predicted, and in fact visualized from knowledge

of the Function-Variable 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 do-loops. As

an example consider Fig. 2-5, 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 do-loops.

If, however, we assign outputs to the Function-Variable 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.

s A / '^ 1 '*" i/;
A 4> (I) A_

IR 2-5
> 0 a X A
/.* ^ ..' >. A t
~ v \ *' ,^ *

^ ^~\ \'\" \ "f
E 1 .,-*.} A A A ; *Vn-




21 2



.i ;-'

tions on output set assignments consider the system of equations from

Fig. 2-5, defined structurally by


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



with IDM shown

index outputs

the following












in Fig. 2-6. Assigning x to f

indicated by the circles on the

incidence matrix:
xll X12 X21 X22 Y11 Y12 Y13 Y91

and y


( 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. 2-6, we have


ii: L




= 2


i2: L = 1

I 2


L 2 3
x o


I1 2

I x x



I 2




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 do-loop type iterative

solution procedure such as:


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)



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 do-loops, 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 non-indexed

equations. Finally, by assigning the outputs to the small sub-matrices

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.


For n=3 the IDM is




1 2 3 4 5 6 7


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


1 2 3 4 5 6



ii 5



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


Define the function index range as:

L = M

U = A.(N)-(A.-M)

A = A.

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.+(l-a)i+ki aE{O,l}

U. = bU.+(1-b)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

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.

= N

The number of variables

U. -L.
n= max min +1
v A.

U3max =bU,+(1-b)imax+k.u
i 3
but imax = U. so,

U = bU +(l-b) U+k .
gmax 2 2 g
= U .+k, .


L = L .+k .'
rmin J 7

U .+k .-(L .+k )
so n +1

U.-L. k.u-k .
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.



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


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 "do-loop" and nest the do-loops. Figure 3-1

illustrates a typical nesting of do-loops. 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 do-loops

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 do-loops 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 /




---Q. I


AT= 1210, 1

12= I,5,

I3= I, 7,v




loop I

loop 2

loop 3

end loop 5

end loop 2

end loop I


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


f(il,i2) = 0 = f[x(j1,j3),y(jl,ji),z(jl)]
g(il,i2) = 0 = g[x(j2,j4),y(j2,j4),z(j2)]

The decomposed representation of these equations is shown in Fig. 3-Z,

which also defines the indices shown in the equations. To solve these

equations, which may be assumed to be non-linear, 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


f ( i,i, )


1, L=


jr: j1 -l


A= I

i2: L=I


j- : L= Li
U= Ui2

SD s


," w C-, !

; v/ 1!



J2: ii

J4: i2


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. 3-3. 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


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 Function-Variable 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

I I, I


- 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 xx xx x
X VX ," V ^7
A ^f\,.f) /A /
,, V XX v /^? v

I. L~ Y





f l



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.
LM. 0 = x..-l 3-2
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 Newton-Raphson 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. 3-4 a). There are three function

types and four variable types. Choosing x as the decision variable

type results in the FVIM in Fig. 3-4 b) and choosing z results in

Fig. 3-4 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 3-1: 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,



f(il ,i2) J, J4 J2 J4 j2 i3

g (1,12) Ja J5 j I 4 J

h(i ) j j5 ij j2
= .



h(i ,)




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



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)+(l-a)L.+kz a{O0,1}

U = U. U = b(i)+(1-b)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 non-null k must be greater than or equal to

k For one row the IDM is output set assignable since it is non-null.

Now assume that the IDM is output set assignable for k rows. The number

of rows can be expressed as:

r i
k A.

and the number of columns as

b(imax)+(l-b)U.+k -a(imin)-(l-a)L -k~
uc A --+ 1 3--3

but imadx = U,

and imin = L.

so eqn. 3-3 becomes

Ui -L +k -ki
k u + 1

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 3-4
= k +

so there is a new column which can be assigned as the output for the new

row. Theorem 3-1 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 3-2: 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 3-2 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,
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


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. 3-5, 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




Ex xxB
L-\ -{4
('A ^-l* 'A*ut* ^^'w<^^^
t: S.? V
6 (1S/ A 5

bs) -,



'A, A-

)r '< i
A' A,
'\ J

S',: ".



7 ,

'A. 'A I 'A 'A 'A

.', 5.- ', s

s s
'A A

\f 'V v' '1

' 1
' X( j
/* S 59 -
- s.

x'x x I. x

6 f '
I i '^ ."t Ii t r

S'V 'P v V .
* -J \.
a '


Cm w. wVZX-I-llW YU- rm

--r-_-^ --~- --Y



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 3-3: 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, .
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. 3-5
i2=1 i= =1

Mi and M2 are the number of rows (and columns) in the IDM's for i1 and

i2 respectively. Equation 3-5 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 3-6

since all functions must be evaluated. Equation 3-6 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
N = MM, I k. 3-7
T i 2

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


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. 3-6 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


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 e-aev *mm "I |e.e
f. '' 4

I I. I
H fi i>
*j :* 'f

Fj I

h F .. Li V
i, t-


tt~~' '"-~



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. 3-6. 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. 3-7. Both incidence matrices can be solved with four


Theorem 3-4: 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
S v B



'.^a J

4 o

XX / "

x x x x

xxx x x xY
X X X X v !



x x":x x:
j ^_^ x x'



b) tJ

L [ /u

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-


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. 3-8. 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. 3-9 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 3-5: 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

2 r x x

Tx x xB
%x X x~

x x
2x xX







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

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


ii. For il=L. all incidences of x in functions other than f which have
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. 3-10.) In this

figure, a) and b) are Index Display Matrices representing a tri-diagonal

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


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. 3-11, each representing one of the one-tear, blocking

factor of 3, solution procedures from Fig. 3-10. 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




l x









x x

X 6)
A A l i

xl x

x x






procedure for a blocking factor remains effective as the IDM expands,

the solution procedure is said to "extend."

Theorem 3-6: Call a block with blocking factor equal to k a k-group.

For the solution of k-groups 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 k-group. An analogous condition exists for

solution of k-groups 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 k-group tear variables by k. When the next k-group is

added, the decisions offset from U. are no longer decision variables.
In order to solve the next to last k-group, 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

k-groups it holds for n+l k-groups. It must hold for one k-group,

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 k-group. When the next

k-group is added, this variable's value must be available to solve the

next to last k-group. It is no longer a decision and is not a tear, but

is an output of a function in the last k-group 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, k-1 tears are made in the first block and k

tears must be made in the second block. Next consider L. offset from L.
and U. offset from i, again by the same factor. A full precedence
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.



The algorithms implemented in the GENDER system for deriving

solution procedures are capable of treating non-indexed 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-


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 function-variable 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


III. A. If the set of unassigned function and variable types does

not fully precedence order, go to IV.A., otherwise,


B. A full precedence ordering has been discovered. Stop.

C. Decoupling will not produce a full precedence ordering.


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 4-1 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.


f(iI ,i2)

g(i, ,i2)










L= Li

L =Li,

A= I


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. 4-2. 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



x x x xI I ,
,x X jx
x x x xx
r V-

xx xxx x
x x x
I x x x

I xxx x@q







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


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 non-linearity of an equation in a certain term. For example,

consider the following set of indexed equations:

f. = 0 = x +2x. 3(lnx -- )+3x 1-4 4-1
11 1- 11 1 x. il+1
For any value of ii, it would be considerably easier to solve eqn. 4-1

for x. or x.1 than for x This can be reflected in the IDM by
11-1 11+1 11

assigning the incidences corresponding to i3 a much higher weight than

the incidences corresponding to either i-i or i1+i. The resulting IDM,

for a blocking factor of 6, would then appear something like the one in

Fig. 4-3 a) (for minimization of the sum of the output weights). The

outputs could then be as in Fig. 4-3 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


4.3 Function-Variable 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) __


It 9G

I 9


S9 1

19 0


1 9 1
9 9

1 9 1
1 9
1 9 1

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. 4-4. For each function-variable partition in the

Jacobian [Fig. 4-4 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. 4-4 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





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




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


(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. 5-1. 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. 5-2. 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,"


Pointer to VS

FPO Pointer
Pointer to IDL
Pointer to VS
L 1

FPO Pointer
Pointer to IDL

- -- 1



FIDV Cell for i1

FIDV Cell for i2



L Flag
L Offset
U Flag
U Offset

No. of Offsets
First Offset
Second Offset




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. 5-3.

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 sub-cells, as shown in Fig. 5-4. 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 sub-cell, the L offset sub-cell, 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







I Length
L Offset 1
L Offset 2
U Offset 1
U Offset 2


L Offset Sub-cell

U Offset Sub-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 Function-Variable Incidence Matrix Representation

It is necessary that the Function-Variable 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. 5-5.

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


Pointer to VIDV
Pointer to VIDV

Pointer to VIDV
Pointer to VIDV

Pointer to VIDV





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


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 5-6 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.

5-6. The first word of MEMORY would have absolute address 1,001, the


Word 1
Word 2

+ Word 1000
- Word 1001
+ Word 1002

Word 2000


Record 1

Record 2


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



1 Number of Characters in Name

2 First 3 Characters in Name

3 Second 3 Characters in Name

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs