OYARO AND TRIVERIO - TURBOMOR

1

TurboMOR: an Efficient Model Order Reduction
Technique for RC Networks with Many Ports
Denis Oyaro, Student Member, IEEE, and Piero Triverio, Member, IEEE

arXiv:1507.00219v1 [cs.CE] 1 Jul 2015

Submitted for publication on the IEEE Transactions on Computer-Aided Design of Integrated Circuits and
Systems on July 1st, 2015

Abstract—Model order reduction (MOR) techniques play a
crucial role in the computer-aided design of modern integrated
circuits, where they are used to reduce the size of parasitic networks. Unfortunately, the efficient reduction of passive networks
with many ports is still an open problem. Existing techniques do
not scale well with the number of ports, and lead to dense reduced
models that burden subsequent simulations. In this paper, we
propose TurboMOR, a novel MOR technique for the efficient reduction of passive RC networks. TurboMOR is based on momentmatching, achieved through efficient congruence transformations
based on Householder reflections. A novel feature of TurboMOR
is the block-diagonal structure of the reduced models, that makes
them more efficient than the dense models produced by existing
techniques. Moreover, the model structure allows for an insightful
interpretation of the reduction process in terms of system theory.
Numerical results show that TurboMOR scales more favourably
than existing techniques in terms of reduction time, simulation
time and memory consumption.
Index Terms—Model order reduction, many ports, moment
matching, parasitics, partitioning.

I. I NTRODUCTION

W

HILE designing VLSI chips, engineers need to take
into account the parasitic resistance, capacitance and
inductance of signal- and power-delivery interconnects, in
order to prevent signal and power integrity issues [1]–[3].
Electromagnetic solvers are used to extract RC or RLC interconnect models, which are then connected to non-linear
devices for system-level simulations. Unfortunately, parasitic
networks can be very large, featuring a huge number of
components, nodes and ports. Direct simulation involving such
large networks is often prohibitive. Model order reduction
(MOR) is frequently used to reduce parasitic models to a
manageable size, and accelerate subsequent simulations.
Several approaches to MOR have been proposed in the last
decades, such as node elimination [4], Krylov subspaces [5],
[6], and balancing [7]. Krylov methods are widely used for
parasitic reduction, since they are more scalable than balancing
methods. Among them, PRIMA [8] is one of the most popular
and widely used Krylov algorithms. PRIMA’s success is due
to its ability to guarantee the passivity of the ROM, a mandatory property to prevent divergent transient simulations [9].
This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada (Discovery Grant program) and in part by the
Canada Research Chairs program.
D. Oyaro and P. Triverio are with the Edward S. Rogers Sr. Department of
Electrical and Computer Engineering, University of Toronto, Toronto, M5S
3G4 Canada (email: piero.triverio@utoronto.ca).

Unfortunately, PRIMA can become very inefficient when
applied to networks with many ports. PRIMA generates the
reduced model through a congruence transformation with an
orthogonal matrix that spans a suitable Krylov subspace. The
orthogonal projection matrix is dense and can become very
large when ports are many. Generating the ROM becomes
very time consuming, since it involves products between large
and dense matrices. In some cases, even storing the projection
matrix can be challenging. Moreover, the obtained reduced
model is dense, large and frequently slower than the original
system. These issues affect most existing techniques and are
an outstanding issue in MOR [10].
A number of techniques have been recently proposed to
address such challenges. Methods like SVDMOR [11], ESVDMOR [12], RECMOR [13] and several others [14], [15] aim
at reducing the number of ports before applying PRIMA. This
is done by exploiting the correlation that may exist between
different ports. However, practical networks with many ports
rarely exhibit a high degree of correlation [16].
In [17]–[19], the problem of reducing networks with many
ports is simplified by clustering inputs into small groups, and
reducing each subsystem individually. These methods generate
accurate and block diagonal ROMs that are sparse. However,
since subsystems are treated independently, passivity is not
always guaranteed.
Another method known as SIP [20] offers a more efficient
approach to moment matching for RC networks. Rather than
explicitly constructing the projection matrix, sparse matrix manipulations are used to generate the reduced matrices directly
using the Schur complement, an idea also used in PACT [21].
This makes SIP more efficient than PRIMA for large networks
with many ports. However, SIP can match only two moments
per expansion point. This level of accuracy is not always
sufficient for practical applications [20], as we will show in
Sec. IV. The authors in [20] suggest using multi-point moment
matching [5], [6], [22] to achieve more accuracy. However, the
obtained reduced matrices can be singular, and avoiding this
issue does not seem to be trivial.
In [23], the SparseRC method is proposed, combining
graph-partitioning techniques [24] with a SIP-like reduction
process. A divide and conquer strategy is used to partition
the original system into smaller subsystems, then reduced
separately with a method similar to SIP [20]. The resulting
ROM has the same partitioned structure as the original system.
Such a reduction strategy is efficient in terms of memory and

OYARO AND TRIVERIO - TURBOMOR

2

cpu time for large networks, since the problem of reducing
the large system simplifies to reducing smaller subsystems
that can be managed efficiently. The generated ROM is also
sparse. SparseRC, however, like SIP, is limited to matching
two moments per expansion point. While PRIMA can be used
to match additional moments, as suggested in [23], this reduces
efficiency, because of the limitations of PRIMA discussed
previously.
In this paper, we propose TurboMOR, a novel MOR technique for RC networks with many ports. TurboMOR achieves
moment-matching without explicitly computing a dense projection matrix as in PRIMA. Efficient and memory-conscious
Householder reflections [25] are used to generate the reduced
model, and match two moments per iteration. Differently from
previous methods such as SIP [20], an arbitrary number of
moments can be matched, providing full control on accuracy.
TurboMOR can be combined with partitioning [23], [24] to
reduce very large networks. A key feature of TurboMOR is the
block-diagonal structure of the reduced models, that addresses
the poor efficiency of the dense models produced by existing
moment-matching techniques. The block diagonal structure
also lends itself to a novel and insightful interpretation of
moment matching in terms of cascaded subsystems. The
reduced models produced by TurboMOR are passive, retain
the input-output structure of the original system, and can be
synthesized into an equivalent RC netlist [26]. Numerical tests
demonstrates the superior scalability of TurboMOR in terms
of reduction time, simulation time, and memory consumption.
The rest of the paper is organized as follows. In Sec. II,
we state the problem and briefly review the foundations
of moment matching. In Sec. III, we discuss the theoretical derivation and practical implementation of TurboMOR.
Sec. IV compares TurboMOR against the state of the art. In
Sec. V we draw our conclusions, and in the Appendix we
provide some mathematical proofs.
II. P ROBLEM F ORMULATION
We consider a passive network made by resistors and capacitors with m nodes and p ports. Using nodal analysis [27],
the network can be described in the Laplace domain by the
systems of equations
(
Gx(s) + sCx(s) = Bu(s)
(1)
y(s) = BT x(s)
where vectors u(s) ∈ Rp and y(s) ∈ Rp collect all port
currents and port voltages, respectively. Vector x(s) ∈ Rm
contains all nodal voltages. Matrices G, C ∈ Rm×m are
conductance and capacitance matrices, respectively. They are
symmetric and non-negative definite. Matrix B ∈ Rm×p maps
input ports to the nodal equations, and T denotes transposition.
The transfer function of (1) reads
H(s) = BT (G + sC)−1 B

(2)

The goal of MOR is to approximate (1) with a model of much
lower order n  m
(
Ĝx̂(s) + sĈx̂(s) = B̂u(s)
(3)
ŷ(s) = B̂T x̂(s)

where Ĝ, Ĉ ∈ Rn×n , B̂ ∈ Rn×p and x̂(s) ∈ Rn . This model
must accurately capture the response of the original system
across the frequency range of interest.
One way of ensuring accuracy is through Padé approximation, also known as moment matching. Around s = 0, the
Taylor series expansion of (2) reads
H(s) = M0 + M1 s + M2 s2 + . . .

(4)

The coefficients Mk are called moments of (1) at DC [5], [6],
[8], and can be related to the systems matrices as
Mk = BT (−G−1 C)k G−1 B ∀k = 0, 1, 2, ...

(5)

The moments of the reduced model are defined similarly, as
the Taylor expansion coefficients of the transfer function
Ĥ(s) = B̂T (Ĝ + sĈ)−1 B̂

(6)

of reduced model (3).
The goal of moment matching is to generate a ROM (3)
that will match the first moments of the original system
Mk = M̂k

∀k = 0, ... , 2q − 1

(7)

up to a given order controlled by q. Since, for RC networks,
moments are typically matched in pairs, we denote the number
of matched moments as 2q. By increasing q the ROM will
become more accurate, but also larger.
In PRIMA, moment matching is performed with a congruence transformation applied to the matrices of the original
system (1)
Ĝ = QT GQ,

Ĉ = QT CQ,

B̂ = QT B

(8)

The columns of Q ∈ Rm×qp span the Krylov subspace
Kq (A, R) = span{R, AR, A2 R, ..., A2q−1 R}

(9)

where A = −G−1 C and R = G−1 B. It can be shown
that ROM (8) matches the first 2q moments of the original
system. The reduced model is of size n = qp, and is passive
by construction since congruence transformation (8) maintains
the non-negative nature of G and C. The projection matrix Q
is constructed numerically with the block Arnoldi process [6],
an orthogonalization procedure similar to the modified GramSchmidt process [25]. Unfortunately, orthogonalization leads
to a dense Q. As a result, when p is high, computing Q and
projection products (8) can be very expensive. For very large
networks, even storing Q becomes an issue, since its size can
easily exceed several Gigabytes. Moreover, transformations (8)
lead to a dense ROM, which will burden any subsequent circuit
simulation. These bottlenecks, which make existing methods
quite inefficient for many-port networks, are tackled by the
proposed method.
III. P ROPOSED M ETHOD
In this section, we discuss the theoretical derivation of
TurboMOR and how it can be implemented for maximum efficiency. The method works recursively, matching two moments
per iteration. We discuss the first two iterations in detail, before
generalizing.

OYARO AND TRIVERIO - TURBOMOR

3

u

x1

(1)

Σ1

The obtained equations lend themselves to a useful interpretation in terms of system theory, depicted in Fig. 1.
System (12a)-(12b) can be seen as the cascade of a system
(1)
Σ1 of order p
(
(1)
(1)
(1)
G11 x1 + sC11 x1 = u1 + B1 u
(1)
Σ1 :
(16)
y = BT1 x1

y

(1)

u1
d
dt

d
dt
(1)

y2

(1)

u2

(1)

Σ2

(1)

Fig. 1. System theory interpretation of (10a) and (10b). The original system
(1)
(1)
has been decomposed into two subsystems Σ1 and Σ2 , decoupled at DC.

and a system Σ2 of order m − p
(
(1)
(1)
(1) (1)
G22 x2 + sC22 x2 = −C21 u2
(1)
Σ2 :
(1)
(1)
(1)
y2 = −(C21 )T x2

(17)

(1)

A. Theoretical Derivation
1) Matching Two Moments: The first iteration of the proposed method is analogous to [20], [21], [23]. Nodes are first
reordered in such a way that port nodes come first, followed
by internal nodes. After reordering, system (1) reads



    
G11
∗
C11
∗
x1
B1
+s
=
u (10a)
G21 G22
C21 C22
x2
0
 

 x1
y = BT1 0
(10b)
x2
where x1 ∈ Rp and x2 ∈ Rm−p denote port and internal
node voltages, respectively. The symbol ∗ is used in symmetric
matrices to denote the transpose of the symmetric block across
the diagonal. For the purpose of shortening our notation, we
do not indicate explicitly the dependency on s for input, output
and state variables. Submatrix G21 describes the resistive couplings present between internal and port nodes. We eliminate
this block through Gaussian elimination, using the congruence
transformation (8) with Q given by


Ip
0
(1)
Q =
(11)
−K−T K−1 G21 Im−p
Matrix K is the Cholesky factor [25] of G22 . For the time
being, we assume G22 to be positive definite (strictly). In
Sec. III-C, we will discuss how a singular G22 can be handled.
Matrix Ip is the identity matrix of size p × p. After the
congruence, equations (10a) and (10b) become
"
#! 
 (1)

  
(1)
x1
C11
∗
B1
G11
0
+s
u
(1) =
(1)
0
0
G22
x2
C21 C22
(12a)


x1
(1)
x2



Only the first subsystem Σ1 is directly connected to the
(1)
input/output ports of the network. Subsystem Σ2 is instead
(1)
(1)
(1)
connected only to Σ1 , through equations u1 = sy2 and
(1)
u2 = sx1 , which define time derivatives. The coupling
between the two subsystems is thus purely dynamical. At DC,
(1)
the second system is completely decoupled from Σ1 and the
network ports, and has no influence on the transfer function
H(s) between input u and output y. At low frequency, the
coupling between the two is weak, and the overall system re(1)
sponse is given mainly by Σ1 . Therefore, the first subsystem
alone can be interpreted as a ROM of order p of the original
system
(
(1)
(1)
G11 x1 + sC11 x1 = B1 u
(18)
ŷ = BT1 x1
In the Appendix, we indeed prove that (18) matches the first
two moments of the original system at s = 0. From an
accuracy standpoint, the proposed ROM is thus equivalent in
size and accuracy to the ROMs generated by other moment
matching techniques. Its computation, however, requires less
effort, since its matrices (13) and (14) can be computed
cheaply using sparse matrix techniques.
2) Matching Four Moments: In order to match more than
(1)
two moments, the presence of Σ2 must be taken into account.
(1)
Instead of applying PRIMA to Σ2 as in [23], loosing
efficiency, we show how additional moments can be efficiently
(1)
matched by further decomposing Σ2 .
First, we apply a congruence transformation to (17) using
Q = K−T in (8)
(
(1)
(1)
(1) (1)
Im−p z2 + sK−1 C22 K−T z2 = −K−1 C21 u2
(1)

(1)

(1)

y2 = −(C21 )T K−T z2

(15)

(19)
(1)
(1)
where x2 = K−T z2 . This step turns G22 into the identity
matrix, and does not require expensive computations since K
is already available from the previous iteration.
Then, with a series of Householder reflections [25], we
compute the QR factorization of the input-to-state matrix
in (19)
 (2) 
R
(2) T
−1 (1)
(Q ) K C21 =
(20)
0

With Gaussian elimination, all resistive couplings between port
nodes and internal nodes have been eliminated, leaving only
capacitive couplings.

where R(2) ∈ Rp×p is upper triangular and Q(2) ∈
R(m−p)×(m−p) is an orthogonal matrix given by the product
of Householder reflectors [25].


y = BT1


0

(12b)

where
(1)

G11 = G11 − GT21 K−T K−1 G21
(1)
C11

(1)
C21

= C11 −
= C21 −

GT21 K−T K−1 C21 − CT21 K−T K−1 G21
+ GT21 K−T K−1 C22 K−T K−1 G21
C22 K−T K−1 G21

(13)
(14)

OYARO AND TRIVERIO - TURBOMOR

4

After Q(2) is applied to (19) with a congruence transformation, the system will read
"
#! "
# 



(2)
(2)
Ip
0
C11
∗
x1
−R(2) (1)
+s
=
u2
(2)
(2)
(2)
0 Im−2p
0
C21 C22
x2
(21a)
#

 "
(2) T x(2)
−R
(1)
1
y2 =
(21b)
(2)
0
x2
where
(2)
C22

(2)


= 0

C21


= 0

(2)
C11


= Ip



Im−2p (Q

(2) T

−1

) K

−T

C22 K

(2)

Q





0

Im−2p

(22)
 
I
Im−2p (Q(2) )T K−1 C22 K−T Q(2) p
(23)
0
 

I
0 (Q(2) )T K−1 C22 K−T Q(2) p
(24)
0

u
(1)

Σ1

x1

y

(1)

u1
d
dt

d
dt
(1)

y2

(1)

(2)

Σ1

(2)

u2

(2)

u1

x1

d
dt

d
dt
(2)

y2

(2)

(2)
Σ2

u2



System (21a)-(21b) is now in the form (12a)-(12b), and the
reduction process used in iteration 1 can be applied again.
System (21a)-(21b) can be seen as the cascade of a first system
(2)
Σ1 of order p
(
(2)
(2) (2)
(2)
(1)
Ip x1 + sC11 x1 = u1 − R(2) u2
(2)
Σ1 :
(25)
(1)
(2)
y2 = −(R(2) )T x1
(2)

and a second system Σ2 of order m − 2p
(
(2)
(2) (2)
(2) (2)
Im−2p x2 + sC22 x2 = −C21 u2
(2)
Σ2 :
(2)
(2)
(2)
y2 = −(C21 )T x2

(26)

The two systems are only dynamically coupled, through equa(2)
(2)
(2)
(2)
tions u1 = sy2 and u2 = sx1 . Overall, the original
system (1) is now decomposed into three blocks, all coupled
dynamically, as shown in Fig. 2. If we retain the first two
(2)
blocks, and neglect Σ2 , we obtained a ROM of order 2p
"
#! 

  
 (1)
(1)
x1
C11
∗
B1
G11 0
+s
u (27a)
(2) =
(2)
0
0
Ip
x1
R(2) C11



 x1
ŷ = BT1 0
(27b)
(2)
x1
As shown in the Appendix, this model matches the first 4
moments of the original system.
3) Matching More Than Four Moments: Additional moments can be matched by iterating the proposed process, and
(2)
further decompose subsystem Σ2 in Fig. 2. This goal can
be achieved by computing, at each iteration j ≥ 3, the QR
decomposition
 (j) 
R
(j−1)
(Q(j) )T C21
=
(29)
0
of the input-to-state matrix of the innermost system (at iter(2)
ation j = 3, matrix C21 in (26)). The QR decomposition
is obtained with a series of Householder reflectors that form
the congruence matrix Q(j) . The obtained system will have
the same structure as (21a)-(21b), and can be seen as the
(j)
cascade of two blocks. The first system Σ1 , of size p, will

Fig. 2. Structure of the system obtained after two iterations of the proposed
method.

add two matched moments to the ROM computed up to that
point. The second system will be further decomposed if j < q.
Otherwise, at the last iteration, it will be discarded. After q
iterations, the obtained ROM will have order pq, and will be
in the form shown in equation (28) at the top of the next
page. In the Appendix, we prove that the obtained model
matches 2q moments of the original network. The proposed
technique therefore leads to a ROM of the same size and
accuracy as PRIMA, but in a more efficient way, which
avoids the explicit construction of a huge and dense projection
matrix. In comparison to SIP [20], that can match only two
moments per frequency point, the proposed method can match
an arbitrary number of moments, and does not suffer from the
singularity issues of multipoint SIP [20]. The use of PRIMA
to match additional moments, advocated in SparseRC [23], is
also avoided.
Another key advantage of the proposed method is the
block-diagonal structure of (28). Unlike PRIMA, that generates dense models, the proposed method naturally leads to
a sparse representation. This reduces the memory footprint
of the ROMs, and accelerates subsequent simulations, as
we shall see in Sec. IV. Although PRIMA models can be
sparsified with an eigenvalue decomposition, this operation
costs extra CPU cycles. The obtained models are stable and
passive by construction, since only congruence transformations
like (8) have been used to generate the ROM matrices. The
positive-definitive nature of G and C in (1) is thus preserved,
which implies passivity and guarantees stable transient simulations [9]. We also note that TurboMOR preserved the matrix
B1 in (10a) and (10b) that maps input ports to state equations.
As discussed in [26], this property facilitates the connection of
the ROM to the surrounding components. Finally, the obtained
ROM can be converted into an RC equivalent circuit using the
procedure in [26], for seamless integration into existing tools
for electronic design automation.
B. Practical Implementation
We now discuss how TurboMOR can be implemented
for maximum efficiency in terms of CPU time and mem-

OYARO AND TRIVERIO - TURBOMOR







5

  x1 

(1)

G11

Ip



(1)



x1 



B1


 

 −R(2) C(2) . . .
 



11
 ..  = 
  ..  + s 
..


.
.




.  
.
..
..
.

∗ 
(q)
(q)
(q)
Ip
x1
x1
−R(q) C11

0
..
.

(2)
  x1

C11

ory consumption. The Cholesky decomposition of G22 can
be obtained using efficient routines for the factorization of
sparse, positive-definitive matrices, such as the supernodal
method [28] available in MATLAB’s chol routine. The QR
decomposition in (20) is computed with the Householder
method. In our MATLAB implementation of TurboMOR, we
used a direct call to the compiled LAPACK routine DGEQRF [29], which returns the orthogonal Q(j) matrix in factored form [25]. Such matrix is never computed explicitly, but
kept in factored form. The LAPACK’s routine DORMQR [29]
can be used to compute products involving Q(j) directly from
(j)
its factorization. Being large and dense, matrix C22 is also
never computed explicitly. Its factored form is always used,
which his given by (22) for j = 2 and by


0
(j)
(j−1)
C22 = [ 0 I(m−jp) ] (Q(j) )T C22 Q(j)
(30)
I(m−jp)
for j > 2.
C. On the Singularity of G22
Throughout the derivation of TurboMOR, we assumed the
block G22 in (10a) to be strictly positive definite, hence
invertible. When this is not the case, we adopt the solution
proposed in [23] for SparseRC. The rows and columns that
make G22 singular are promoted into the first set of equations,
and not eliminated. Since the number of such rows is typically
very low, this does not significantly increase the size of the
obtained ROMs.
D. TurboMOR with partitioning
Graph partitioning techniques can be integrated into TurboMOR to reduce very large networks, such as the power
grid models that we will consider in Sec. IV. A possible
partitioning strategy, used in [23] and [24], is to partition
the given network into subnetworks that interact only through
a limited set of nodes, called separator nodes. An optimal
partitioning can be found with the nested dissection algorithm
nesdid from the SuiteSparse package [28]. Once the network
nodes are reordered according to the partitions identified by
nesdis, the matrices in (1) assume a bordered block diagonal form [30]. To illustrate this, consider a three-component
partitioning of (1)




G1
0
∗
C1
0
∗
 0
G2
∗  + s 0
C2
∗ 
G31 G32 G3
C31 C32 C3
   
x1
B1
× x2  = B2  u (31)
x3
B3

∗

  x(2)
1





u


(28)

0

Blocks G1 , C1 and G2 , C2 correspond to two decoupled
subsystems, that interact only through a set of separator nodes
associated to G3 , via coupling matrices G31 , C31 , G32 ,
C32 . Subsystems 1 and 2 can be reduced individually. The
coupling matrices are then updated accordingly. For instance,
for reducing subsystem 1, we first form its nodal equations



    
G1
∗
C1
∗
x1
B1
+s
=
u
(32)
G31 G3
C31 C3
x3
B3
and then reorder its nodes such that
• port nodes and separator nodes come first, and form the
state vector x1 in (10a);
• internal nodes come second, forming x2 in (10a).
Then, we perform the reduction as in Sec. III-A. After all
subsystems have been reduced, the obtained ROM will read

 
   
Ĝ1
0
∗
Ĉ1
0
∗
x̂1
B̂1
 0
Ĝ2
∗ +s  0
Ĉ2
∗ x̂2 =B̂2  u
x3
B3
Ĝ31 Ĝ32 G̃3
Ĉ31 Ĉ32 C̃3
(33)
As numerical results will show, partitioning reduces the overall
cost of the reduction, since TurboMOR is applied to subsystems of smaller size. Additionally, it reduces the number of fillins in the ROM, since the zero blocks in (31) are maintained
in (33).
IV. N UMERICAL R ESULTS
The proposed TurboMOR algorithm has been implemented
in MATLAB, with direct calls to compiled LAPACK libraries
for a few key operations, namely the QR decomposition
of (20), and the computation of the products with the Householder matrices Q(j) . In this section, we compare the performance of TurboMOR against PRIMA [8] and SparseRC [23].
Computations were performed on a 3.40 GHz Intel i7 CPU,
with 16 GB of memory and MATLAB R2013b.
A. Reduction Time
Table I shows the time needed by the different methods to
reduce various test networks. Example 1 is an on-chip bus
consisting of 128 signal lines. The bus was modelled with
lumped RC segments, and has the characteristics of a global
interconnect in the 65nm technology node [31]. Examples 2 6 are power grid benchmarks obtained from [32]. The original
benchmarks include some inductors, which were neglected. A
variable number of input current sources has been considered
to investigate the scalability of the MOR methods with respect
to port count.
We first compare the proposed method without partitioning
against PRIMA, in order to assess its intrinsic efficiency in

OYARO AND TRIVERIO - TURBOMOR

6

TABLE I
R EDUCTION TIME FOR THE DIFFERENT METHODS ON VARIOUS TEST NETWORKS . A LL TIMES ARE IN SECONDS .
Examples

q

PRIMA
cpu time

SparseRC
cpu time

Speedup w.r.t
PRIMA

TurboMOR

TurboMOR with partitioning

cpu time

Speedup w.r.t
PRIMA

cpu time

Speedup w.r.t
PRIMA

1. On-chip bus
p = 256
m = 38, 528

1

0.94

0.40

2.35×

0.28

3.36×

0.37

2.54×

2

2.84

1.83

1.55×

1.48

1.92×

1.58

1.80×

3

4.78

3.63

1.32×

3.00

1.59×

2.94

1.63×

2. ibmpg1t (RC)
p = 200
m = 25, 195

1

0.41

0.16

2.56×

0.19

2.16×

0.18

2.28×

2

1.10

0.51

2.16×

0.66

1.67×

0.45

2.44×

3

1.98

0.94

2.11×

1.38

1.43×

0.85

2.33×

3. ibmpg2t (RC)
p = 800
m = 163, 697

1

22.00

6.24

3.53×

10.84

2.03×

6.28

3.50×

2

65.55

20.89

3.14×

37.82

1.73×

18.89

3.47×

3

118.64

39.48

3.01×

78.38

1.51×

35.28

3.36×

4. ibmpg2t (RC)
p = 1200
m = 163, 697

1

35.51

9.53

3.73×

16.64

2.13×

9.52

3.73×

2

109.41

32.63

3.35×

60.60

1.81×

29.19

3.75×

3

224.06

64.41

3.48×

132.55

1.69×

56.58

3.96×

5. ibmpg2t (RC)
p = 1500
m = 163, 697

1

49.50

11.66

4.25×

21.29

2.33×

11.76

4.21×

2

152.81

43.06

3.55×

83.66

1.83×

37.58

4.07×

3

729.92

83.76

8.71×

186.95

3.90×

72.32

10.09×

6. ibmpg2t (RC)
p = 2000
m = 163, 697

1

73.17

16.29

4.49×

31.23

2.34×

16.29

4.49×

2

340.18

62.12

5.48×

228.32

1.49×

54.76

6.21×

3

9807.11

122.12

80.31×

1051.72

9.32×

115.71

84.76×

B. Accuracy of the Reduced Models
In this section we demonstrate that, from an accuracy
standpoint, TurboMOR is equivalent to PRIMA. For this purpose, we consider the power grid “ibmpg1t” from [3], which
corresponds to example 2 in Table I. A transient simulation
is performed to calculate the voltage at one of the supply
ports of the power grid, when switching currents are drawn
by the different blocks of the integrated circuit. Fig. 3 shows

Voltage [V]

1.8

1.798

1.796
Original
PRIMA
TurboMOR

1.794
0

2

4

6

8

10

Time [ns]

Fig. 3. Transient response of the original system and the reduced models
obtained with TurboMOR and PRIMA. The reduced models match two
moments (q = 1).
−3

x 10
6
Maximum Absolute Error [V]

matching moments. For each test case, reduced order models
have been generated to match 2, 4, and 6 moments. From the
results in Table I, we observe that TurboMOR is consistently
faster than PRIMA, up to 9.32 times. Savings are particularly
high when order and port count are high, as in example 6.
While PRIMA takes 2 hours and 43 minutes (9807 s) to
match 6 moments, TurboMOR achieves the same result in only
17.5 minutes (1051 s). This speed-up is due to the fact that
TurboMOR achieves moment matching without computing
and storing a large projection matrix as PRIMA does.
Then, we compare TurboMOR with partitioning against the
recently-proposed SparseRC method [23]. From Table I, we
observe that partitioning improves reduction time substantially,
especially for large networks (examples 3, 4, 5 and 6). Comparing the proposed method and SparseRC, we see that for two
moments matched (q = 1), both methods have almost the same
reduction time. This is expected since, in this case, the methods
perform the same operations. However, when additional moments are matched (q = 2 and q = 3), the proposed method
is always faster than SparseRC, which employs PRIMA to
match additional moments, losing some efficiency. This result
shows how, with the Householder transformations proposed in
Sec. III-A, additional moments can be efficiently matched.

PRIMA
TurboMOR

5
4
3
2
1
0
0

2

4

6

8

10

Time [ns]

Fig. 4. Error between the response of the original system and the response
of the reduced models computed with PRIMA and the proposed method. The
reduced models match two moments (q = 1).

OYARO AND TRIVERIO - TURBOMOR

7

150

1.8

TurboMOR
PRIMA
Reduction Time [s]

1.799
Voltage [V]

1.798
1.797
1.796
1.795

Original
PRIMA
TurboMOR

1.794
1.793

100

50

0

0

2

4

6

8

0

10

200

400

Time [ns]

1200

100

−4

x 10

TurboMOR with partitioning
sparseRC

6
PRIMA
TurboMOR

5

Reduction Time [s]

Maximum Absolute Error [V]

1000

Fig. 7. Reduction time for PRIMA and TurboMOR without partitioning vs
number of ports. Both methods match six moments (q = 3).

As in Fig. 3, but for four moments matched (q = 4).

Fig. 5.

600
800
Number of Ports

4
3
2

80
60
40
20

1
0

0
0

2

4

6

8

0

10

200

400

600
800
Number of Ports

1000

1200

Time [ns]

Fig. 6.

As in Fig. 4, but for four moments matched (q = 4).

the time response obtained with the original system and the
reduced models from TurboMOR and PRIMA, for the case of
two moments matched (q = 1). Both methods provide similar
results. This confirms that the proposed method is as accurate
as PRIMA, but more efficient.
In Fig. 4, the maximum error for the two ROMs is depicted.
Figures show that a ROM with only two moments matched is
not suitable for an accurate assessment of the voltage drop
across the power grid. Indeed, the ROMs underestimate the
voltage drop, by as much as 5 mV. In Fig. 5, we show the
transient results obtained with PRIMA and TurboMOR models
that match four moments (q = 2). Now, both models lead to
a very accurate prediction of the original system response.
The worst case transient error is indeed below 1 mV, as
shown by Fig. 6. This example shows that matching only
two moments as in SIP [20] is not accurate enough for
some applications. TurboMOR can instead match an arbitrary
number of moments, and meet any accuracy requirement set
by the user.
C. Efficiency of the Reduced Models
We now evaluate the efficiency of the ROMs generated by
the proposed method, PRIMA, and SparseRC. In Table II, the
simulation time for the original network and the various ROMs
is reported.
Without partitioning, TurboMOR produces ROMs that are
consistently faster than PRIMA models. This is attributed to
the block diagonal structure of the reduced models, which

Fig. 8. Reduction time for SparseRC and TurboMOR with partitioning vs
number of ports. Both methods match six moments (q = 3).

reduces the cost of the LU factorizations used to perform
subsequent transient simulations. TurboMOR models are faster
by up to five times.
Comparing now the simulation times for the methods with
partitioning (SparseRC and TurboMOR with partitioning), we
observe that when two moments are matched (q = 1), the simulation times are essentially the same, which is expected since
both methods adopt the same reduction strategy. However,
when additional moments are matched, TurboMOR delivers
models that are always faster than those from SparseRC,
because of higher sparsity. SparseRC uses PRIMA to match
additional moments, which introduces some large and dense
blocks in the ROM.
D. Scalability
Finally, we investigate the scalability of TurboMOR and
existing methods with respect to network order and number
of ports. Tests are performed on the first example (on-chip
bus) for the case of six moments matched.
1) Varying Number of Ports, Constant Node-to-Port Ratio:
In the first test, we vary the number of signal lines and,
consequently, ports. Since bus length is kept constant, the
network order increases linearly with the number of ports.
The node-to-port ratio remains constant at 150.5.
Fig. 7 depicts the reduction time for TurboMOR (without
partitioning) and PRIMA versus the number of ports. We
observe that TurboMOR scales better than PRIMA, and time
savings grow as port count increases. In Fig. 8, the analysis

OYARO AND TRIVERIO - TURBOMOR

8

TABLE II
S IMULATION TIME FOR THE ROM S OBTAINED WITH THE DIFFERENT METHODS . A LL TIMES IN SECONDS .
Examples

Original
Sim.Time

1. On-chip bus
p = 256
m = 38, 528

3.00

2. ibmpg1t (RC)
p = 200
m = 25, 195

2.26

3. ibmpg2t (RC)
p = 800
m = 163, 697

29.84

4. ibmpg2t (RC)
p = 1200
m = 163, 697

30.24

5. ibmpg2t (RC)
p = 1500
m = 163, 697

30.69

6. ibmpg2t (RC)
p = 2000
m = 163, 697

30.80

q

PRIMA

SparseRC
Sim.Time

Speedup

Sim.Time

Speedup

Sim.Time

Speedup

1

0.07

42.86×

0.06

50.00×

0.07

42.86×

0.06

50.00×

2

0.76

3.95×

0.80

3.75×

0.24

12.50×

0.60

5.00×

3

3.23

0.93×

2.38

1.26×

0.54

5.56×

1.48

2.03×

1

0.23

9.83×

0.14

16.14×

0.13

17.38×

0.14

16.14×

2

0.76

2.97×

0.32

7.06×

0.33

6.85×

0.27

8.37×

3

1.92

1.18×

1.19

1.90×

0.61

3.70×

0.58

3.90×

1

4.32

6.91×

1.58

18.89×

1.60

18.65×

1.61

18.53×

2

12.98

2.30×

5.94

5.02×

8.11

3.68×

4.78

6.24×

3

27.28

1.09×

12.68

2.35×

13.63

2.19×

7.55

3.95×

1

8.58

3.52×

3.56

8.49×

3.67

8.24×

3.56

8.49×

2

28.75

1.05×

12.98

2.33×

20.30

1.49×

10.10

2.99×

3

62.51

0.48×

28.40

1.06×

34.98

0.86×

16.21

1.87×

1

13.31

2.31×

5.55

5.53×

5.98

5.13×

5.60

5.48×

2

45.32

0.68×

20.27

1.51×

35.23

0.87×

15.55

1.97×

3

104.12

0.29×

42.93

0.71×

60.40

0.51×

24.21

1.27×

1

23.08

1.33×

9.45

3.26×

10.89

2.83×

9.58

3.22×

2

81.00

0.38×

34.84

0.88×

73.17

0.42×

26.80

1.15×

3

173.67

0.18×

77.28

0.40×

121.54

0.25×

43.18

0.71×

TurboMOR
PRIMA

Reduction Time [s]

2000
1500
1000
500
0
0

50

100
150
200
Node to Port Ratio

250

300

Fig. 9. Reduction time for PRIMA and proposed method without partitioning,
as a function of the ratio of network order and number of ports.
200
TurboMOR with partitioning
sparseRC
Reduction Time [s]

TurboMOR with partitioning

Speedup

3000
2500

TurboMOR

Sim.Time

150

100

2) Varying Node-to-Port Ratio, Constant Number of Ports:
In the second test, we keep the number of ports constant to
1024, which corresponds to 512 lines. We increase the number
of nodes and, consequently, order by making the bus longer.
Fig. 9 shows the reduction time for the two methods
without partitioning (proposed and PRIMA) as a function of
the number of nodes. Beyond a certain point, the reduction
time for PRIMA increases dramatically, because the projection
matrix becomes larger than the 16 GB of memory available on
the machine. PRIMA starts resorting to slow swap memory,
and becomes very inefficient. With TurboMOR, large projection matrices are avoided. The matrices used to perform
the congruence transformations are either sparse (Cholesky
factor K) or stored in efficient factored form (Householder
reflectors in Q(i) ). This results in lower memory consumption,
and allows TurboMOR to achieve high scalability even for
very large port counts. In Fig. 10, the analysis is repeated
for TurboMOR with partitioning and SparseRC. The figure
confirms the efficiency of the proposed models, which are
faster than those generated by SparseRC especially for large
systems with many ports.

50

V. C ONCLUSION
0

0

50

100
150
200
Node to Port Ratio

250

300

Fig. 10.
As in Fig. 9, but for SparseRC and proposed method with
partitioning.

is repeated for the proposed method with partitioning and
SparseRC. Also in this case, TurboMOR scales better than
existing methods.

We introduced TurboMOR, a new model order reduction
method for large RC networks with many ports. TurboMOR
achieves moment matching via efficient Householder transformations, sparse matrix factorizations, and graph partitioning
techniques. Differently from popular methods such as PRIMA,
no large and dense projection matrices need to be computed
nor stored. This feature makes TurboMOR more efficient than
existing methods in terms of both CPU time and memory
consumption. A key novelty of the proposed method is the
sparse and block-diagonal structure of the generated models,

OYARO AND TRIVERIO - TURBOMOR

9

which makes them faster at run-time. Based on this structure,
we provide a nice interpretation of moment matching in
terms of system theory. TurboMOR models are passive by
construction, and can be cast into an equivalent RC circuit, for
seamless integration into electronic design automation tools.
Numerical results demonstrate the superior performance of
TurboMOR in reducing large passive networks with many
ports, that arise more and more frequently in practice.
A PPENDIX A
P ROOF OF MOMENT MATCHING
We prove that the reduced model (28), obtained after q
iterations of the proposed method, matches 2q moments. We
assume G invertible, since otherwise moments (4) are not
defined. If G is singular, the proposed method will still work,
but one cannot speak of moment matching.
The starting point of the proof is realization (12a)-(12b),
which is obtained from the original system (10a)-(10b) by
means of congruence transformation (11). Since (11) is invertible by construction, the transformation does not change
the transfer function nor the system moments.
The key argument of the proposed proof is the derivation
of the relation between the moments Mk of the original
system (12a)-(12b) and the moments of the inner subsystem (17) extracted by TurboMOR after one iteration. The
transfer function of the original system (12a)-(12b) can be
written as [21], [23]
H(s)

=


−1
(1)
(1)
BT1 G11 + sC11 − s2 H1 (s)
B1

(34)

(35)
(1)

is the transfer function of the inner subsystem Σ2 . The
moments of this subsystem are denoted with Nl , so we have
H1 (s) =

+∞
X

Nl sl

(36)

l=0

After substituting (4) and (36) into (34), we obtain
+∞
X

k

Mk s =

BT1

−1

+∞
X
(1)
(1)
l+2
G11 + sC11 −
Nl s
B1 (37)
l=0

k=0

For circuits, matrix B1 is typically a permutation of the identity matrix, and is thus invertible1 . We can thus rewrite (37)
as

X
+∞
+∞
X
(1)
(1)
l+2
k
G11 + sC11 −
Nl s
B−T
(38)
1 Mk s = B1
l=0

+∞ 
X
(1)
(1) −T
k
k+1
G11 B−T
1 Mk s + C11 B1 Mk s
k=0

−

+∞
X

k+l+2
Nl B−T
1 Mk s


= B1

(39)

l=0

Both sides of (39) are polynomials in s that, in order to be
equal, must have the same coefficients. Imposing the equality
between the coefficients of s0 we obtain

−1
(1)
(1)
T
G11 B−T
B1
(40)
1 M0 = B1 ⇒ M0 = B1 G11
(1)

The inverse of G11 exists since we G is non-singular. By
equating the coefficients of s1 , we have
(1)

(1)

(1)

M1 = −BT1 (G11 )−1 C11 (G11 )−1 B1

(41)

Equations (40) and (41) show that the first two moments of
(1)
the original system just depend on the matrices B1 , G11 and
(1)
C11 . Such matrices are preserved in reduced model (18),
which thus matches the first two moments of the original
system. By equating the coefficients of a generic power sr
in (39) for r ≥ 2, we obtain the recursive relation

−1
(1)
(1)
Mr = −BT1 G11
C11 B−T
1 Mr−1
r−2

−1 X
(1)
+ BT1 G11
Nl B−T
1 Mr−l−2

(42)

l=0

where

T
(1)
(1)
H1 (s) = C21
(G22 + sC22 )−1 C21

where superscript −T denotes the inverse of the transpose.
After exchanging the two series, we have

k=0

1 If B is not full rank, a correlation between some inputs exists, which can
1
be extracted before the reduction [11], making the ROM smaller and leading
to a full-rank B1 .

Equation (42) shows that the moment Mr of order r of the
original system (12a)-(12b) depends on:
(1)
(1)
1) the matrices B1 , G11 and C11 of the outer subsystem (16), which are always preserved in the reduced
model (28);
2) the moments Nl of the inner subsystem (17) up to order
r − 2.
Therefore, if one replaces the nested subsystem (17) with a
reduced model that preserves its first r − 2 moments, then the
overall model will match r moments of the original system.
By iterating this argument, it is straightforward to prove that
ROM (28) matches 2q moments of the original system.
The developed relation between the moments of the original
system and the moments of its inner subsystem (17) plays a
fundamental role in the proposed method. It allows us to match
moments recursively, two at a time, by iterative application of
the same transformation to subsystems of decreasing size. The
proposed proof is also applicable to the ROMs obtained from
other techniques such as SparseRC [23]. The main differences
between our proof and the one in [23] are two. First, the
proof in [23] considers only the first two moments, while
ours is general. Second, [23] proves moment matching for
the moments of the network admittance. Our proof is instead
based on the original impedance representation of network (1).
Our contribution therefore establishes the equivalence, from
a moment-matching perspective, of fast MOR methods (proposed, SparseRC) and PRIMA.

OYARO AND TRIVERIO - TURBOMOR

R EFERENCES
[1] H. H. Chen and J. S. Neely, “Interconnect and circuit modeling techniques for full-chip power noise analysis,” IEEE Trans. Adv. Packag.,
vol. 21, no. 3, pp. 209–215, 1998.
[2] J. M. Silva and L. M. Silveira, “Issues in model reduction of power
grids,” in Vlsi-Soc: From Systems To Silicon. Springer, 2007, pp. 127–
144.
[3] S. R. Nassif, “Power grid analysis benchmarks,” in Proceedings of the
2008 Asia and South Pacific Design Automation Conference. IEEE
Computer Society Press, 2008, pp. 376–381.
[4] B. N. Sheehan, “TICER: Realizable reduction of extracted RC circuits,”
in Proceedings of the 1999 IEEE/ACM international conference on
Computer-aided design. IEEE Press, 1999, pp. 200–203.
[5] E. J. Grimme, “Krylov projection methods for model reduction,” Ph.D.
dissertation, University of Illinois at Urbana-Champaign, 1997.
[6] M. Celik, L. Pileggi, and A. Odabasioglu, IC interconnect analysis.
Springer Science & Business Media, 2002.
[7] A. C. Antoulas, Approximation of large-scale dynamical systems.
SIAM, 2005, vol. 6.
[8] A. Odabasioglu, M. Celik, and L. T. Pileggi, “Prima: passive reducedorder interconnect macromodeling algorithm,” in Proceedings of the
1997 IEEE/ACM international conference on Computer-aided design.
IEEE Computer Society, 1997, pp. 58–65.
[9] P. Triverio, S. Grivet-Talocia, M. S. Nakhla, F. G. Canavero, and
R. Achar, “Stability, causality, and passivity in electrical interconnect
models,” IEEE Trans. Adv. Packag., vol. 30, no. 4, pp. 795–808, 2007.
[10] J. M. Silva, J. F. Villena, P. Flores, and L. M. Silveira, “Outstanding
issues in model order reduction,” in Scientific Computing in Electrical
Engineering. Springer, 2007, pp. 139–152.
[11] P. Feldmann, “Model order reduction techniques for linear systems with
large numbers of terminals,” in Proceedings of the conference on Design,
automation and test in Europe-Volume 2. IEEE Computer Society, 2004,
p. 20944.
[12] P. Liu, S. X.-D. Tan, B. Yan, and B. McGaughy, “An efficient terminal
and model order reduction algorithm,” Integration, the VLSI journal,
vol. 41, no. 2, pp. 210–218, 2008.
[13] P. Feldmann and T. Liu, “Sparse and efficient reduced order modeling
of linear subcircuits with large number of terminals,” in IEEE/ACM
International Conference on Computer Aided Design (ICCAD). IEEE,
2004, pp. 88–92.
[14] P. Benner and A. Schneider, “Model order and terminal reduction
approaches via matrix decomposition and low rank approximation,” in
Scientific Computing in Electrical Engineering SCEE 2008. Springer,
2010, pp. 523–530.
[15] P. Li and W. Shi, “Model order reduction of linear networks with massive
ports via frequency-dependent port packing,” in Proceedings of the 43rd
annual Design Automation Conference. ACM, 2006, pp. 267–272.
[16] B. Yan, S.-D. Tan, L. Zhou, J. Chen, and R. Shen, “Decentralized and
passive model order reduction of linear networks with massive ports,”
IEEE Trans. VLSI Syst., vol. 20, no. 5, pp. 865–877, 2012.
[17] P. Benner, L. Feng, and E. B. Rudnyi, “Using the superposition property
for model reduction of linear systems with a large number of inputs,”
in Proceedings of the 18th International Symposium on Mathematical
Theory of Networks & Systems, 2008.
[18] Z. Zhang, X. Hu, C.-K. Cheng, and N. Wong, “A block-diagonal
structured model reduction scheme for power grid networks,” in Design,
Automation & Test in Europe Conference & Exhibition (DATE). IEEE,
2011, pp. 1–6.
[19] B. Nouri, M. S. Nakhla, and R. Achar, “Efficient reduced-order macromodels of massively coupled interconnect structures via clustering,”
IEEE Trans. Compon., Packag., Manuf. Technol., vol. 3, no. 5, pp. 826–
840, 2013.
[20] Z. Ye, D. Vasilyev, Z. Zhu, and J. R. Phillips, “Sparse implicit projection
(SIP) for reduction of general many-terminal networks,” in Proceedings
of the 2008 IEEE/ACM International Conference on Computer-Aided
Design. IEEE Press, 2008, pp. 736–743.
[21] K. J. Kerns and A. T. Yang, “Stable and efficient reduction of large, multiport RC networks by pole analysis via congruence transformations,”
IEEE Trans. Comput.-Aided Design Integr. Circuits Syst., vol. 16, no. 7,
pp. 734–744, 1997.
[22] S. Tan and L. He, Advanced model order reduction techniques in VLSI
design. Cambridge University Press, 2007.
[23] R. Ionuţiu, J. Rommes, and W. H. Schilders, “SparseRC: sparsity
preserving model reduction for RC circuits with many terminals,” IEEE
Trans. Comput.-Aided Design Integr. Circuits Syst., vol. 30, no. 12, pp.
1828–1841, 2011.

10

[24] P. Miettinen, M. Honkala, J. Roos, and M. Valtonen, “PartMOR:
Partitioning-based realizable model-order reduction method for RLC
circuits,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst.,
vol. 30, no. 3, pp. 374–387, 2011.
[25] G. H. Golub and C. F. Van Loan, Matrix computations. JHU Press,
2012, vol. 3.
[26] F. Yang, X. Zeng, Y. Su, and D. Zhou, “RLCSYN: RLC equivalent
circuit synthesis for structure-preserved reduced-order model of interconnect,” in International Symposium on Circuits and Systems (ISCAS).
IEEE, 2007, pp. 2710–2713.
[27] C.-W. Ho, A. E. Ruehli, and P. A. Brennan, “The modified nodal
approach to network analysis,” IEEE Trans. Circuits Syst., vol. 22, no. 6,
pp. 504–509, 1975.
[28] T. Davis, “SuiteSparse,” http://faculty.cse.tamu.edu/davis/suitesparse.
html.
[29] Netlib, “LAPACK,” http://www.netlib.org/lapack/index.html.
[30] A. Zecevic and D. Siljak, “Balanced decompositions of sparse systems
for multilevel parallel processing,” IEEE Trans. Circuits Syst., vol. 41,
no. 3, pp. 220–233, 1994.
[31] Nanoscale Integration and Modeling (NIMO) Group, “Predictive technology model,” http://ptm.asu.edu/.
[32] Z. Li, P. Li, and S. R. Nassif, “IBM Power Grid Benchmarks,” http:
//dropzone.tamu.edu/∼pli/PGBench/.

