1
PRIMA: Passive Reduced-order Interconnect Macromodeling Algorithm
*
Altan Odabasioglu, Mustafa Celik
+
, and Lawrence T. Pileggi
Department of Electrical and Computer Engineering
Carnegie Mellon University
Pittsburgh, PA 15213
Abstract
This paper describes PRIMA, an algorithm for generating
provably passive reduced order N-port models for RLC
interconnect circuits. It is demonstrated that, in addition to
requiring macromodel stability, macromodel passivity is
needed to guarantee the overall circuit stability once the
active and passive driver/load models are connected.
PRIMA extends the block Arnoldi technique to include guar-
anteed passivity. Moreover, it is empirically observed that
the accuracy is superior to existing block Arnoldi methods.
While the same passivity extension is not possible for MPVL,
we observed comparable accuracy in the frequency domain
for all examples considered. Additionally, a path tracing
algorithm is used to calculate the reduced order macromodel
with the utmost efficiency for generalized RLC interconnects.
1. Introduction
As integrated circuits and systems are designed with
smaller feature sizes and for faster operation, RLC intercon-
nect effects have a more dominant impact on signal propaga-
tion than ever before. In addition, parasitic coupling effects
and reduced power supply voltage levels make interconnect
modeling increasingly important. Since these interconnect
models can contain thousands of tightly coupled R-L-C com-
ponents, reduced order macromodels are imperative
[1][2][3][4]. Ideally, a simulator would isolate the large lin-
ear portions of the circuit from the nonlinear elements (e.g.,
transistor models) and preprocess them into reduced order
multiport macromodels.
It is well known that an N-port can be fully represented by
its admittance parameters in the Laplace domain, however,
the objective is to apply model order reduction to produce
low order rational approximations for each entry in Y(s) (see
Fig.1). To find Y(s), voltage sources are connected to the
ports and the currents into the ports are measured. The volt-
age sources are the inputs to the system and the port currents
are the outputs. A single-input single-output (SISO) N-port
model approach would perform model order reduction on
+
He is now with Motorola, Austin, TX 78721.
*This work was supported in part by the Defense Advanced Research
Projects Agency, sponsored by the Air Force Office of Scientific Research
under grant number F49620-96-1-0069, and by grants from Intel Corpora-
tion and SGS-Thomson Microelectronics.
each term Y
ij
individually. Both Asymptotic Waveform Eval-
uation (AWE) [1] and Padé via Lanczos (PVL) [2], which
are Padé approximations, can perform SISO reduction by
matching 2q moments for a q
th
order approximation of each
Y
ij
term. The Arnoldi Algorithm [4] can also be used to
obtain SISO approximations, however it matches only q
moments for a q
th
order approximation. MPVL (Matrix Padé
via Lanczos) [5] and Block Arnoldi [6] are multi-input
multi-output (MIMO) versions of PVL and Arnoldi respec-
tively. In the block techniques, the system Modified Nodal
Analysis (MNA) matrices are directly reduced by matrix
transformations.
Regardless of the reduction method used in all of the
approaches cited above, the reduced order model of an RLC
circuit can have unstable poles. It is always possible to
obtain an asymptotically stable model by simply discarding
the unstable poles, however, passivity is not guaranteed. In
addition, discarding unstable poles requires re-adjustment of
the residues to improve the quality of the approximation.
Passivity uncertainty is problematic since even the test for N-
port passivity can be very costly for a large number of ports
[7]. The coordinate transformed Arnoldi Algorithm [8] was
introduced as a remedy for the instability problem, but it
does not guarantee passivity. The passivity extension of this
stable Arnoldi algorithm was recently developed in [9], how-
ever its applicability is limited to RC circuits only. The
PACT algorithm [3] proposed a new direction for passive
reduced-order model for RC circuits based on congruence
transformations. The same authors proposed Split Congru-
ence Transformations [10] for passive reductions of RLC cir-
cuits, producing equivalent circuit realizations. In [10],
however, the extra steps required to split the transformation
matrix can result in a decrease in accuracy and efficiency.
Moreover, the passivity proof is somewhat controversial, and
we will consider a more complete proof in this paper.
Linear circuit with
N ports, defined by N
voltage-current
pairs (v
k
, i
k
)
i
1
i
N
....
+
-
+
-
v
1
v
N
I
1
s()
I
N
s()
Y
11
s()…Y
1N
s()
Y
N1
s()…Y
NN
s()
V
1
s()
V
N
s()
=
...
...
...
...
Y(s)
FIGURE 1: The multiport representation of a linear circuit.
0-89791-993-9/97 $10.00 1997 IEEE
2
A passive system denotes a system that is incapable of
generating energy, and hence one that can only absorb
energy from the sources used to excite it [11]. As we will
show in Section 2.2., passivity is an important property to
satisfy because stable, but not passive macromodels can pro-
duce unstable systems when connected to other stable, even
passive, loads. A property in classical circuit theory states
that: interconnections of stable systems may not necessarily
be stable; but (strictly) passive circuits are (asymptotically)
stable; and arbitrary interconnections of (strictly) passive cir-
cuits are (strictly) passive, and, therefore, (asymptotically)
stable [12].
In this paper, we propose a Passive Reduced-order Inter-
connect Macromodeling Algorithm, PRIMA, based on the
Block Arnoldi Algorithm but with congruence transforma-
tions that produce provably passive reduced order macro-
models for arbitrary RLC circuits. PRIMA has accuracy
comparable to MPVL and superior to Block Arnoldi. Fur-
thermore, the block Arnoldi vectors are generated with the
utmost efficiency following the algorithms in RICE [13] that
are used to calculate moments. This includes efficient han-
dling of interconnect trees and meshes, as in RICE, but with
renewed focus on efficient handling of large problems with a
huge number of mutual inductances.
2. Background
To obtain the admittance matrix of a multiport, voltage
sources are connected to the ports. The multiport, along with
these sources, constitutes the Modified Nodal Analysis
(MNA) equations:
(1)
The i
p
and u
p
vectors denote the port currents and voltages re-
spectively and
(2)
where v and i are the MNA variables corresponding to the
node voltages, inductor and voltage source currents respec-
tively. The n x n matrices G and C represent the conductance
and susceptance matrices (except that the rows corresponding
to the current variables are negated as in [7]). N, Q and H are
the matrices containing the stamps for resistors, capacitors
and inductors respectively. E consists of ones, minus ones
and zeros, which represent the current variables in KCL equa-
tions. Provided that the original N-port is composed of pas-
sive linear elements only, Q, H and N are symmetric
nonnegative definite matrices. This implies C is also sym-
metric and nonnegative definite. Since this is an N-port for-
mulation, whereby the only sources are the voltage sources at
the N port nodes, B=L. But we maintain the separate B and L
notation for generality of the equations.
Cx
n
˙
Gx
n
– Bu
p
+=
i
p
L
T
x
n
=
G
NE
E
T
– 0
≡ C
Q0
0H
≡ x
n
v
i
≡
Returning to equation (1), following the notation in [2]
we define
and . (3)
With unit voltages at the ports, taking the Laplace transfor-
mation of (1) and solving for the port current variables, the
y-parameter matrix is given as
(4)
where I
n
is the n x n identity matrix. It is apparent from (4)
that the eigenvalues of A represent the inverses of the poles
of Y(s).
Using any of the aforementioned model-order reduction
techniques, we can find reduced order rational approxima-
tions to Y
jk
(s) terms, for all j, k ≤ N. The reduced-order
Y(s) can then be simulated along with other nonlinear and
linear portions of the complete circuit using a simulator that
employs either recursive convolution [14] or state-space
realization [7], both of which have linear complexity. If the
reduction is block, the reduced order multi-input multi-out-
put circuit can also be realized using linear circuit elements.
2.1.Block Arnoldi Algorithm
The Block Arnoldi algorithm reduces the system matrix
A in (3) to a small block upper Hessenberg matrix H
q
. To
do so requires an orthonormal basis, X, for the correspond-
ing Krylov space which satisfies the following:*
(5)
where N is the number of ports and I
q
is a q x q identity ma-
trix. The Krylov space is defined as
. (6)
Finding the reduced order admittance matrix can be
explained by a change of variable,
. (7)
where z
q
is now the reduced order system variable, which re-
duces the number of unknowns in the system (q is generally
much smaller than n). Substituting (7) into (1), then multi-
plying the first equation by yields
(8)
Therefore, in the Laplace domain,
(9)
where I
q
is a q x q identity matrix.
The reduced order system equations and admittance
matrix are given by (8) and (9) respectively. The poles of
*The operator is the truncation to the nearest integer towards zero.
AG
1–
C–≡ RG
1–
B≡
Ys() L
T
I
n
sA–()
1–
R=
.
colsp X() Kr AR
q
N
----
,,
=
X
T
AX H
q
=
X
T
XI
q
=
Kr ARk,,()colsp RARA
2
R…A
k
R,, ,,[]≡
x
n
X
nxq{}
z
q
=
X
T
G
1–
H
q
z
q
˙
z
q
X
T
Ru
p
–=
i
p
L
T
Xz
q
=
Y
ˆ
s() L
T
XI
q
sH
q
–()
1–
X
T
R=
3
the reduced order system are the reciprocal eigenvalues of
H
q
. A complete pole/residue decomposition can be obtained
by eigendecomposing H
q
. Using the information in [6], it
can be shown that the first moments of in (9)
match those of in (4).
2.2.Importance of Passivity
It is always possible to come up with stable reduced order
macromodels by utilizing a number of heuristics, however,
none of these tricks can be used to obtain provably passive
approximations. Moreover, in [7] it was shown that the test
for passivity of an AWE-reduced N-port macromodel is pro-
hibitive in terms of CPU run time cost. Fig.2 is a numerical
example generated in [7] that demonstrates the passivity
problem. Y
1
(s) in this figure represents a reduced order
transfer function which has all poles and zeros in the left half
plane. Y
dr
(s) represents a capacitor and resistor in parallel. If
we drive this circuit with a current source, it will oscillate at
2.5/π Hz. To show that it is also a practical problem, we took
a simple interconnect and connected the load and the nonlin-
ear driver as shown in the examples section in Fig.5. The
interconnect is represented by a fifth order approximation
obtained by PVL [2]. The figure clearly shows the growing
oscillations at the output (instability) although all of the
poles obtained from PVL were stable. A Thevenin equivalent
linear driver with a resistance of 2 ohms generates a similar
instability for this 2-port example.
3. PRIMA: Passive Reduced-order
Interconnect Macromodeling Algorithm
The Block Arnoldi Algorithm is employed in PRIMA to
generate the orthonormal basis for a congruence transforma-
tion matrix. After (the extra step is not necessary
when is an integer) iterations of PRIMA, the n x q matrix
X is found such that:
(10)
q
N
----
Y
ˆ
s()
Ys()
Y
1
s()
s 4+
s
2
2s 5++
--------------------------=
Y
dr
s() 0.06 0.056s+=
Y
dr
(s)
Y
1
(s)
Z
c
s()
s
2
2s 5++
0.056s
3
0.172s
2
1.4s 4.3+++
--------------------------------------------------------------------------
⇒=
Poles at: 3.074 , 5 j±–
FIGURE 2: A non-passive system example demonstrating
potential instability.
q
N
----
1+
q
N
----
colsp X() Kr AR
q
N
----
,,
=
X
T
XI
q
=
In the classical Arnoldi approach [4], the reduced order
Y(s) is calculated using the qxq upper Hessenberg matrix in
(5) as shown in (9):
(11)
In our variation, the conductance and susceptance matrices
are directly reduced so that passivity is preserved during re-
duction. Applying the change of variable in
(1), and multiplying the first row by from (10) yields
(12)
The reduced order MNA matrices are, therefore,
(13)
❍ Connect voltage sources to the multiport & obtain the
MNA matrices as in (2).
❍ Set [b
1
| b
2
| ... | b
p
] = B and [l
1
| l
2
| ... | l
p
] = L
❍ Solve for R
❍ ; qr factorization of R
❍ Set
❍ For k=1, 2, ..., n
Set
Solve for
For j=1, .., k
; qr factorization of
❍ Set and truncate X so that it has q
columns only.
❍ Compute ,
❍ Find eigendecomposition of : *
❍ To find poles and residues for :
Solve for w
*
Set and
❍ Set
GR B=
X
0
K,()qr R()=
nint
q
N
----
1+=
VCX
k1–
=
GX
k
0()
V= X
k
0()
HX
kj–
T
X
k
j1–()
=
X
k
j()
X
k
j 1–()
X
kj–
H–=
X
k
K,()qr X
k
k()
()= X
k
k()
X
X
0
X
1
…X
k 1–
=
C
˜
X
T
CX= G
˜
X
T
GX=
G
˜
1–
C
˜
G
˜
1–
C
˜
SΛS
1–
=
Λ diag λ
1
λ
2
…λ
q
,,,()=
Y
ˆ
ij,
s()
G
˜
wX
T
b
j
=
µS
T
X
T
l
i
= νS
1–
w=
Y
ˆ
ij,
s()
µ
n
ν
n
1 sλ
n
+
-----------------
n 1=
q
∑
=
Y
ˆ
s()
Y
ˆ
11,
…Y
ˆ
1p,
Y
ˆ
p1,
…Y
ˆ
pp,
=
...
...
*Inversion of can be avoided.G
˜
FIGURE 3: The Passive Reduction Algorithm
Y
ˆ
s() L
T
XI sH–()
1–
X
T
R=
x
n
X
nxq{}
x
˜
q
=
X
T
X
T
CX()x
˜
˙
q
X
T
GX()x
˜
˙
q
– X
T
B()u
p
+=
i
p
L
T
X()x
˜
˙
q
=
C
˜
X
T
CX= G
˜
X
T
GX=
B
˜
X
T
B= L
˜
X
T
L=
4
These types of transformations are known as congruence
transformations. Congruence transformations were first in-
troduced by [3] for order reduction of circuits. From (12) and
(13), the reduced Y(s), namely , is now
(14)
Since the size of and is typically very small, it is
easy to find the poles and zeros of by eigendecomposi-
tion. The complete algorithm is given in Fig.3. It employs
the Block Arnoldi Algorithm using modified Gram-Schmidt
orthogonalization [6], which is mathematically equivalent to
ordinary Gram-Schmidt process, but behaves better numeri-
cally [15]. In addition, it is possible to avoid the inversion of
to find the poles and residues by using a generalized
eigendecomposition. In this case, the computation of
can be avoided by using (31) and replacing it by . It is
observed that this scheme is numerically much better.
The complexity of the algorithm to produce q poles for an
N-port is slightly less than AWE, PVL and MPVL. It
requires 1 LU factorization (or path tracing equivalent as
explained in Section 4.) of the G (MNA conductance)
matrix, which dominates all the other computational costs
and is common in all reduction techniques. However, to find
q poles, only q backward-forward substitutions are needed,
whereas in MPVL, PVL and AWE, twice as many are
required. As in MPVL, there will be only one eigendecom-
position to find the poles and residues, whereas PVL requires
N
2
eigendecompositions, since for each Y
ij
(s), there will be a
different T
q
. AWE will solve the N
2
different Hankel matri-
ces to get to the poles.
3.1.Preservation of Passivity
If the system described by (1) and (2) is reduced by the
transformations in (13), it can be shown that the reduced sys-
tem is always passive. In [16], necessary and sufficient con-
ditions for the system admittance matrix (eqn. (14)) to
be passive are:
1. for all complex s, where
*
is the complex
conjugate operator.
2. is a positive matrix, that is
for all complex values of s satisfying Re(s) > 0 and for any
complex vector z.
The second condition also implies the analyticity of for
Re(s) > 0, since is a rational function of s (details in
[16]). Therefore, the test of analyticity is unnecessary.
Due to the fact that the reduced matrices, , , and
are all real since the transformation matrix, X, is real, condi-
tion 1 is automatically satisfied. To show that condition 2 is
Y
ˆ
s()
Y
ˆ
s() L
˜
T
G
˜
sC
˜
+()
1–
B
˜
=
G
˜
C
˜
Y
ˆ
s()
G
˜
G
˜
1–
B
˜
X
T
R
Y
ˆ
s()
Y
ˆ
s
∗
() Y
ˆ
∗
s()=
Y
ˆ
s() z
∗
T
Y
ˆ
s() Y
ˆ
∗
T
s()+()z0≥
Y
ˆ
s()
Y
ˆ
s()
G
˜
C
˜
B
˜
L
˜
satisfied, we first set and use the
property (since B=L in our formulation when there
are no sources inside the N-port, ) and some
algebra to obtain,
(15)
Setting and yields,
(16)
Similarly, let to get
(17)
Since C is symmetric, . C is known to be non-
negative definite (since we negate the rows corresponding to
current variables as in (2)), so
(18)
for any complex vector y and σ=Re(s) > 0. N (the resistor
stamps) is a symmetric nonnegative definite matrix, there-
fore
(19)
is also nonnegative definite for any complex vector y. From
(17), (18) and (19), it follows that the second passivity con-
dition is satisfied.
3.2.Preservation of Moments
The transformation in (13) preserves moments of the
original system, which is the same as the classical Block
Arnoldi reduction and half of that in MPVL. The proof is as
follows. The exact (block) moments, M
i
, of the circuit are
given as:
(20)
where , and G, C, B, L are the system
matrices as defined in (1).
Likewise, the moments of the PRIMA reduced order sys-
tem are given by
(21)
where , and , , , are as de-
fined in (13). Substitution of (13) in (21) yields:
Y
h
s() Y
ˆ
s() Y
ˆ
T
s
∗
()+=
B
˜
L
˜
=
X
T
BX
T
L=
z
∗
T
Y
h
s()zz
∗
T
B
˜
T
G
˜
sC
˜
+()
1–
B
˜
B
˜
T
G
˜
s
∗
C
˜
+()
T–
B
˜
+()z=
z
∗
T
B
˜
T
G
˜
sC
˜
+()
1–
G
˜
sC
˜
+()G
˜
s
∗
C
˜
+()
T
+[]G
˜
s
∗
C
˜
+()
T–
B
˜
z=
wG
˜
s
∗
C
˜
+()
T–
B
˜
z= sjωσ+=
z
∗
T
Y
h
s()zw
∗
T
G
˜
σjω+()C
˜
+()G
˜
σjω–()C
˜
+()
T
+[]w=
w
∗
T
G
˜
G
˜
T
σC
˜
C
˜
T
+()++[]w=
w
∗
T
X
T
GG
T
σCC
T
+()++[]Xw=
yXw=
z
∗
T
Y
h
s()zy
∗
T
GG
T
σCC
T
+()++[]y=
C
T
C+2C=
y
∗
T
σC
T
C+()y2σy
∗
T
C y 0≥=
y
∗
T
G
T
G+()yy
∗
T
NE
E
T
– 0
T
NE
E
T
– 0
+
y=
y
∗
T
2N0
00
y 0≥=
q
N
----
M
i
L
T
A
i
R=
AG
1–
C–≡ RG
1–
B≡
M
ˆ
i
L
˜
T
A
˜
i
R
˜
=
A
˜
G
˜
1–
C
˜
–≡ R
˜
G
˜
1–
B
˜
≡ G
˜
C
˜
B
˜
L
˜
5
(22)
It is shown in [6] that the Arnoldi algorithm yields
. (23)
Rearranging the terms and using the definitions from (13):
(24)
(25)
Inserting (23) in (25) results in:
(26)
where
. (27)
From (26), it can be shown by recursion that
(28)
Therefore, using (27) it follows that
(29)
Replacing in (22) with
yields
(30)
Evaluating (26) when i=0 gives
(31)
Then from (30) and (31),
, (32)
Finally, combining (32) and (28) with (20), it follows that
(33)
Note that the number of poles in each entry of Y(s) is q,
and we have matched the first moments at all N ports,
yielding a total of q moments. The number of moments
matched in PRIMA is, therefore, the same as that for the
Block Arnoldi algorithm and half as many as matched by
MPVL.
4. Integration of PRIMA within RICE
For all of the model order reduction schemes, the LU
decomposition of the MNA conductance matrix (G in (2))
dominates the run time. In [13], RICE (Rapid Interconnect
Circuit Evaluation) was described as a general path tracing
M
ˆ
i
L
T
XX
T
G X()–
1–
X
T
CX()[]
i
X
T
G X()
1–
X
T
B=
A
i
RXH
q
i
X
T
R=,0i
q
N
----
<≤
AA
i 1–
RXH
q
i
X
T
R=
G–
1–
CA
i 1–
RXH
q
i
X
T
R=
C–A
i1–
RGXH
q
i
X
T
R=
X
T
–CA
i 1–
RX
T
GX H
q
i
X
T
R=
XX
T
GX()
1–
X
T
CA
i 1–
R– XH
q
i
X
T
R=
KA
i 1–
RA
i
R,=0i
q
N
----
<≤
KXX
T
GX()
1–
X
T
C–=
K
i
RA
i
R,=0i
q
N
----
<≤
XX
T
G X()–
1–
X
T
CX()[]
i
K
i
X=
XX
T
G X()–
1–
X
T
CX()[]
i
K
i
X
M
ˆ
i
L
T
K
i
XX
T
G X()
1–
X
T
B=
XX
T
GX()
1–
X
T
BR=
M
ˆ
i
L
T
K
i
R=,0i
q
N
----
<≤
M
ˆ
i
M
i
,=0i
q
N
----
<≤
q
N
----
algorithm to obtain moments with optimal efficiency for
interconnect trees and mesh structures. Using RICE to cal-
culate moments, the explicit construction and inversion of
G is avoided, and the moments are more accurate than those
obtained via matrix factorization.
The moments of the circuit can be obtained recursively
from:
(34)
where the matrices G, C and B are as defined in (2). As
shown in [1], this can be viewed as recursive dc circuit solu-
tions, when capacitors and inductors are replaced by current
and voltage sources respectively, with the values derived
from the columns of . The Krylov vectors, which
can be viewed as well conditioned moments, can be obtained
from a very similar recursive scheme:
The “orth” operator can be implemented as a simple Gram-
Schmidt orthonormalization procedure. The space spanned
by the block Krylov terms is called the
Krylov space. Therefore, the Krylov vectors can be obtained
via a path tracing procedure using RICE-like routines to
solve for equations (35) and (37).
The Krylov space constitutes the congruence transforma-
tion matrix, X in PRIMA. The reduced MNA matrices
and are
. (40)
Note, however, that the matrices and are obtained
using RICE without explicitly constructing G and C. The
columns of are the values of current and voltage sourc-
es that are used to replace capacitors and inductors at each
moment computation stage. This information is easily ob-
tained during a path trace [13]. The k
th
block of (i.e.
) is a function of previous blocks of and
since from (38),
(41)
and using ,
M
0
G
1–
B=
M
k
G
1–
CM
k 1–
= k 0>
CM
k 1–
1. Obtain zero
th
moment and orthonormalize it:
Solve from (35)
(36)
2. Recursively obtain higher order Krylov vectors:
Solve from (37)
(38)
(39)
M
0
GM
0
B=
X
0
orth M
0
()=
M
k
GM
k
CX
k 1–
=
X
k
ζ
M
k
X
k 1–
X
k 1–
T
M
k
()– …–X
0
X
0
T
M
k
()–=
X
k
orth X
k
ζ
()=
X
k
X
k1–
…X
0
,,,()
G
˜
C
˜
C
˜
X
T
CX= G
˜
X
T
GX=
C X GX
C X
GX
GX
k
GX CX
k 1–
GX
k
ζ
GM
k
GX
k 1–
X
k 1–
T
M
k
()– …–GX
0
X
0
T
M
k
()–=
GM
k
CX
k 1–
=
