IEEE TRANSACTIONS ON CIRCUIT THEORY, VOL. CT- 18, NO.

1, JANUARY 1971 101

The Sparse Tableau Approach to

Network Analysis and Design
GARY D. HACHTEL, MEMBER, ROBERT
IEEE, K. BRAYTON,
AND FRED G. GUSTAVSON

Abstract-The tableau approach to automated network design LkP,)

optimization via implicit, variable order, variable time-step integra-
!FF--l+
tion, and adjoint sensitivity computation is described. In this ap- E(t) 7 y 43’ *(%~“z)
preach. the only matrix operation required is that of repeatedly solv-
ing linear algebraic equations of fixed sparsity structure. Required Fig. 1. Nonlinear network.
partial derivatives and numerical integration is done at the branch
level leading to a simple input language, complete generality and
maximum sparsity of the characteristic coefficient matrix. The bulk
of computation and program complexity is thus located in the sparse i5z+(i3,pI)
matrix routines; described herein are the routines OPTORD and l-2-3 wqqsvq
GNSO. These routines account for variability type of the matrix ele- n “,;+ jJk pyy
ments in producing a machine code for solution of,Ax=b in nested
iterations for which a weighted sum of total operations count and
T IT
round-off error incurred in the optimization is minimized. OGZGG PARA;TOR

Fig. 2. Augmented network.

I. INTRODUCTION
OMPUTER-AIDED design of electrical networks
has recently been undergoing a rapid evolution due
to the development of highly efficient methods for
solving systems of linear equations where the matrix of
Fig. 3. Node-branch labeling for the directed graph
coefficients is sparse [l]-[4]. These new techniques, when
fully absorbed in a computer-aided design program, dras-
tically affect even the initial formulation of the problem. of a tableau which includes all network information in a
In this paper we describe a research investigation aimed nonreduced form.
at designing and implementing a computer network design For illustration consider the automatic design of the
program (NDP) for automated network optimization network of Fig. 1. The nonlinear elements were chosen to
which fully incorporates the sparse matrix methods. The demonstrate that there are no restrictions on element type.
formulation is unorthodox yet simple and is applicable to In general, we consider design objectives of the class
any system described mathematically by a set of algebraic T

and differential equations. A general nonlinear electrical minimize Q = $+J> P, t) dt (1)

network is such a system, and it follows that the basic pro- PWI~P~PM s0
cedure for solving for the currents and voltages in a network where @ represents a generalized scalar objective functional
can be viewed in terms of a simple mathematical procedure : and 4 an objective function of the unknown vector w and the
Gaussian elimination. This has both a simplifying and vector of designable parameters p. In Fig. 1, $=4(i3, pi).
unifying effect and makes possible a simple yet truly gen- It is convenient in automated network design to add these
eral-purpose computer program. By combining the concept quantities directly to the tableau of unknowns. This is done
of variability type with sparse matrix techniques, this pro- (Fig. 2) by augmenting the original network with a unit
gram can achieve practically optimum efficiency by choos- capacitance driven by a dependent current source of value
ing the order of elimination so as to minimize the total op- 4 (the integrand in (1)) and another unit capacitance con-
erations count required for simulation .and/or automated nected only to the reference node. The reason for this last
optimization. augmentation has to do with adjoint sensitivity computa-
In this section we give an overview of our method for tion [5] and will be made clear in Section II. The first aug-
orientation purposes. Basic to our approach is the concept mentation (cf., Fig. 2) is called herein an objector and the
second a parametor. The directed graph of the augmented
network is shown in Fig. 3 with integer node and branch
Maduscript received May 27, 1970; revised August 10, 1970.
The authors are with the IBM T. J. Watson Research Center, York- labels. Fig. 4 shows the tableau of algebraic and differential
town Heights, N. Y. 10598. equations which describe the network, where I/, i, v stand
102 IEEE TRANSACTIONS ON CIRCUIT THEORY, JANUARY 1911

= fl V. i, Y.4. 0, p, tl = 0

Tableau Operator

Tableau error vector

Fig. 4. Tableau operator, unknown vector, and tableau vector.

for the vectors’ of node voltages, branch currents, and

branch voltages, respectively. Energy storage is denoted by
the generic term 4 which stands, in network applications,
either for magnetic flux or electronic charge. However,
since the objector and parametor charges, are of special Fig. 5. Tableau matrix.
interest, the labels Q and p are retained for these elements.
The tableau operator F(. , . , . , . , . , . , t) (Fig. 4) is an alge-
0 A’ 0 0
braic-differential operator which operates to the right on the
unknown vector w c co1 (V, i, u, ~7,c$,p) to form the tableau
A 0 -1 0
error vector f(w, t). The network equations may thus be
stated in the form
0 .z Y a

0 E ‘%+
a
The tableau matrix F,(w, t) is the linearization of the tableau Fig. 6. Form of general tableau matrix
operator F around the operating point w, i.e.,
ql(=“K w+I+J=%+,tJ

+N+l+J
(3) T K=I,N J=l.“p
ihl.~=+(q.P.t)
EQUATOR O= PARAGTOR

Fig. 7. Network analog of general ODE optimization problem.

The tableau matrix for the network of Fig. 1 is illustrated in
Fig. 5. Figs. 4 and 5 show that (for this case) the first ten
substitution:
rows represent the Kirchoffs current (A’i=O) and voltage
(Al/- u=O) laws. The matrix A is the branch versus node
incidence matrix [6], which has the advantage of being
directly obtainable from the network description input data
cards (which also establish the labeling of Fig. 3.) Rows
11-16 of f (Fig. 4) state the branch constitutive relations
(one per branch) of the network, while rows 17-20 give the
dynamic relations describing the generalized energy storage The network analog of the general ODE optimization
elements (e.g., i, - cj3 = 0 (row 18)). problem is shown in Fig. 7, wherein each of the N com-
For general networks, the tableau matrix Fl has the form ponents gK, K = 1, N, of the system g(a, 4, p, t)=O is repre-
of Fig. 6, wherein for network applications submatrices sented by an equator, i.e., a unit capacitance driven by a
A, B, and E have all elements + 1 or 0. Note FI is very sparse, current source, i, = uK, and a floating current source of
averaging in general about three nonzero elements per row. value i,,=g,(u, 4, p, t). The objector and parametors (one
Note that the optimization problem (1) for any system of for each of the np designable parameters) are as previously
ordinary differential equations (ODE) of the form defined for the network problem.
g(Q,q, p, t) = 0 can be put in the form of Fig. 6 by the following The tableau approach to time-domain simulation is to
discretize, at the branch level, the derivative operator d/dt.
’ A lower case letter is a column vector unless otherwise specified, The nth discrete time point is denoted by t, and _.w,= w(t,).
” ’ ” denotes the transoose oeeration. Thus u’ is a row vector. Thus the last nq + 1 + np rows (nq charges, 1 objector, np
HACHTEL t?t al. : SPARSE TABLEAU APPROACH TO ANALYSIS AND DESIGN 103

parametors) have d/d replaced by the &h-order backward differentiation formulas (5) to algebraic-differential sys-
differentiation formula [7]. For example, for charges tems. This approach is similar to that developed by Gear
[8] for strictly differential systems except that -we store
backward differences instead of the Nordsieck vector. Ex-
cept for systems with only a small number of differential
where h = t, - t,- r. This transforms the tableau operator F, equations, the use of backward differences is more efficient
the tableau error vector f, and the tableau matrix F,, into computationally. This can be established by a simple op-
their discrete time forms Fd, fd, and FdP Here fd is obtained erations count. Section IV discusses the sparse matrix
from f by replacing 4,6, and Ij by sums of the form (5), and algorithms, and in Section V we show how the tableau fd,
Fdl is obtained from Fl (Figs. 5 and 6) by replacing d/dt & for the network optimization may be straightforwardly
by -cc,/h. set up from a SCEPTRE-likeinput language [9].
Thus the time-domain simulation problem (2) is trans- Sparse matrix packages OPTORDand t-2-3 GNSOare used
formed into that of solving the difference equations to invert F,, (actually the L/U factorization method is em-
ployed rather than direct inversion). Section IV gives the
w() = w(t,) = w(0) specific algorithms employed. The output of OPTORDis a
Fd(.3P, t”)w” = j&J”, P, t) = 0 (6) pivoting order, i.e., the order in which the rows and columns
of Fdl are eliminated in inverting &. The ordering selected
at each time point t, in the interval [0, T]. (Actually, fd gives an optimized compromise between 1) minimum
depends on w,- r,. . . , w,-~ because of (5), but this de- roundoff error in the solution Aw, and 2) minimum total
pendence is suppressed here for convenience of notation.) operations count in performing the double iteration (7).
A Newton iteration is used to solve the nonlinear algebraic 1-2-3 GNSO produces the required machine code SOLVEin
difference equations so that the overall computational task three partitions : C-SOLVE,to be executed at the beginning of
of time-domain simulation is the following double iteration. each optimization step, i.e., before entering the time loop of
t, = 0
n=O
w0 = initial condition
4 n=n+l
t, = t,.e 1 + h,
K=O
I wp = W”_ 1 (or wk”) is obtained by prediction, cf., Section II).
(7)
K=K+l
time iteration t Newton iteration FdI(w;K-l))Aw;K-l) = - fd(w;K-l’)
w;KIK)
= w;K-l) + Awf-I)
r
no Converged?
J yes
yes et,, I T?d

With only trivial modifications (7) applies to dc and ac (7) ; T-SOLVE,to be executed in the time loop of (7) ; and
cases also. In Section II we show how the routine for in- X-SOLVE, to be executed in the Newton loop of (7). All
verting Fd, may be used also to obtain the performance elimination of variables is done by the Gaussian elimination
gradient &D/iTp by the efficient adjoint approach [5]. routines. For example, the program does not require the
(With this approach the cost of computing the np-dimen- determination of state variables, trees, nodal admittance
sional vector X@p is roughly equal to that involved in (7)). matrices, capacitor loops, or any of the various special-
Thus the tableau approach achieves dc, ac, and transient ized network concepts which tend to complicate a basically
optimization with one matrix operation, namely, the in- simple procedure. The sparse matrix packages OPTORDand
version of Fd, in (7). Note (7) is essentially a two-step algo- l-2-3 GNSOallow us to dispense with these procedures. In
rithm. The first step is to evaluatef, and F,,. A virtue of the this sense the tableau matrix-Gaussian elimination ap-
tableau approach is that fd and Fdl may be read directly proach tends to have a simplifying effect on network anal-
from a simple input language using one input card per ysis not only computationally but to some extent theoretic-
branch. The second computationally critical step is the in- ally as well. In Section VI we demonstrate that nodal
version of the sparse matrix Fdl. This operation is handled by analysis [6] and state-variable analysis [lo] may be re-
package programs for sparse Gaussian elimination. The garded as special methods of performing Gaussian elimina-
tableau format allows these packages to be general-purpose tion with different pivoting strategies. From this point of
programs not special to the network application. In Section view there are many methods of analysis, and the most ap-
III we describe prediction and control of truncation error propriate one for purposes of computation is the one which
in a variable-order variable-time-step application of the gives the best combination of accuracy and efficiency.
 104 IEEE TRANSACTIONS ON CIRCUIT THEORY, JANUARY 1971

Comparisons of unconstrained OPTORD analysis with

OPTORD constrained to eliminate state variables or nodal
variables last are among the numerical results given in
Section VII. Also, results showing how storage require-
ments grow with problem size are given for a sequence of
examples.
II. THE TABLEAU AND ADJOINT SENSITIVITY COMPUTATION
The main computation in minimizing O(p) is to compute,
given p, Q’(p) and dQ/~?p. With these quantities, nonlinear
programming algorithms (e.g., Davidon-Fletcher-Powell
[ 111) can be utilized in minimizing Q@). In this section we
.,J 1
Fig. 8. Algebraic-differential format of tableau.
present a theorem which demonstrates the role of the tableau
matrix Fdl in the adjoint-method computation of &D/ap. It is
convenient for this purpose to group the algebraic network and 7 satisfies the adjoint algebraic equations
equations and variables so that
---- 1110) af,(o) o ~aa4
[a’, q’, @,p’] da ~’
a = col(l/,i,v) all ap
E 0 0 0
and a40 w(o) _ 1 w(o)
aa aq ap
da, 42 P, t)
Ea -4
e, 4, P, t) - 6
1> 0 I t I T. (8) 0 0 0 -1
Note that except for the d/dt entries, the coefficient

Thus
-P
matrices F,(t)
matrices in (11) and (13) are identical to the tableau operator
and F,(O), respectively. Thus the columns of
F,, operating to the left on the adjoint unknown vector
[a’, 4’, 6, $‘I give the adjoint sensitivity equations just as the
@ = (D(T) = ‘+(a, q, p, t) dt (9) rows of F, (cf., Fig. 4) operating to the right on co1 [a’, q’, Q, p’],
s0 give the original differential equations (8). Thus to solve (11)
and the format of the tableau matrix Fdl is as shown in and (13) as well as (8) involves Gaussian elimination applied
Fig. 8. only to Fdl (of Fig. 8), modified as shown in the d/dt entries.
In these terms we may now state, without proof, the fol- This theorem establishes the actual flow of the compu-

a(t,) t,
lowing theorem. tations leading to aQ/ap. First (7) is solved yielding
a(t,), q(t,), cD(t,) at the time points in the interval [0, T].
Theorem 1 Those components of and q(t,) which cause Fdl to vary
Assume that the system (8) starts at equilibrium, i.e.,
4(O)= 0, and that T is independent of p, i.e., aT/ap = 0. Then t,t,
are stored on disk. These are returned to core during
solution of (11). Since the times
do not generally coincide with the
involved in solving (11)
involved in (7), kth-
- = p(0) + 2
ap obtain a(t,) and q(tJ. The final value of Q’(t),
order interpolation formulas, similar to (5), are used to
i.e., g’(O), is
then passed into the right-hand side of (13) which is solved to
where f(t), 0 <t I T, satisfies the adjoint differential equa-
tions’ give p’. Finally, f(O) and 7 are added to give dQ/dp.
At this point both Q(J) and SD/ap are available, and an
age)
__ ~ado o am optimization step may be taken. It is important to emphasize
[ii’, q, 6, f] = [O,0, 0, 01
da aq ap that in the iterations involved in solving (7) and (11) only
E dldt 0 0 one matrix structure, that of Fdl, is subject to sparse Gaussian
elimination, and no other matrix operations, e.g., addition,
ada ada -d,dt w(t)
transposition, and multiplication, are necessary.
da a4 ap
0 0 0 - dJdt (11) III. DISCRETIZATION OF THE TIME DERIVATIVE OPERATOR
Two basic assumptions have guided the development of
with initial conditions
the tableau approach to network simulation: for most net-
d’(T) = q(T) = p’(T) = 0, a)(T) = - 1 (12) works the tableau matrix 1) is sparse and 2) has widely
separate eigenvalues (i.e., it describes a stiff system of ODE).
* ag(f)/&z, etc., denotes a/&z{g(a(f), q(t), p, t)}, where a(t), q(t) is the These assumptions dictate that for efficiency, implicit
solution of (8). quasi-A-stable integration formulas be used [7]. The back-
HACHTEL et al. : SPARSE TABLEAU APPROACH TO ANALYSIS AND DESIGN 105

ward differentiation formulas (5) of order k(1 I k< 6) tive starting value. Experience shows that this prediction
satisfy these requirements. Using the Nordsieck vector 1121, step, if efficiently executed, is usually worthwhile. The kth-
Gear [8] has demonstrated their advantage in providing a order prediction formula
variable-order variable-time-step approach to minimizing
kfl
the number of time steps subject to specified allowable x;+1 E - izl Yixn+l-i (18)
integration errors. In the present context of simulation and
optimization such an approach is especially necessary be-
has the same order of accuracy as (17). Again, the yi are
cause the computation of &D/ap (cf., discussion of Theorem
chosen to annihilate polynomials of degree <k. The coeffi-
2) requires storage of all arguments of the nonlinearities at
cients yi are also used in computing the truncation error
every time step.
Gear’s basic strategy, with minor modifications, of step h(G n+l)true - (-%+l)amxoximate) 3 Ek t o(hk+2) (I91
size and order selection based on prediction and control of
truncation error have been incorporated into the tableau where CC,+lhpproximatedenotes the value given by (14). E, is
obtained by expanding the terms x,+ r-i, i = 0, 1, . . . , k, in
approach. However, the Nordsieck vector has been replaced
(14) in a Taylor series around t,, 1. Using the fact that (18)
with stored backward differences to reduce operations
and (14) annihilate polynomials of degree I k, one can show
counts when integrating a large number of equations. This Al.,+
LlldL
replacement gives rise to new difference formulas for pre-
diction and truncation error, and these will be presented E = H hkf ’ dCk+‘)x,+ 1
along with the description of the discretization of the d/dt k. kk+l dtk+’
operator. Since this operator can act on true charges (q), k+l
the objector (a), or parametors (p), we state the develop- + C Yixn+l-i w-9
ment in terms of a generic vector x. Thus restated (5) reads i=l

where
(14) 1 k tn+l-i- tn+i k+i
Hk =j&gocri h >
where x,+l-i~x(t,+l-i), h;zti+,-ti, and h=h,. Note the
hi depend on i in a variable-step method. The coefficients ai =L fi ttn+l - tn+l-i)

are determined by the requirement that (14) annihilate k! i=l (tn+l - tn)
polynomials of degree Sk. In particular, this means that
(14) must be satisfied for each of the k + 1 substitutions and where the extra coefficient 6, is determined by requiring
xi=(tn+l-tn+l-i)j,j=O, l,**., k,thus yielding k+l equa- the right-hand equation of (20) to be satisfied identically for
tions in the k-t- 1 unknowns CZ~,CI~,. . *, elk These can be polynomials of degree k+ 1. As in Gear’s work on auto-
solved in closed form, yielding an efficient method for matic step size and order control [7], the new time step h and
computing the [Link] in practice one solves for the vector order k are chosen so that h is maximum, consistent with
Axi instead of xi, higher precision is obtained by expressing
(14) in terms of the Axi. To this end, let (21)
X n+l-i E x,ei + Ax,vi, i = 0, 1, . . ..k where T is the total integration time and ET a specified total
and absolute. error. In practice, to reduce roundoff error, the
prediction and truncation error formulas (18) and (20) are
rewritten in terms of the AXi of (16), just as (17) was obtained
(15)
from (14).
Then (14) becomes The Lemma presented demonstrates that it is unnecessary
to predict any components of the w vector (cf., (7)) other
k-l
than those which cause the tableau matrix Fdl to vary with w.
izo P&n-i + hi-,+1 = 0. (16) In the Lemma u,+ 1 represents the subset of w,+ r, on which
depend the nonlinear components of the tableau error
Thus in (7) the effect of d/dt operating on q in forming fd is
vector f(w, t). (For the example of Fig. 4, u = co1 (ug, u2, is)
given by
=d (w12, Wll, WI&
da,+1 - $%I.+1- 4.1- ky$’B&n-i*
-= (17) Lemma
dt I 1
Let w%: I be the result of predicting (cf., (18)) all com-
‘Correspondingly, the d/dt elements in Fl are replaced by ponents of w,+ i, and WC: I be the result of predicting only
-PO/h to obtain Fdl (note PO= CI,,). the subset u,,+~ of w,+ i. Then w::(i) (obtained from the
In (7), one has the option of starting the Newton iteration first Newton iteration step of (7) with wiy 1= wf: J is equal
with wlpl i = w,, or of predicting from past history an alterna- to w::(j) (obtained from (7) with wlp! I = wff r).
106 IEEE TRANSACTIONS ON CIRCUIT THEORY, JANUARY 1971

Thus in contrast to Gear’s method [7], it is not necessary TABLE I

to predict all quantities operated upon by d/dt. Since u is
generally a very small subset of w (3 components of 20 in Type
Type 2
1 elements which
elements which
are +
are - 1
1 called topological type

Fig. (4) the computational advantage of predicting u instead Type 3 elements which never change, called c type
of predicting either 1) all the variables, or 2) all the non- Type 4 elements which change with design parameters, called p type
linearities plus all the charges [13] is substantial. The proof Type 5 elements which change with time, called t type
Type 6 elements which change with the unknown, called x type
of the Lemma, which will be omitted here, is based on the
fact that the tableau matrix is identical for the two modes of
prediction, i.e., a p-type element should properly be executed in the time
loop of’(7) but not in the inner Newton loop. To this end
there is associated with each multiplication of the term
In summary, three difference operators are used, one each S[k, L]*S[K, L], a variability type
for kc-,+1 (as in (17)) xIp! 1 (as in (18)), and ZZ,#Z~+~/(~+ 1))
.(d(k+UX ,,+ l/dtk+l) (as in (20)). Each of these formulas use VTM = max [VT[k, L], VT[K, L]]
stored backward differences Ax,,+iFi. The coefficients cli, which identifies when the multiplication must be done.
yi, 6,, are computed at each time step. The time step and In order that OPTORD may reduce the total operations
order (k, 1 I k I 6) are determined as recommended by Gear count in a given design problem, a set of weights wt[i],
using the inequality (21). Prediction is only done for the i= 1, . . ,6, are specified, one for each of the six variability
variables U, + I which cause Fdl to vary. types. Thus the weight wt [i] can be an estimate of relatively
IV. SPARSE MATRIXMETHODS how often i-type operations must be executed during an
optimization run.
There are basically two sparse matrix programs used.
The permutation vectors rp and cp are used to specify the
The first program OPTORD deals with finding an optimal
pivoting order. Thus the ith pivot is in row rp[i] E Z and
pivot order, i.e., an order of elimination of rows and col-
column cp [i] = J. The task of OPTORD is to select Z and J so
umns so that “fill-in” and “operations count” are mini-
as to minimize the ith weighted pivot cost
mized at the same time numerical stability is retained. Quali-
tatively, in terms of reducing fill-in, the program OPTORD PC[Z, J] = pM[Z, J] + vR[Z, J] (22)
uses a strategy due to Markowitz [4]. Similar techniques
where ,u,v 2 0 are specified input parameters, M is a weighted
can be found in Tinney-Walker [3] and Dantzig et al. [4].
multiplication count, and R is a roundoff error factor. We
OPTORD is novel, however, in its introduction of variability
define
type to identify and avoid redundant operations in cyclic
computations. The use of a numerical test in OPTORD is M[Z, J] = c wtj (23)
also not usually incorporated in these methods. j=l

The second program 1-2-j GNSO takes the permutation with wcj standing for the weight associated with the highest
generated by OPTORD and, taking into account the variabil- variability type involved in thejth of the m,, multiplications
ity types of the matrix elements, writes five machine code required in eliminating row Z and column J. Also,
programs (called C-SOLVE, T-SOLVE, X-SOLVE, B-SOLVE, and
A-SOLVE) which execute only the nontrivial arithmetic 1 C NM
operations necessary in the Crout factorization algorithm (24)
for Gaussian elimination [15]. This program is based on a
R[zyJl = ij js[z,J]l
similar program GNSO [l 1, which generated only a single where j runs over the rI nonzero column indices in row I.
FORTRAN program, and did not recognize variability type. The ith of the n formally identical steps in the process of
We discuss in turn OPTORLI and t-2-3 GNSO, both of which selecting rp and cp is described by the following algorithm.
recognize three basic properties of the nz nonzero elements
OPTORD Algorithm
of the ns x ns matrix S, namely, 1) location Z, .Z; 2) numerical
value S[Z, J] ; and 3) variability type VZ’[Z, J]. The six At the end of the (i- 1) step, i- 1 eliminations have been
variability types are listed in Table I. carried out. Thus in the ith step the (n + 1 - i) x (n + 1 - i)
These types are represented in the example of the matrix matrix S@)is to be considered (note S(‘)=S) in which the
FL of Fig. 6. In this example, the 1, 5 and 5, 11 elements of Fl nonzero locations, the values, and the variability types of
are the topological type (I/T[l, 5]=1, I/T[5, 11]=2), the its elements are known.
12, 5 element is p type (I/T[12, 5]=4), the 17, 17 element is 1) For each row K determine the subset Z[i],j= 1, . . . , i,
t type (because dldt is replaced in Fdl by -cl,/h, which of the rK nonzero columns of the Kth row of SC’)for which
changes with time in a variable-step method, I/T[17, 171
T 0 L. R[K, l[j]] 51 (25)
= 5), and element 15, 7 is of x type (I/T[15, 7]= 6). Note
that in a time-domain optimization, topological-type ele- where 0 I T 0 L I 1 is a threshold factor. The remaining
ments are invariant, p-type elements change every optimiza- rK - P, elements in row K have numerical values which are
tion step, t-type elements every time step, and x-type ele- considered too small to be candidate pivots.
ments every Newton iteration (cf. (7)). 2) Determine the P, values of M[K, lb]] (cf., (23)), j= 1,
Ideally, an arithmetic operation is done only when and 3,, associated with the elimination of row K and column
if it must be done. Thus multinlication of a t-type element bv 2Til (element K, lfil is the candidate pivot). In Gaussian
HACHTEL et al. : SPARSE TABLEAU APPROACH TO ANALYSIS AND DESIGN 107

elimination a multiple of the pivot row K is added to each The choice of minimizing multiplications was in keeping
row k in Sci),which has a nonzero element in column lb], in with minimizing the length of the SOLVE code and maximiz-
such a way that the k, Z[i] element becomes zero. Thus for ing executi’on speed.
each such row k we compute l-2-3 GNSO
S(‘)[k, L] - S(‘)[k, l[j]]*(S@‘[K, L]]/S”‘[K, l[j]]) (26) l-2-3 GNSO creates an executable nonlooping program
for solving both b=Sx=LUx=Ly and 6’=1z’S-FLU
where L runs over the union of column indices in rows k and =yU. 1-2-sGNSO differs from GNSO [l ] mainly in its use of
K. Associate with each nontrivial multiplication (repre- variability type in creating five distinct lists of machine in-
sented by * in (26)) its variability type structions : C-SOLVE ; T-SOLVE ; X-SOLVE ; B-SOLVE ; and
A-soLvn (GNSO creates a single list of FORTRAN instructions).
VTM[k, L]
C-SOLVE, T-SOLVE, and X-SOLVE compute the factorization of
=
- max {J’T[k, lb]], JJ’T[K, l[j]], VT[K, L]} S into the product of upper and lower triangular matrices L
0, if S”‘[K, L] = 0. (27) and U. C-SOLVE, T-SOLVE, and X-SOLVE need be executed only
once per optimization step, time step, and Newton step,
If there are cj nonzero elements in the pivot column lb] of
respectively (cf. (7)). B-SOLVE carries out the back substitu-
Sti), then (cj- l)(rK- 1) such multiplications would be re-
tions Ly= b, Ux=y required in simulation (cf., (7)).
quired to use element K, lb] as the next pivot. The weighted
A-SOLVE carries out the adjoint back substitution 6’ = jj’U,
multiplication count for each of the 0, pivots K, Z[i],
jY = YL required in the sensitivity computation (11) or (13).
j= 1,2,. * . ) BK, is therefore given by3 The logic of l-2-3 GNSOis based on the Crout method of
M[Kl[j]] = wt[VT[K,$]]] factorizing S into LU. For purposes of presentation here it
will be assumed that S has been rearranged so that the
w@‘TM[k
L-j] (28) Pivots S[rpKl, cp[K]] appear on the diagonal of S. Thus
> the Crout formulas are for K= 1, 2, . . . , n,
where k runs over the nonzero row indices in column Z[i],
X-l
and L runs over the union of column indices in rows k and
Z = K;..,n
K. The first term accounts for the computation of LIK = SIK - c LIJUJK ’
J=l ,>
(S”‘[K, Z[i]])-’ and for k=K; the sum accounts for the K-l
multiplications (S”‘[K, l[i]])-‘*S(“[K, L].
&I = &I - c
3) Compute PC[K, lb]], j= 1, . . . , P,, using (24) and find .I=1
the minimum pivot cost column JC[K] for each row Z = K + 1,. . , n. (31)
K, K = 1, 2, * * . , n + 1 - i, i.e.,
PC[K, JC[K]] = mjn {PC[K $I]>. (29) (Note U,,= 1.) The index ordering indicates how the Kth
column of L and Kth row of U are to be computed se-
4) Find the minimum operations-count row I, i.e., quently and assures that all right-hand-side terms of (31)
are known.
M[Z, JC[Z]] = m$ {M[K, JC[K]]}. (30) As in OPTORD the variability type of each new L,, or U,,
term is equated to the highest type of any term in (31).
5) Execute the numerical Gaussian elimination of row Similarly, each product in (3 1) is given the highest type of its
ZE rp(i) and column J= JC(Z)? q(i), creating the (n-i) two factors. These variability types are used to decide to
x (n - i) matrix S@+I). which lists, C-SOLVE, T-SOLVE, or X-SOLVE, the instructions
6) Update the variability types of each nonzero element
created to compute each L,, (or U,,) are to be added. For a
in Sci+ ‘) (cf., (26)), according to
multiplication, this is done by placing the required machine
VT[k, L] = max { VTM[K, L], VT[k, L]}. instructions (LOAD, MULTIPLY, ADD) at the end of the
appropriate list. For example, instructions for x-type
Note that if the operation count weight is set to zero
multiplications are added to X-SOLVE. If a multiplication
(i.e., ZJ= 0 in (22)), step 3 selects the largest pivot element in
involves a topological type factor (+ l), the LOAD and MULTI-
row K. Thus partial pivoting [17] is effected in each row,
PLY instructions are suppressed. If a number of topological
and the row is chosen (cf., (30)) to minimize the weighted
type terms are to be added, the result must be topological
multiplication count. The decision to have OPTORD mini-
(f 1 or 0) (zero if the number of such terms is even).” In this
mize only the number of multiplications required at each
‘step of Gaussian elimination was made on the basis of some
4 That the sum of two or more terms of topological type in a Crout
comparison with other methods. Roughly, the result of this formula is also of topological type ( f 1 or 0) is due to the fact that the
comparison was that it was important to base the next pivot tableau with all nontopological type elements replaced by zero is a totally
choice on the result of executing the previous pivot. How- unimodular matrix [18]. This is becausein the tableau we have made the
restriction that not more than one topological element may appear in a
ever, how one then chose the best pivot was not so critical branch-constitutive-relation row. Therefore, since the incidence matrix
and many methods seemed to give the same benefits. These A is totally unimodular and any totally unimodular matrix augmented by
conclusions are in agreement with those reported in [15]. rows or columns of the identity matrix is totally unimodular, it follows
that the topological part of the tableau matrix is also totally unimodular.
The result then follows from the fact that the elements of L and U can be
3 We define wt [0] = 0 since V TM can be zero (cf., (27)). expressedin terms of subdeterminants of the tableau matrix.
108 IEEE TRANSACTIONS ON CIRCUIT THEORY, JANUARY 1971

case the only instruction generated is a LOAD of & 1 (or no vector s is used to store the nonzero nontopological values
instruction if the result is 0). in the (row-wise) order in which they appear in E;. A vector
If a given formula (31) involves multiplications of various SP, which points to s, and a column index vector CZ are as-
types, we have sociated with the nonzero elements of l$. Finally, a row
pointer vector RP is used to identify which column indices
L,, (or UK,) = TC(IC) + TT(ZT) + TX (32) belong to a given row. For the example of Figs. 1-5, the
where TC(IC) is a temporary storage location for the sum arrays and the components of SP, CI, and RP, belonging to
of terms of f 1 p type or c type and TT(IT) a temporary rows 11-17 of E; (cf., Fig. 5), are

SP(21-36) = (1, 3, 2, 4, 5, 6, 2, 1, 2, 7, 2, 8, 1, 2, 1, 9)
t t t tt tt
C1(21-36) = (11, 6, 17, 20, 12, 13, 18, 14, 19, 7, 9, 20, 16, 20, 12, 17)
t t t tt tt
RP(ll-17) = (21,22, 25, 28, 30, 33, 35, 1. (33)
location for the sum of t-type terms. TX represents the sum The row pointers point to the first nonzero column index
of the remaining terms, but no temporary storage location in a given row. Thus RP and CI combine to define a unique
is required. The instructions required for computing I, J location in Fl for each component of the pointer SP. This
TC(ZC), TT(IT), and TX are placed at the end of the pointer in turn identifies the component of s to be associated
C-SOLVE,T-SOLVE, and X-SOLVE lists, respectively. The de- with (F,),.
tailed procedure for TT[ZT] (the procedure is similar for 1-z-3 GNSO requires knowledge of SP, CI, and RP, as well
TC and TX) is given in the following list. We add the follow- as a similarly compacted variability-type array VT. Never-
ing instructions to the T-SOLVE list. theless, the SOLVEcode generated refers only to the com-
pacted s array (which thus replaces the S elements in (31)).
1) LOAD, ADD S,, (if S,, is t type and nonzero) Thus, although in the present example, location, value, and
2) LOAD, MULTIPLY, ADD (for each t-type product in (3 1)) variability type must be known for 41 of the 400 elements in
3) ADD TC(ZT) (if present in (32)) F,, only 10 quantities need be in core during the simulation-
4) STORE L,, (if TX absent) optimization phase of NDP. Note also (cf. Fig. 4) that in
5) STORETT(ZT) (if TX present). (33) only 4 quantities, v3, i,, v2, and i,, need be stored on
A development similar to that based on (3 1) and (32) de- disk for the purpose of recreating s (i.e., E;) for the back-
scribes the creation of the instruction lists for the back solve ward time adjoint integration (11).
programs B-SOLVE (for solving Sx=b) and A-SOLVE (for In the NDP language the input description of the time-
solving R’S=&). Like X-SOLVE, B-SOLVE (or A-SOLVE) is domain network optimization problem of Figs. l-5 is
executed within the Newton iteration since the right-hand- shown in Table II. The first data card defines the reference
side vector b and the auxiliary vector y change whenever the node (to be node 0 of Fig. 3). The next nb cards, one for
unknown variable x changes. each branch in the network, are free format with successive
Thus by creating the five lists of instructions to be exe- data fields giving the branch label, the start node, the finish
cuted in the appropriate loops of the computation, l-2-3 node, the element type, and the name and argument list of a
GNSO achieves near-ultimate efficiency : multiplications are user defined function subprogram. The fifth card gives the
performed only if they are absolutely necessary. name of the subprogram which defines the objective func-
tion to appear as integrand in (1) and signals NDP to add
V. AN INPUT LANGUAGE FORTHE TABLEAU APPROACH the objector and parametor to the network. The sixth and
An automated network design program (NDP) has been seventh cards state that there is one designable parameter
written, incorporating the approach to automated network and that in the constrained optimization problem (1) the
optimization previously described. NDP interfaces the parameter labeled 1 has a lower bound of 0.01, a starting
sparse matrix, implicit integration, and adjoint sensitivity value of 1, and an upper bound of 100. The eighth card
routines through a simple SCEPTRE-like [9] input language. identifies the unknown variables required to recreate the
As described in the following paragraphs, the input lan- matrix E; in (11). According to the Lemma of Section III,
guage processor passes a complete description of the those are also the variables which are predicted in the for-
tableau matrix 4, which is identified as the S matrix of ward time integration (7). Note NDP translates bracketed
OPTORDand l-2-3 ~~~0.1-2-3 GNSOthen produces a machine function arguments into unknown vector components. In
code program, SOLVE(cf., Section IV). To conserve storage, the example, /!‘2+i,-tx(6), /A3/+v,+x(13), and /Pl/
the SOLVEcode is designed to work only on the nonzero -‘pl +x(20). Other examples of notation are Nl (voltage of
nontopological entries of E;. node 1) and Q3 (charge on the capacitor in branch 3).
Before describing the input language in detail, the Given the card input of Table II, NDP creates a sub-
threaded list structure used to compact E; is described. A routine VARQ and input to two permanent subroutines
HACHTEL et d. : SPARSE TABLEAU APPROACH TO ANALYSIS AND DESIGN 109

TABLE II varying inductor is then computed and stored in QL. Then

the derivative aQL/c?p, (dummy argument S2) is computed
NETWORK DESCRIPTION, NREF= 0
1, l-0, 3 = E(t) in terms of aL/ap, (dummy argument S3) and stored in S(4).
2, 1-2, 12 = QL(IT2IPlIt) The third executable program QTVARQ creates the Ea-4
3,2X)), 11 = PSI (lA2/A3/) (cf. (8)) rows of the tableau vector f and the d/dt elements of
OPTIMIZE PHI (/73/Pl/)
DESIGN PARAMETERS, NP = 1 the tableau matrix. QTVARQ computes Ea-4 by replacing 4 by
Pl (0.01, 1, 100) a summation as described by (17).
NONLINEARITIES /A3/A2/T3/T2/ Note that FTOP, VARQ, and QTVARQ operate directly on
the input data of Table II with a minimum of intervening
data processing. There are no instructions in NDP which
TABLE III refer to excess or link capacitors or tree inductors. Thus the
tableau approach imposes no topological restrictions what-
f(l1) = x(11) - FTIME(t) soever on the type or connectivity of the network elements.
fV2) = QW(f9, s(3), x(201, s(4), t) -x(17) Note also that functions like those in Table IV are entirely
f(13) = PWx(l2), 4% x(13), ~(6))
f(l5) = PHW(7), s(7), x(201, ~(8)) -x(9) user defined. Any function of any group of variables may
thus be defined by the user so that there is no restriction on
the type of elements treatable by the tableau approach.

TABLE IV VI. RELATIONBETWBENPIVOTINGORDERAND

NETWORK METHODSOFANALYSIS
FUNCTION QL(X1, il, X2, S2, T)
Sl = FL(X2, s3, q When the time discretization procedure (cf., Section III)
QL = Xl * Sl (zi,*L) is carried out before any reduction is made on the network
equations, the only constraint on pivoting order is that the
aL
s2 = Xl * S3

RETURN
( >
=i2*-
aPI
numerical value of the pivot element should not be too
small. From the point of view of pivot order selection, we
examine two of the principal methods of network analysis :
END
nodal analysis and state-variable analysis. For ease of
illustration consider the tableau of Fig. 9 which describes a
FTOP and QTVARQ. These subroutines compute the tableau network of only capacitances and other admittance-type
error vector f and the compacted version s of the tableau elements but no forcing terms. Using the pivoting order
matrix E;. The input to FTOP is the (suitably compacted) (4,6), (2,4), (3,5), (5, 3), (6,2) we obtain the nodal equations
incidence matrix A. Note there is a + 1 (- 1) in the column
of A corresponding to the start (finish) node. Thus A is essen- A;YA,+A;CA,; (34)
tially read directly from the input cards. The result of FTOP
is the evaluation of the first nn + nb components off.
The generated subprogram VARQ for the present example where the matrix operating on node voltages I/ is known as
comprises the list of FORTRAN instructions shown in Table the nodal admittance matrix.
III. There is a one-to-one correspondence between these The state equations are more difficult to describe in terms
FORTRAN instructions an4 input cards 2-5 of Table II, except of a pivoting order since it requires finer groupings to be
that the argument lists in the user-defined function sub- made on the variables. Thus a proper tree must be obtained,5
programs, QL, PSI, and PHI, differ from those on the input the tree variables separated from the link variables, and all
cards. This difference is for user convenience. The conven- variables except the capacitor voltages in the proper tree
tion is that given the unknown vector x the subprograms and the inductor currents in its links eliminated. For a
must return not only the function value (e.g., QL) needed simple example, it is now assumed that the capacitor
for computingf, but also the partial derivatives (needed for branches of Fig. 9 form a tree ; hence this is a proper tree.
evaluating 4) of the function value with respect to members Thus for the tableau of Fig. 9, the pivoting order (2, l),
of the argument list which are components of the unknown
vector x, e.g., s(3) = dQL/ax(6) = L, s(4) = 8QL/ax(20) s ' To illustrate that in the general casethe state-variable method corre-
sponds to a particular set of pivots and this set gives rise to a proper tree,
iz(8L/8p,). Thus unknown vector components in the argu- consider the following constraints on the pivot order.
ment list are always followed by a component of the com- 1) Select the first nb (number of branches) pivots from A V- Iv=0
pacted tableau matrix array s. For user convenience these rows of Fd,(Fig. 6). If possible, let the kth pivot fall in the first nonzero
column in the kth row of A. In this case, which will occur nn (number of
extra arguments are left out of the input cards. independent nodes) times, the kth branch is in the tree. If after k- 1
The user-specified FORTRAN function subprogram QL is pivots, the kth row of A has been annihilated, then pivot on
described as an example in Table IV. In this function sub- Ikk=(FdJ”“+k,nn+“b+[Link] this case the kth branch is a link.
2) Select the next nn pivots to fall in locations in the A’i=O rows which
program, another FORTRAN function subprogram FL is called are the transpose of the nn pivot locations in A.
to evaluate the inductance L-FL (dummy argument SI) These constraints causea tree to be selectedwhose identity dependson
as a function of the time T and the design parameter p1 the ord& of branches on the input cards. A proper tree [lo] can be ob-
tained by rearranging the tableau so that voltage sources come first in the
(dummy argument X2). The value of L is stored in S(3) list of branches followed by capacitors, conductors, resistors, inductors,
(cf., Table III). The magnetic flux for this linear time- and current sources.
110 IEEE TRANSACTIONS ON CIRCUIT THEORY, JANUARY 1911

tions on the input data and therefore are not precise. Let E
-
be the error in S and Ebe the error in b2. Thus in the machine
AT, AT2 the reduced problem is represented by
-I
9 (S + E)y = b”, + 1. (38)
I! !-I 1 To solve (36) or (38) precisely (to the accuracy of the

:;I
m1
machine) one can use iterative refinement [17]. For (36) this
means the iteration

Fig. 9. Simplified tableau.

(3,5),(1,2), (4,6), ($3) yields the equations Y(k+ 1) = yW + Ay’W.

For (38) this iteration is
c% = - (A;)-lA;YA,A;'?l, (35) (S + E)Ayck’ = & + E - (S + E)ytk’
(k+ 1) = y’k’ + Ay’W.
Y
which are the state equations in this case.
In both nodal and state-variable analysis the equations According to Wilkinson [17], in the iterative refinement of
are reduced to their standard form and then the dldt the problem Mz= b, if the matrix M being inverted can be
operator is replaced by a difference operator by some pro- factored into LU, where LU= M(I+F), with
cedure analogous to that described in Section III. The dif-
ference equations thus obtained are identical to the differ- IFIlm 5 3 (39)
ence equations obtained by first replacing the d/dt operator
and if the error vector b - Mzck’ is computed with accumula-
and then pivoting in the same sequence as was done to
tion, then [Iz(~)-- z*II m <2-’ for k large enough. (Here z* is
obtain the state-variable or nodal-variable equations.
the true solution and t is the number ofdigits in the machine.)
Both cases (34) and (35) represent a reduced form of the
Thus to within the accuracy of the machine and subject to
original tableau matrix FL. Thus the classical network
(39), equations (36) and (38) can be solved exactly. In the
methods may be regarded as special cases of the tableau
case of equation (36) we then have an accurate answer.
approach. From this point of view the particular methods
However, in the case of (38) we only have an accurate answer
are identified by the constraints they impose on the pivoting
to equation (38) but not necessarily equation (36) which was
order in Gaussian elimination. Since these constraints on
our original purpose. Thus the initial errors E and Emade in
the pivoting order are not in general consistent with objec-
manipulating the input data to form the reduced equations
tives of minimizing operations count, it follows that these
cannot be recovered.
classical network methods, being primarily motivated by
analytical simplicity, are intrinsically slower than the sparse VII. NUMERICAL RESULTS FOR A PRACTICALEXAMPLE
tableau approach.6 The sparse tableau approach also has
The current switch-emitter follower logic chain of Fig.
an advantage in terms of the possibility for higher precision
,lO is used to demonstrate optimal design and large-scale
answers. To illustrate this, consider the following system of
application capabilities of the tableau approach and to
linear equations
compare its computational efficiency with competitive
methods. Each five-terminal stored model in Fig. 10(c)
(2I)(;)=(ii). (36) represents one logic stage (Fig. 10(b)) in the chain. Each of
these logic circuit models use three transistors, which in turn
are represented by the high-speed transistor model of
Let this represent the unreduced tableau and the system Fig. 10(a).
The automated design problem selected has eight logic
(D - CA-lB)y = b, - CA-lb, (37) stages and is characterized by the parameter set p’ = (R,, R,,
R,) (for the lower and upper parameter bounds we use
represent the reduced equations, e.g., these would corre-
pk-(0.05, 0.05, 0.05) and ph-(0.15, 0.15, 0.15)). The design
spond to state-variable or nodal analysis. Since the unre-
objective is to realize a specified output voltage waveform
duced tableau is taken directly from the input data, the
Qt(t) for a given input voltage waveform Q(t) (Fig. 11). The
matrices A, B, C, D and vectors b,, b, can be assumed to be
specified output is characterized by delay, rise time, and
precise. The matrix S-D- CA- ’ B and the the vector
up and down levels. Thus the objective instantaneous func-
b2 = b,- CA-’ b, are the result of some arithmetic opera-
tion 4 = 3( I/realized - I/specified)’ vanishes when the shaded
area of Fig. 11 vanishes. For an eight-circuit logic chain the
6 It can be seenthat the symbolic elimination of topological equations
is formally equivalent to the l-2-3 GNSO process of converting multiplica- result of six Fletcher-Powell optimization steps (bottom of
tions by + 1 into additions or subtractions. Fig. 11) is that Q(p) decreases from 0.7 to 0.0003.
HACHTEL et cd. : SPARSE TABLEAU APPROACH TO ANALYSIS AND DESIGN

I ‘RISE 1

m(p) RI R2 R3

INITIAL .7 75n loon 150a

FINAL .3x10-3 II451 loon Illa

E RG Fig. 11. Three-parameter optimization of eight-stage logic chain

TABLE V
(b)
Problem complexity (8 circuit optimization)
nodes 180
“OUT branches 460
nonlinear elements 112
. .. differential equations 192
total equations 1300
... Optimization (Fletcher-Powell with quadratic search)
... designable parameters 3
gradient computations (&D/lap) 6
function evaluations (@(JI)) 36
Fig. 10. Current switch-emitter follower circuit logic chain. (a) Two- time steps 30 000
dimensional transistor model. (b) Current switch-emitter follower newton iterations 85 000
logic stage. (c) Logic chain.

Various aspects of problem complexity are given in so the original tableau has only 446 nonzero elements) and
Table 4. Note that for this example, p-type elements serve to quantify the comparison of the tableau approach
(variability type = 4) of the tableau matrix F1varied 36 times using unconstrained OPTORD with OPTORD constrained to
and t- and x-type elements varied 30 000 and 85 000 times, correspond to competitive network analysis methods given
respectively. in Section VI. In each case the ordering parameters T 0 L, p,
To show how code length varies with problem size, we v, and wt[i], i= 1, 2, . . . , 6, (cf., (22)-(26)) were chosen by
give in Table V the results of applying OPTORD/I-2-3 GNSO trial and error to minimize code length irrespective of ac-
to logic chains of 1, 8, and 22 circuits. curacy considerations and are tabulated in Table VIII.
For the data given in Table VI we made no attempt to try For the nodal-analysis method, the last nn ‘(number of
to minimize the length of the code but chose partial pivoting independent nodes) pivots were constrained to come from
(p= 0, v= 1 in (22)) to insure accurate answers. The data the upper-left-hand partition of Fl (originally zero in Figs. 6
indicate that for the one-circuit case, I;; has 462 nonzero and 9). For the state-variable method, the last ns (number of
elements, of which 92 are stored in the threaded list (33). state variables is ns = 21) pivots were constrained to be in the
The L/U factorization adds another 450 fill-ins so that there state-variable columns and in the corresponding d/dt rows.
are 912 nonzero elements in L/U, of which 360 must be Within these constraints, the full OPTORD algorithm was
stored during execution of the generated SOLVE code. The used to select the pivots.
sum of the lengths of the data, and of the c-, T-, x-, B-, and One basis for comparing the performance of the gen-
A-SOLVE codes is 27 316 bytes. The central processing unit erated codes is the sum of the X-SOLVE and B-SOLVE code
IBM Model 91 times required for function evaluation lengths7 since these are executed in the innermost loop of the
(time-domain simulation of 24 ODE) and gradient compu- program (85 000 times in the example of Table V). In Table
tation are 5 s and 3 s, respectively. Note that the various VII the data (compare cases 1 and 2 to 5 and 6) demonstrate
data in Table VI (in particular note the total code length) our general experience that the pivot order constraint im-
increase approximately linearly with the number of circuits posed by the nodal-analysis method did not significantly in-
analyzed and, therefore, with the order of l$. (Such linear crease SOLVE code lengths over unconstrained OPTORD. The
growth may be typical for network tableau matrices.) How- pivot order constraint (cases 3 and 4) imposed by the state-
ever, the growth of simulation time (@D(p)column of Table variable analysis method caused a more substantial increase
VI) is supralinear. This is because both the number of time in code length. In all cases observed to date, unconstrained
steps and the number of Newton iterations per time step OPTORD has given the best performance, with nodal analysis
generally increase with problem size. a fairly close second. (An explanation of this phenomenon
The results of Table VII illustrate the effects of varying will be given in a subsequent paper.)
the ordering strategy to minimize the length of x-type code
(X-SOLVE and B-SOLVE). The data are for the one-circuit case ’ Since the generated codes are relatively homogeneous and have no
of Fig. IO(b) (this time without objectors and parametors loops, the code length is roughly proportional to execution time.
112 IEEE TRANSACTIONS ON CIRCUIT THEORY,JANUARY 1971

TABLE VI
PARTIAL PIVOTING (p=O, U=~)OPTORD/I-z-3 GNSO FOR LOGIC CHAIN OF 1 CIRCUIT, 8 CIRCUITS,
~~CIRCUITS(PARAMETORSAND OBJECTORINCLUDED)

4 -Vu Model 91 Time(s)

Total Code Number of
Nodes/Branches Nonzero Stored Fill Stored in Bytes WP) amp ODE

16/61 462 92 450 360 27 316 5 3 24

180/460 3472 680 3735 3430 276 404 120 30 192
49211260 8632 1736 9332 8617 713 504 900 500 928

TABLE VII TABLE VIII

COMPARISONOF SPARSETABLEAU APPROACH(NDP) ORDERINGPARAMETERSFORCASESOFTABLEVII
WITH CLASSICAL NETWORK METHODS
OPTORD/I-?-3GNS0 Input
OPTORD/I-2-3 GNSOOUtpUt
Case TOL /I u ~(1) ~(2) ~(3) wt(4) ~(5) it
LIU Code Length (bytes)”
Analysis 1 0 0.5 1 100 100 500 400 450 1500
Case Mode Fill Stored T-SOL X-SOL B-SOL 2 0 1 0 100 100 150 200 150 300
3 0 0.5 1 100 100 500 400 450 1500
1 Nodal 215 4 0 1 1 100 100 500 600 450 900
-0PTORD 5 0 1 0.5 100 100 90 101 90 150
6 0 1 0 100 100 102 101 102 150
2 Nodal 194
-0PTORD
time invariant rows of the tableau vectorf(e.g., the first 11
3 S-V 472 rows in Fig. (4), the right-hand side would stay zero if started
-0PTORD
at zero. This is not the case. Instead it has been found neces-
4 s-v 500 sary to evaluate the full right-hand side at every Newton
-0PToRD
iteration to obtain rapid convergence.
5 NDP The preceding results and discussion, combined with the
191
-0PTORD error discussion of Section VI, suggest that nodal or state-
NDP
variable constraints on the pivoting order offer no appreci-
6 204
-0PTORD able storage or execution speed advantages and create the
possibility of accuracy problems.
a Upper number in bracketed pairs represents output of current
1-1-3GNSO compiler. Lower number represents estimates based on hypo- VIII. CONCLUSIONS
thetical optimum compiler.
The advantages claimed for the tableau approach are
Note in Table VII a correlation between generated code 1) generality, 2) simplicity, 3) speed, and 4) problem size.
lengths and the number of stored quantities in L/U. This Generality is claimed because any system of differential
suggests that minimizing the fill in of non f 1 quantities and algebraic equations can be handled. In particular, for
may be the best criterion for optimizing code lengths. electrical networks no topological restrictions are necessary
When circuits become much larger, as in the example of and there are no restrictions on the type of elements allowed.
Table V, accuracy requirements rather than code length This allows, for example, the inclusion of the objective
can become the dominating factor in selecting ordering function in the tableau by introducing the objector element.
parameters. The nodal and state-variable column con- Simplicity is achieved with the tableau approach by
straints on the pivoting order prevent the implementation of basing the main computational task of the program on
a true partial pivoting strategy (p=O, v= 1) and thus con- Gaussian elimination. All other computations have a simple
tribute to a reduction in the accuracy of the answers. This direct relation to the input data for the program. Although
consideration may be the cause of problems observed in the programs OPTORD and 1-2-3 GNSO are not simple, they
NODAL OPTORD code for large values of the time step h and are general in nature and relate to the more general problem
with s-v OPTORD for small values of h. Another accuracy of solving fix, a)=0 by Newton’s method for different
consideration is the fact that the number of Newton itera- values c1where LJf/ax is sparse and has elements with different
tions per time step increases with problem size for large net- variability types. It therefore seems wise to develop such
works. The evidence suggests that the Newton iteration is programs to a high degree and base more special programs
not only solving the nonlinear difference equations in (7), (such as network design programs) around these develop-
but also serving to refine the solutions of the linear equa- ments.
tions. This consideration has thus far prevented successful The claim of speed is based on the fact that operations
runs on very large networks with parameters set for other count in the basic Newton loop has been minimized so
than partial pivoting. This aspect of the accuracy problem that only necessary operations are executed. This has been
may also prevent extensive exploitation of possible sparsity done by using variability type and near-optimal ordering
of the right-hand-side vector. One might think that for linear algorithms in OPTORD. In addition, the number of time steps
HACHTEL et al. : SPARSETABLEAU APPROACH TO ANALYSIS AND DESIGN 113

required for integration has been made small by using im- If transient analyses for even larger networks are desired,
plicit methods and Gear’s error control algorithms. the l-2-3 GNSO SOLVE code could be buffered into core. This
Large problems can be handled under this program be- would save a factor of 2-3 in storage, but cause an unknown
cause storage is saved by not storing O’s and 1’s. In addition, degradation in execution speed. An alternative approach
OPTORD saves space by causing little fill in. Since the SOLVE would be to not generate the SOLVE code but to use a scheme
code contains no loops, it can be stored in auxiliary storage similar to that suggested by Chang [9]. This would save a
and executed as it is brought into core thus requiring no factor of 3 in storage, and perhaps a factor of 2 in speed
storage. (However, for maximum speed, it should remain would be lost. (However, with this method, variability type
in core.) cannot be exploited.) Thus a 7500 branch network transient
As was shown in Section VII, the s-v OPTORD and NODAL- analysis would be possible and a rough lower bound on
OPTORD methods offer no advantages over the uncon- execution time would be about 4 h on the 360/91.
strained OPTORD methods. In fact, possibly serious ac-
curacy problems may occur. Note that the NODAL-• PTORD ACKNOWLEDGMENT
column constraints cause differencing errors in the back The authors would like to thank T. E. Grapes, who was
substitution which may be serious when arguments of non- responsible for programming the compiler which converts
linearities are branch voltages. Thus if such a branch voltage GNSO to l-2-3 GNSO, and D. A. Calahan, for extensive dis-
is small in comparison with its terminal node voltages, it will cussions and the interchange of ideas.
be computed to fewer significant digits.
Note that some speed and storage advantages may accrue REFERENCES
to the nodal or state-variable methods if the elimination [l] F. G. Gustavson,W. Liniger, and R. A. Willoughby, “Symbolic
leading to the appropriate final variable set is done sym- generation of an optimal Crout algorithm for sparse systems of
linear equations,” J. Ass. Comput. Mach., vol. 17, pp. 87-109,
bolically at the input processor level. Since variables which January 1970.
enter linearly and are neither energy storage nor involved [2] R. K. Brayton, F. G. Gustavson, and R. A. Willoyghby, “Some
in output can be eliminated, the amount of B-SOLVEcode can results on sparsematrices,” to be published in Math. Comput.
[3] W. F. Tinney and J. W. Walker, “Direct solutions of sparsenetwork
be reduced. However, the excess B-SOLVE code (which equations by optimally ordered triangular factorization,” Proc.
corresponds roughly to solving Ax= (6, -By) in (36)) IEEE, vol. 55, pp. 1801-1809, November 1967.
could, with some programming effort, be removed from the 141 Proc. Sparse Matrix Symp., IBM Rep. RA-1, March 1969.
i5j G. D. Hachtel and R. k. Rohrer, “Techniques for the optimal design
output of l-2-3 [Link] this excess code is a deficiency of and synthesis of switching circuits,” Proc. IEEE, vol. 55, pp. 1864
the implementation but not of the sparse tableau approach. 1877, November 1967.
Note that a priori symbolic reduction introduces topolog- [6] F. H. Branin, Jr., “Computer methods of network analysis,” Proc.
IEEE, vol. 55, pp. 1787-1801,November 1967.
ical matrix-vector operations and complicates the forming [7] C. W. Gear, “The automatic integration of stiff O.D.E.‘s,” Proc.
of the Jacobian matrix and right-hand side vector. (These 1968 IFIPS Congr., pp. A81-A85.
additional operations correspond roughly to the forming of [S] -, “The control of parameters in the automatic integration of
ordinary differential equations, ” J. Ass. Comput. Mach.; Department
D-AC’B and b2- CAplb, in equation (37).) In con- of Computer Science, University of Illinois, Urbana, Internal Rep.,
trast, the analogous operations for OPTORD analysis are p. 14, May 1968.
executed by the highly efficient SOLVE code. [91 S. Sedore et al., “SCEPTRE-an automated digital computer pro-
gram for determining the responseof electronic systemsto transient
In NDP we have chosen to emphasize speed at some cost nuclear radiation.” vol. 11. IBM Snace Guidance Center, Oswego,
to problem size by ‘keeping the SOLVE code in core. On a N. Y., Tech. Rep: AFWL-TR-66-li6.
large machine such as 360191 with 2 x lo6 bytes of memory, 1101E. S. Kuh and R. A. Rohrer, “The state-variable approach to network
analvsis. Proc. IEEE, vol. 53, pp. 672486, July 1965.
this turns out to be a good choice. On this machine the
[Ill R. Fletcher and M. J. D. Pdwell, “A rapidly convergent descent
capacity limitation of NDP seems to be about 2500 method for minimization.” Comaut. J., vol. 6, pp. 163-168, 1965.
branches. A transient analysis of a network with 1200 [I21 A. Nordsieck, “On num&ical integration of o&nary differential
equations,” Math. Comput. vol. 16, pp. 22-49, 1962.
branches took about 20 min. An entire optimization in-
u31 D. A. Calahan, “Optimization of switching circuits,” presented at
volving three parameters of a 460 branch network and re- the 2nd Biennial Cornell Conf. Engineering Applications of Elec-
quiring 36 function evaluations and 6 gradient evaluations tronic Phenomena, 1969.
took 1 h. Extrapolating linearly, we calculate that optimiza- [I41 H. M. Markowitz. “The elimination form of the inverse and its
application to linear programming,” Management Sci., vol. 3, pp.
tion for a 2500 branch network would require 7 h. However, 255-269, April 1957.
this is low since larger networks usually require more time u51 G. B. Dantzig, R. P. Harvey, R. D. McKnight, and S. S. Smith,
“Sparse matrix techniques in two mathematical programming
steps per transient integration and more Newton iterations codes,” Proc. Sparse Matrix Symp., IBM Rep. RA-1, pp. 85-99,
per time step. Further, for a large circuit more than three March 1969.
parameters will normally be used which surely will increase [I61 P. D. Grout, “A short method for evaluating determinants and
solving svstemsof linear equations with real or complex coefficients,”
the number of function and gradient evaluations required. AIEE-Tians. (Supplementj, vol. 60, pp. 1235-1240,-December1941.
Thus all indications are that 7 h is a very rough lower bound [I71 J. H. Wilkinson, The Algebraic Eigenvulue Problem. New York:
for optimization of a 2500 branch network. It would not be Oxford University Press, 1965.
VI __ Rounding Errors in Algebraic Processes. Englewood Cliffs,
surprising if as many as 20-60 h were required. On the N. J.‘: Prentice-Hall, 1963, pp. 121-126.
other hand, no serious attempt has been made yet to op- [19] S. Seshu, “Topological considerations in the design ot driving pomt
timize the various parameter settings in NDP to cut down functions,” Proc. AMS, pp. 1068~1073,October 1965.
[20] A. Chang, “Applications of sparsematrix methods in electric power
the number of time steps, Newton iterations, etc., so that system analysis,” Proc. Sparse Matrix Symp., IBM Rep. RA-1,
this figure may be cut down. pp. 113-121, March 1969.