1

Simulation Algorithms with Exponential Integration
for Time-Domain Analysis of Large-Scale Power
Delivery Networks

arXiv:1505.06699v3 [cs.CE] 2 Feb 2016

Hao Zhuang, Student Member, IEEE, Wenjian Yu, Senior Member, IEEE, Shih-Hung Weng, Ilgweon Kang,
Student Member, IEEE, Jeng-Hau Lin, Student Member, IEEE, Xiang Zhang, Member, IEEE, Ryan Coutts,
and Chung-Kuan Cheng, Fellow, IEEE

Abstract—We design an algorithmic framework using matrix
exponentials for time-domain simulation of power delivery network (PDN). Our framework can reuse factorized matrices to
simulate the large-scale linear PDN system with variable stepsizes. In contrast, current conventional PDN simulation solvers
have to use fixed step-size approach in order to reuse factorized
matrices generated by the expensive matrix decomposition. Based
on the proposed exponential integration framework, we design a
PDN solver R-MATEX with the flexible time-stepping capability.
The key operation of matrix exponential and vector product
(MEVP) is computed by the rational Krylov subspace method.
To further improve the runtime, we also propose a distributed
computing framework DR-MATEX. DR-MATEX reduces Krylov
subspace generations caused by frequent breakpoints from a
large number of current sources during simulation. By virtue
of the superposition property of linear system and scaling
invariance property of Krylov subspace, DR-MATEX can divide
the whole simulation task into subtasks based on the alignments
of breakpoints among those sources. The subtasks are processed
in parallel at different computing nodes without any communication during the computation of transient simulation. The
final result is obtained by summing up the partial results among
all the computing nodes after they finish the assigned subtasks.
Therefore, our computation model belongs to the category known
as Embarrassingly Parallel model.
Experimental results show R-MATEX and DR-MATEX can
achieve up to around 14.4× and 98.0× runtime speedups over
traditional trapezoidal integration based solver with fixed timestep approach.
Index Terms—Circuit simulation, power delivery/distribution
networks, power grid, time-domain simulation, transient simulation, matrix exponential, Krylov subspace method, parallel
processing.

Manuscript received in May 2015; Revised in August, 2015 and November
2015. Accepted in Jan, 2016.
This work was supported in part by nsf-ccf 1017864, NSFC 61422402 and
NSFC 61274033.
H. Zhuang, I. Kang, J.-H. Lin, and C. K. Cheng are with the Department of
Computer Science and Engineering, University of California, San Diego, CA.
(e-mail: hao.zhuang@cs.ucsd.edu, igkang@ucsd.edu, jel252@eng.ucsd.edu,
ckcheng@ucsd.edu)
W. Yu is with Tsinghua National Laboratory for Information Science and
Technology, the Department of Computer Science and Technology, Tsinghua
University, Beijing, China. (e-mail: yu-wj@tsinghua.edu.cn)
X. Zhang and R. Coutts are with Department of Electrical and
Computer Engineering, University of California, San Diego, CA. (email:xiz110@eng.ucsd.edu, rmcoutts@eng.ucsd.edu)
S.-H. Weng was with the Department of Computer Science and Engineering, University of California, San Diego, CA. He is now with Facebook, Inc.,
Menlo Park, CA, USA. (e-mail: shweng@fb.com)

I. I NTRODUCTION

M

ODERN VLSI design verification relies heavily on
the analysis of power delivery network (PDN) to
estimate power supply noises [1]–[8]. The performance of
power delivery network highly impacts on the quality of
global, detailed and mixed-size placement [9]–[11], clock tree
synthesis [12], global and detailed routing [13], as well as timing [14] and power optimization. Lowering supply voltages,
increasing current densities as well as tight design margins
demand more accurate large-scale PDN simulation. Advanced
technologies [15], [16], three dimensional (3D) IC structures
[17]–[19], and increasing complexities of system designs all
make VLSI PDNs extremely huge and the simulation tasks
time-consuming and computationally challenging. Due to the
enormous size of modern designs and long simulation runtime
of many cycles, instead of general nonlinear circuit simulation
[20], [21], PDN is often modeled as a large-scale linear circuit
with voltage supplies and time-varying current sources [22]–
[24]. Those linear matrices are obtained by parasitic extraction
process [4], [25]–[28]. After those processes, we need timedomain large-scale linear circuit simulation to obtain the
transient behavior of PDN with above inputs.
Traditional methods in linear circuit simulation solve differential algebra equations (DAE) numerically in explicit ways,
e.g., forward Euler (FE), or implicit ways, e.g., backward
Euler (BE) and trapezoidal (TR), which are all based on low
order polynomial approximations for DAEs [29]. Due to the
stiffness of systems, which comes from a wide range of time
constants of a circuit, the explicit methods require extremely
small time step sizes to ensure the stability. In contrast, implicit
methods can handle this problem with relatively large time
steps because of their larger stability regions. However, at each
time step, these methods have to solve a linear system, which
is sparse and often ill-conditioned. Due to the requirement of
a robust solution, compared to iterative methods [30], direct
methods [31] are often favored for VLSI circuit simulation,
and thus adopted by state-of-the-art power grid (PG) solvers
in TAU PG simulation contest [32]–[34]. Those solvers only
require one matrix factorization (LU or Cholesky factorization)
at the beginning of the transient simulation. Then, at each fixed
time step, the following transient computation requires only
pairs of forward and backward substitutions, which achieves
better efficiency over adaptive stepping methods by reusing the

factorization matrix [24], [32], [34] in their implicit numerical
integration framework. However, the maximum of step size
choice is limited by the smallest distance hupper among the
breakpoints [35]. Some engineering efforts are spent to break
this limitation by sacrificing the accuracy. In our work, we
always obey the upper limit hupper of time step to maintain
the fidelity of model, which means the fixed time step h cannot
go beyond hupper in case of missing breakpoints.
Beyond traditional methods, a class of methods called
matrix exponential time integration has been embraced by
MEXP [36]. The major complexity is caused by matrix
exponential computations. MEXP utilizes standard Krylov
subspace method [37] to approximate matrix exponential and
vector product. MEXP can solve the DAEs with much higher
order polynomial approximations than traditional ones [36],
[37]. Nevertheless, when simulating stiff circuits with standard
Krylov subspace method, it requires the large dimension of
subspace in order to preserve the accuracy of MEXP approximation and poses memory bottleneck and degrade the adaptive
stepping performance of MEXP.
Nowadays, the emerging multi-core and many-core platforms bring powerful computing resources and opportunities
for parallel computing. Even more, cloud computing techniques [38] drive distributed systems scaling to thousands
of computing nodes [39]–[41], etc. Distributed computing
systems have been incorporated into products of many leading
EDA companies and in-house simulators [42]–[46]. However, building scalable and efficient distributed algorithmic
framework for transient linear circuit simulation is still a
challenge to leverage these powerful computing tools. The
papers [47], [48] show great potentials by parallelizing matrix
exponential based method to achieve the runtime performance
improvement and maintain high accuracy.
In this work, we develop a transient simulation framework
using matrix exponential integration scheme, MATEX, for
PDN simulation. Following are the challenges we need to
address. First, when the circuit is stiff, the standard Krylov
subspace has convergence problem and slows down the computation of MEVP. Second, the frequent time breakpoints due
to the transitions of PDN current sources modeling triggers
the generations of Krylov subspace. Therefore, we might gain
performance where we leverage the large time stepping, but
we also lose runtime for the small step size. Our contributions
are listed as below:
• MEVP in MATEX is efficiently computed by rational or
invert Krylov subspace method. Compared to the commonly adopted framework using TR with fixed time step
(TR-FTS), the proposed MATEX can reuse factorized
matrix at the beginning of transient simulation to perform
flexible adaptive time stepping.
• Among different Krylov subspace methods, we find rational Krylov subspace is the best strategy for MEVP in
PDN simulation. Therefore, we design R-MATEX based
on that and achieve up to around 15× runtime speedup
against the benchmarks over the traditional method TRFTS with good accuracy.
• Furthermore, DR-MATEX is designed to improve RMATEX with distributed computing resources.

– First, PDN’s current sources are partitioned into
groups based on their alignments. They are assigned
to different computing nodes. Each node runs its
corresponding PDN transient simulation task and has
no communication overhead with other nodes.
– After all nodes finish the simulation computations,
the results are summed up based on the linear
superposition property of the PDN system.
– Proposed current source partition can reduce the
chances of generating Krylov subspaces and prolong
the time periods of reusing computed subspace at
each node, which brings huge computational advantage and achieves up to 98× speedup over traditional
method TR-FTS.
The rest of this paper is organized as follows. Sec. II
introduces the background of linear circuit simulation and
matrix exponential formulations. Sec. III illustrates the Krylov
techniques to accelerate matrix exponential and vector product
computation. Sec. IV presents MATEX circuit solver and
the parallel framework DR-MATEX. Sec. V shows numerical
results and Sec. VI concludes this paper.
II. BACKGROUND
A. Transient Simulation of Linear Circuit
Transient simulation of linear circuit is the foundation
of modern PDN simulation. It is formulated as DAEs via
modified nodal analysis (MNA),
Cẋ(t) = −Gx(t) + Bu(t),

(1)

where C is the matrix for capacitive and inductive elements.
G is the matrix for conductance and resistance, and B is
the input selector matrix. x(t) is the vector of time-varying
node voltages and branch currents. u(t) is the vector of supply
voltage and current sources. In PDN, such current sources are
often characterized as pulse or piecewise-linear inputs [22],
[24] to represent the activities under the networks. To solve
Eq. (1) numerically, the system is discretized with time step h
and transformed to a linear algebraic system. Given an initial
condition x(0) from DC analysis or previous time step x(t)
and a time step h, x(t + h) can be obtained by traditional low
order approximation methods [29].
B. Traditional Low Order Time Integration Schemes
1) BE: Backward Euler based time integration scheme
(Eq.(2)) is a robust implicit first-order method.
!
C
C
+ G x(t + h) = x(t) + Bu(t + h).
(2)
h
h
2) TR: Trapezoidal based time integration scheme (Eq.(3))
is a popular implicit second-order method.
!
!
C G
C G
x(t + h) =
x(t)
(3)
+
−
h
2
h
2
u(t) + u(t + h)
.
2
It is probably the most commonly used strategy for large-scale
circuit simulation, which has higher accuracy than BE.
+ B

2

ically, turning the solution with matrix exponential operator:


−1
−2 b(t + h) − b(t)
+
x(t + h) = − A b(t + h) + A
h


b(t + h) − b(t)
. (6)
ehA x(t) + A−1 b(t) + A−2
h

3) BE-FTS and TR-FTS: Methods BE and TR with fixed
time step (FTS) h are efficient approaches, which were adopted
by the top PG solvers in 2012 TAU PG simulation contest [24],
[32]–[34]. If only one h is used for the entire simulation, the
choice is limited by the minimum breakpoint [35] distance
hupper among all the input sources. Fig. 1 (a) has 10ps as
the upper limit for h in BE-FTS and TR-FTS. When the
alignments of inputs change (as shown in Fig. 1 (b)) shift
by 5ps, the resulting upper limit for h becomes 5ps for the
approaches with fixed step size. If h is larger than the limit,
it is impossible to guarantee the accuracy since we may skip
pivot points of the inputs.

For the time step choice, breakpoints (also known as input
transition spots (TS) [47]) are the time points where slopes of
input vector change. Therefore, for Eq. (6), the maximum time
step starting from t is (ts − t), where ts is the smallest one
in T S larger than t. In matrix exponential based framework,
the limitation of time step size is not the local truncation error
(LTE), but the activities among all input sources.
III. K RYLOV S UBSPACE M ETHODS FOR M ATRIX
E XPONENTIAL AND V ECTOR P RODUCT (MEVP)
In PDN simulation, A is usually above millions and makes
the direct computation of matrix exponential eA infeasible.
The alternative way to compute the product is through Krylov
subspace method [37]. In this section, we first introduce the
background of standard Krylov subspace for MEVP. Then,
we discuss invert (I-MATEX) and rational Krylov subspace
(R-MATEX) methods, which highly improve the runtime
performance for MEVP.
A. MEXP: MEVP Computation via Standard Krylov Subspace
Method
The complexity of eA v can be reduced using Krylov
subspace method and still maintained in a high order polynomial approximation [37]. MEXP [36] uses standard Krylov
subspace, which uses A directly to generate subspace basis
through Arnoldi process (Algorithm 1). First, we reformulate
Eq. (6) into

Fig. 1. Example: Interleave two input sources to create smaller transition time.
( a) Before interleaving, the smallest transition time of the input sources is
hupper = 10ps; (b) After interleaving, the smallest transition time of the
input sources is hupper = 5ps.

C. Matrix Exponential Time Integration Scheme

(7)

b(t + h) − b(t)
h

(8)

where

The solution of Eq. (1) can be obtained analytically [29].
For a simple illustration, we convert Eq. (1) into
ẋ(t) = Ax(t) + b(t),

x(t + h) = ehA (x(t) + F(t, h)) − P(t, h),

F(t, h) = A−1 b(t) + A−2
and

(4)

when C is not singular1 ,

P(t, h) = A−1 b(t + h) + A−2

A = −C−1 G , and b(t) = C−1 Bu(t).

b(t + h) − b(t)
.
h

(9)

The standard Krylov subspace
Km (A, v) := span{v, Av, · · · , Am−1 v}

Given a solution at time t and a time step h, the solution at
t + h is
Z h
x(t + h) = ehA x(t) +
e(h−τ )A b(t + τ )dτ.
(5)

(10)

obtained by Arnoldi process has the relation
AVm = Vm Hm + hm+1,m vm+1 eT
m,

0

where Hm is the upper Hessenberg matrix


h1,1 h1,2 · · · h1,m−1 h1,m
h2,1 h2,2 · · · h2,m−1 h2,m 


h3,2 · · · h3,m−1 h3,m 
Hm = 
,
 0
··· ··· ···
···
··· 
0
0
· · · hm,m−1 hm,m

Assuming that the input u(t) is a piecewise linear (PWL)
function of t, we can integrate the last term of Eq. (5) analyt1 The assumption is to simplify the explanation in this section. After Sec.
III-B, we use I-MATEX, R-MATEX and DR-MATEX to compute the solution
of DAE without inversion of C. Therefore, the methods are suitable for
general DAE system, i.e., Eq. (1) without the assumption here.

3

(11)

(12)

Vm is a n × m matrix by (v1 , v2 , · · · , vm ), and em is the
m-th unit vector with dimension n × 1. MEVP is computed
via
ehA v ≈ βVm ehHm e1 .

reason is that the Hessenberg matrix Hm of standard Krylov
subspace tends to approximate the large magnitude eigenvalues
of A [50]. Due to the exponential decay of higher order
terms in Taylor’s expansion, such components are not the
crux of circuit system’s behavior [50], [51]. Therefore, to
simulate stiff circuit, we need to gather more vectors into
subspace basis and increase the size of Hm to fetch more
useful components, which results in both memory overhead
and computational complexity to Krylov subspace generations
for each time step. In the following subsections, we adopt
ideas from spectral transformation [50], [51] to effectively
capture small magnitude eigenvalues in A, leading to a fast
and accurate MEVP computation.

(13)

The posterior error term is
rm (h) = kβhm+1,m vm+1 eT ehHm e1 k,

(14)

where β = kvk. However, for an autonomous system Cẋ(t) =
−Gx(t) in circuit simulation, we consider the residual between Cẋ(t) and −Gx(t), which is
Cẋ(t) + Gx(t),
instead of

B. I-MATEX: MEVP Computation via Invert Krylov Subspace
Method

ẋ(t) − Ax(t)
in ẋ(t) = Ax(t). This leads to
T hHm

r(m, h) = kβhm+1,m Cvm+1 e e

e1 k

Instead of A, we use A−1 (or G−1 C) as our target matrix
to form

(15)

Km (A−1 , v) := span{v, A−1 v, · · · , A−(m−1) v}.

and helps us mitigate the overestimation of the error bound.
To generate x(t + h) by Algorithm 1, we use
[L, U] = LU Decompose(X1 ),

(17)

Intuitively, by inverting A, the small magnitude eigenvalues
become the large ones of A−1 . The resulting Hm is likely to
capture these eigenvalues first. Based on Arnoldi algorithm,
the invert Krylov subspace has the relation

(16)

where
X1 = C and X2 = G

A−1 Vm = Vm Hm + hm+1,m vm+1 eT
m.

as inputs for standard Krylov subspace. The error budget ǫ and
Eq. (14) are used to determine the convergence when time step
is h and Krylov subspace dimension is j (from line 11 to line
14 in Algorithm 1).

(18)

The matrix exponential eA v is calculated as
eA v ≈ βVm ehHm e1 .
−1

(19)

To put this method into Algorithm 1 is just by modifying the
inputs X1 = G for the LU decomposition in Eq. (16), and
X2 = C. In the line 16 of Algorithm 1,

Algorithm 1: MATEX Arnoldi
Input: L, U, X2 , h, t, x(t), ǫ, P(t, h), F(t, h)
Output: x(t + h), Vm , H, v
1 v = x(t) + F(t, h);
v
2 v1 = kvk ;
3 for j = 1 : m do
4
w = U\(L\(X2 vj )) ; /* a pair of forward
and backward substitutions.
*/
5
for i = 1 : j do
6
hi,j = wT vi ;
7
w = w − hi,j vi ;
8
end
9
hj+1,j = kwk;
w
10
vj+1 = hj+1,j
;
11
if r(j, h) < ǫ then
12
m = j;
13
break;
14
end
15 end
hH
16 x(t + h) = kvkVm e
e1 − P(t, h);

H = H−1
m
for the invert Krylov version. The posterior error approximation [47] is
−1 hH
rm (h) = kβhm+1,m Avm+1 eT
e1 k,
m Hm e

(20)

which is derived from residual based error approximation
in [51]. However, as mentioned in Sec. III-A, we consider
the residual of (Cẋ(t) + Gx(t)), instead of (ẋ(t) − Ax(t)),
which leads to
−1 hHm
e1 k.
r(m, h) = kβhm+1,m Gvm+1 eT
m Hm e

(21)

We use Eq. (21) for the line 11 of Alg. 1.
C. R-MATEX: MEVP Computation via Rational Krylov Subspace Method
The shift-and-invert Krylov subspace basis [50] is designed
to confine the spectrum of A. Then, we generate Krylov
subspace via
Km ((I − γA)−1 , v) :=
(22)
−1
−(m−1)
span{v, (I − γA) v, · · · , (I − γA)
v},

The standard Krylov subspace may not be efficient when
simulating stiff circuits [36], [49]. For the accuracy of approximation of eA v, a large dimension of Krylov subspace
basis is required, which not only brings the computational
complexity but also consumes huge amount of memory. The

where γ is a predefined parameter. With this shift, all the
eigenvalues’ magnitudes are larger than one. Then the inverse
4

Proof. 2 If λ and u are an eigenvalue and eigenvector of a
generalized eigenvalue problem

limits the magnitudes smaller than one. According to [50],
[51], the shift-and-invert basis for matrix exponential based
transient simulation is not very sensitive to γ, once it is
set to around the order near time steps used in transient
simulation. The similar idea has been applied to simple power
grid simulation with matrix exponential method [52]. Here,
we generalize this technique and integrate into MATEX. The
Arnoldi process constructs Vm and Hm . We have
(I − γA)−1 Vm = Vm Hm + hm+1,m vm+1 eT
m.

−Gu = λCu.
Then, x = etλ u is the solution of Cẋ = −Gx.
Because we do not need to compute C−1 explicitly during
Krylov subspace generation, I-MATEX and R-MATEX are
regularization-free. Instead, we factorize G for invert Krylov
subspace basis generation (I-MATEX), or (C + γG) for
rational Krylov subspace basis (R-MATEX).3 Besides, their
Hessenberg matrices Eq. (12) are invertible, which contain
corresponding important generalized eigenvalues/eigenvectors
from matrix pencil (−G, C), and define the behavior of linear
dynamic system in Eq. (1) of interest.

(23)

We can project the eA onto the rational Krylov subspace as
follows.
eAh v ≈ βVm eh

I−H−1
m
γ

e1 .

(24)

In the line 16 of Algorithm 1,
H=

E. Comparisons among Different Krylov Subspace Algorithms
for MEVP Computation

I − H−1
m
.
γ

In order to observe the error distribution versus dimensions
of standard, invert, and rational Krylov subspace methods for
MEVP, we construct a RC circuit with stiffness

Following the same procedure [47], [51], the posterior error
approximation is derived as
rm (h) = kβhm+1,m

I − γAm
−1 hH
vm+1 eT
e1 k.
m Hm e
γ

Re(λmin )
= 4.7 × 106 ,
Re(λmax )

(25)

Note that in practice, instead of computing (I−γA)−1 directly,
(C + γG)−1 C is utilized. The corresponding Arnoldi process
shares the same skeleton of Algorithm 1 with input matrices

where λmax = −8.49 × 1010 and λmin = −3.98 × 1017 are
the maximum and minimum eigenvalues of A = −C−1 G.
Fig. 2 shows the relative error reductions along the increasing
Krylov subspace dimension. The error reduction rate of rational Krylov subspace is the best, while the one of standard
Krylov subspace requires huge dimension to capture the same
level of error. For example, it costs almost 10× of the size
to achieve around relative error 1% compared to Invert and
Rational Krylov subspace methods. The relative error is

X1 = (C + γG)
for the LU decomposition Eq. (16), and
X2 = C.
The residual estimation is

||ehA v − βVm ehHm e1 ||
,
||ehA v||

C + γG
−1 hHm
e1 k. (26)
vm+1 eT
r(m, h) = kβhm+1,m
m Hm e
γ

where h = 0.4ps, γ = 10−13 . The matrix A is a relatively
small matrix and computed by MATLAB expm function. The
result of ehA v serves as the baseline for accuracy. The relative
error is the real relative difference compared to the analytical
solution ehA v of the ODE
dx
= Ax
dt
with an initial vector v, which is generated by MATLAB rand
function.
The error reduction rate of standard Krylov subspace is the
worst, while the rational Krylov subspace is the best. It is the
reason that we prefer rational Krylov subspace (R-MATEX).
The relative errors of BE, TR and FE are 0.0594, 0.4628, and
2.0701 × 104 , respectively. The large error of FE is due to
the instability issue of its low order explicit time integration
scheme. In Fig. 2, when m = 3, standard, invert and rational
Krylov subspace methods have 0.8465, 0.0175, and 0.0065,
respectively. It illustrates the power of matrix exponential

Then, we plug Eq. (26) into the line 11 of Algorithm 1.
D. Regularization-Free MEVP Computation
When C is a singular matrix, MEXP [36] needs the regularization process [53] to remove the singularity of DAE in Eq.
(1). It is because MEXP needs factorize C directly to form
the input X1 for Algorithm 1. This brings extra computational
overhead when the case is large [53]. It is not necessary if
we can obtain the generalized eigenvalues and corresponding
eigenvectors for matrix pencil (−G, C). Based on [54], we
derive the following lemma,
Lemma 1. Considering a homogeneous system
Cẋ = −Gx.
u and λ are the eigenvector and eigenvalue of matrix pencil
(−G, C), then
x = etλ u

2 We repeat the proof from [54] with some modifications for our formulation.
3 It is also applied to the later work of DR-MATEX in Sec. IV-B.

is a solution of the system.
5

10 0
Standard Krylov
Invert Krylov
Rational Krylov

10 0

10 -2

10 -5

Relative Error

Relative Error

10 -1

10 -3

10 -4

10 -10

10 -15

10 -5

10

20

30

40

50

0

10 -20
10 -20

60

m

50
10 -15

10 -10

h

Fig. 2. The relative error vs. dimensional m of different Krylov subspace
||ehA v−βVm ehHm e1 ||
, where h = 0.4ps,
methods. The relative error is
||ehA v||
γ = 10−13 . Note: The relative error is the difference compared to analytical
= Ax with an initial vector v, which is
solution ehA v of the ODE dx
dt
generated by MATLAB rand function, and its entries are positive numbers
in (0, 1].

10 0

4.

The

error

of MEVP via standard Krylov Subspace:
vs. time step h and dimension of standard
Krylov subspace basis (m). The standard Krylov subspace approximates
the solution well in extremely small h, since it captures the important
eigenvalues and eigenvectors of A at that region. However, the small h is
not useful for the circuit simulation. For large h, it costs large m to reduce
the error.
||ehA v−βVm ehHm ||
||ehA v||

Standard Krylov (4.7e6)
Invert Krylov (4.7e6)
Rational Krylov (4.7e6)
Standard Krylov (4.7e10)
Invert Krylov (4.7e10)
Rational Krylov (4.7e10)

10 -1

Relative Error

Fig.

m

100

10 -2

10 0
10

-3

Relative Error

10 -4
10 -5

10

20

30

40

50

10 -5

10 -10

60

m

0

10 -15
10 -20

Fig. 3. The relative error vs. dimension m of different Krylov subspace
hA
hHm
e1 ||
me
methods. The relative error is ||e v−βV
, where h = 0.4ps,
||ehA v||
−13
γ = 10
. The rational Krylov subspace has very stable error reduction
rate. The number in the bracket represents the stiffness value of the system.

20

m

10 -15

10 -10

h

Fig.

5.

The

error

of

MEVP

via

40

invert

Krylov

Subspace:

||ehA v−βVm ehHm ||
||ehA v||
−1

vs. time step h and dimension of invert Krylov
subspace basis (m). Compared to Fig. 4, invert Krylov subspace method
reduces the errors for large h.

method. Our proposed methods are all stable and can achieve
improved error numbers.
In order to observe the different stiffness effects on Krylov
subspace methods, we change the entries in C and G to make
the different stiffness value 4.7 × 1010 . Fig. 3 illustrates the
stable reduction rate of rational method. The stiffness degrades
the performance of standard Krylov subspace method. Both invert and rational Krylov subspace methods are good candidates
for stiff circuit system.
Regarding the relative error distributions vs. time step h and
dimension m, Fig. 4, Fig. 5, and Fig. 6 are computed by standard, invert, and rational Krylov subspaces (γ = 5 × 10−13 ),
respectively. Fig. 4 shows that the errors generated by standard
Krylov subspace method has flat region with high error values
in time-step range of interests. The small (‘unrealistic’) time
step range has small error values. Compared to Fig. 4, invert
(Fig. 5) and rational (Fig. 6) Krylov subspace methods reduce
errors quickly for large h. The explanation is that a relatively

small portion of the eigenvalues and corresponding invariant
subspaces determines the final result (vector) when time step
h is larger [50], which are efficiently captured by invert and
rational Krylov subspace methods.
The error of rational Krylov subspace is relatively insensitive to γ when it is selected between the time-step range of
interests (Fig. 7). Above all, rational Krylov (R-MATEX) and
invert Krylov (I-MATEX) subspace methods have much better
performance than standard version. When we deal with stiff
cases, standard Krylov subspace is not a feasible choice due
to the large dimension m of Krylov subspace, which causes
huge memory consumption and poor runtime performance.
6

Algorithm 2: MATEX Circuit Solver
Input: C, G, B, u, ǫ, and time span T .
Output: The set of x from [0, T ].
1 Set X1 , X2 ;
2 t = 0;
3 x(t) =DC analysis;
4 [L, U] = LU Decompose(X1 );
5 while t < T do
6
Compute maximum allowed step size h;
7
Update P(t, h), F(t, h);
8
Obtain x(t + h) by Algorithm 1 with inputs
[L, U, X2 , h, t, x(t), ǫ, P(t, h), F(t, h)];
9
t = t + h;
10 end

Relative Error

10 0

10 -5

10 -10

0

10 -15
10 -20

20

m

10 -15

10 -10

h

Fig. 6.
hA

||e

The
h

error

of MEVP via

40

rational

Krylov

Subspace:

I−H−1
m
γ
e1 ||

v−βVm e
||ehA v||

, where γ = 5 × 10−13 , vs. time step h and
dimension of rational Krylov subspace basis (m). Compared to Fig. 4, rational
Krylov subspace method reduces the errors for large h as Fig. 5.

B. DR-MATEX (Distributed R-MATEX Framework) by Decomposition of Input Sources, Linear Superposition, and Parallel Computation Model
1) Motivation:
There are usually many input sources in PDNs as well as
their transition activities, which might narrow the regions for
the stepping of matrix exponential based method due to the
unaligned breakpoints. In other words, the region before the
next transition ts may be shortened when there are a lot of
activities from the input sources. It leads to more chances of
generating new Krylov subspace bases. We want to reduce
the number of subspace generations and improve the runtime
performance.4

Relative Error

0.02
0.015
0.01
4
9

0.005

2) Treatment and Methodology:
In matrix exponential based integration framework, we can
choose any time spot t + h ∈ [t, ts ] with computed Krylov
subspace basis. The solution of x(t+h) is computed by scaling
the existing Hessenberg matrix H with the time step h as
below

14
19 m
24

0
10 5

Fig.
hA

||e

7.

10 0

The
h

error

10 -5

of

γ

MEVP

10 -10

via

10 -15

rational

Krylov

Subspace

I−H−1
m
γ
e1 ||

x(t + h) = kvkVm ehH e1 − P(t, h).

v−βVm e
||ehA v||

, where h = 4ps. The flat region shows
the error is actually relatively insensitive to γ, when γ is in the range of step
size h of interests.

This is an important feature for computing the solutions
at intermediate time points without generating the Krylov
subspace basis, when there is no current transition. Besides,
since the PDN is linear dynamical system, we can utilize
the well-known superposition property of linear system and
distributed computing model to tackle this challenge.
To illustrate our distributed version of MATEX framework,
we first define three terms to categorize the breakpoints of
input sources:

IV. MATEX F RAMEWORK
A. MATEX Circuit Solver
We incorporate matrix exponential based integration scheme
with Krylov subspace method into our MATEX framework,
which is summarized in Algorithm 2. We set X1 and X2 in
Line 1 based on the choice of Krylov subspace method as
follows,
•
•

(27)

•
•

I-MATEX: X1 = G, X2 = C
R-MATEX: X1 = C + γG, X2 = C

•

Local Transition Spot (LT S): the set of T S at an input
source to the PDN.
Global Transition Spot (GT S): the union of LT S among
all the input sources to the PDN.
Snapshot: a set GT S \ LT S at one input source.

If we simulate the PDN with respect to all the input sources,
the points in the set of GT S are the places where generations

For linear system of PDN, the matrix factorization in line
4 is only performed once, and the matrices L and U are
reused in the while loop from line 5 to line 10. Line 8 uses
Arnoldi process with corresponding inputs to construct Krylov
subspace as shown in Algorithm 1.

4 The breakpoints also put the same constraint on TR-FTS and BE-FTS.
However, their time steps are fixed already, which refrains them from reaching
this problem in the first place.

7

Global Transition Spots (GTS)

0

#%

#$

Local Transition Spots
(LTS) at #1.1 in Group 1

Local Transition Spots
(LTS) at #1.2 in Group 4

!""

Local Transition Spots
(LTS) at #2.1 in Group 2

#&

Local Transition Spots
(LTS) at #2.2 in Group 3

Fig. 8. Part of a PDN model with input sources from Fig. 9.
Local Transition Spots
(LTS) at #3 in Group 4

of Krylov subspace cannot be avoided. For example, there are
three input sources in a PDN (Fig. 8). The input waveforms are
shown in Fig. 9. The first line is GT S, which is contributed by
the union of LT S in input sources #1, #2 and #3. However,
we can partition the task into subtasks by simulating each
input source individually. Then, each subtask generates Krylov
subspaces based on its own LT S and keeps track of Snapshot
for the later usage of summation via linear superposition.
Between two LTS points t and t + h, the Snapshot points

t_width
!"delay

Fig. 10. Grouping of “Bump” shape transitions for sub-task simulation.
The matrix exponential based method can utilize adaptive stepping in each
LT S and reuse the Krylov subspace basis generated at the latest point in
LT S. However, traditional methods (TR, BE, etc.) still need to do time
marching, either by pairs of forward and backward substitutions and proceed
with fixed time step, or by re-factorizing matrix and solving linear system for
adaptive stepping. (Pulse input information: tdelay : delay time; trise : rise
time; twidth : pulse width; tf all : fall time; and tperiod : period).

can reuse the Krylov subspace generated at t. For each node,
the chances of generation of Krylov subspaces are reduced.
The time periods of reusing latest Krylov subspaces are
enlarged locally and bring the runtime improvement. Besides,
when subtasks are assigned, there is no communication among
the computing nodes, which leads to so-called Embarrassingly
Parallel computation model.

Global Transition Spots (GTS)

into one set. For example, the input source #1 of Fig. 9 is
divided to #1.1 and #1.2 in Fig. 10. The input source #2
in Fig. 9 is divided to #2.1 and #2.2 in Fig. 10. Therefore,
there are four groups in Fig. 10, Group 1 contains LT S#1.1.
Group 2 contains LT S#2.1. Group 3 contains LT S#2.2.
Group 4 contains LT S#1.2 and #3. Our proposed framework
MATEX is shown in Fig. 11. After pre-computing GT S and
decomposing LT S based on “bump” shape (Fig. 10), we group
them and form LT S #1 ∼ #K.6

part of snapshots at LTS #1

4) MATEX Scheduler in DR-MATEX:
In DR-MATEX, the role of MATEX scheduler is just to
send out GT S and LT S to different MATEX slave nodes
and collect final results after all the subtasks of transient
simulation are finished. The node number is based on the
total number of subtasks, which is the group number after
PDN source decomposition. Then the simulation computations
are performed in parallel. Each node has its own inputs. For
example, Node#k has GT S, LT S#k, Pk and Fk , which
contain the corresponding b for node k. Scheduler does not
need to do anything during the transient simulation, since there
are no communications among nodes before the stage of “write
back” (in Fig. 11), by when all nodes complete their transient
simulations.
Within each slave node, the circuit solver (Algorithm 3)
computes transient response with varied time steps. Solutions
are obtained without re-factorizing matrix during the computation of transient simulation. The computing nodes write back
the results and inform the MATEX scheduler after finishing
their own transient simulation.

Input Source #1
Local Transition Spots
(LTS) at #1
Input Source #2
Local Transition Spots
(LTS) at #2
Input Source #3
Local Transition Spots
(LTS) at #3

part of snapshots at LTS #3

Fig. 9. Illustration of input transitions. GT S: Global Transition Spots; LT S:
Local Transition Spots; Snapshots: the crossing positions by dash lines and
LTS #k without solid points.

3) More Aggressive Tasks Decomposition: We divide the
simulation task based on the alignments of input sources.
More aggressively, we can decompose the task according to
the “bump” shapes of the input sources.5 We group the input
sources, which have the same
(tdelay , trise , tf all , twidth )
5 IBM

t_fall
!"period

t + h1 < t + h2 < · · · < t + hl ∈ (t, t + h]

0

t_rise

6 There are alternative decomposition strategies. It is also easy to extend
the work to deal with different input waveforms. We try to keep this part as
simple as possible to emphasize our framework.

power grid benchmarks provide the pulse input model in SPICE

format.

8

MATEX Scheduler:
Circuit
Global Transitions Spots
(GTS)

MATEX Slave Node #1:
Circuit (Local Copy)
Global Transitions Spots
(Local Copy of GTS)
Local Transition Spots
(LTS) #1

…

Local Transition Spots
(LTS) #1
Local Transition Spots
(LTS) #2

Circuit Solver

MATEX Slave Node #K:
Circuit (Local Copy)

...

Local Results

write back

Local Transition Spots
(LTS) #K

Global Transitions Spots
(Local Copy of GTS)

Local Results Collection

Local Transition Spots
(LTS) #K

MATEX Slave Node #2:
Circuit (Local Copy)
Global Transitions Spots
(Local Copy of GTS)
Local Transition Spots
(LTS) #2
Circuit Solver

Circuit Solver
Superposition based on
GTS

write
back

Local Results

write
back

Local Results

Global Results

…
Fig. 11. DR-MATEX: The distributed MATEX framework using R-MATEX circuit solver.

C. Runtime Analysis of MATEX PDN Solver
Suppose we have the dimension of Krylov subspace basis
m on average for each time step and one pair of forward and
backward substitutions consumes runtime Tbs . The total time
of serial parts is Tserial , which includes matrix factorizations,
result collection, etc. For x(t + h), the evaluation of matrix
exponential with ehHm is TH , which is in proportion to the
time complexity O(m3 ). Besides, we need extra Te to form
x(t + h), which is proportional to O(nm2 ) by βVm ehHm e1 .
Given K points of GT S, without decomposition of input
sources, the runtime is

Algorithm 3: DR-MATEX: The distributed MATEX
framework using R-MATEX at Node#k.
Input: LT S#k, GT S, Pk , Fk , error tolerance Etol , and
simulation time span T .
Output: Local solution x along GT S in node
k ∈ [1, · · · , S], where S is the number of nodes
1 t = 0, X1 = C + γG, and X2 = C;
2 x(t) = Local Initial Solution;
3 [L, U] = LU Decompose(X1 );
4 while t < T do
5
Compute maximum allowed step size h based on
GT S;
6
if t ∈ LTS#k then
/* Generate Krylov subspace for the

KmTbs + K(TH + Te ) + Tserial .

After dividing the input transitions and sending to enough
computing nodes, we have k points of LT S for each node
based on feature extraction and grouping (e.g., k = 4 for one
“bump” shape feature). The total computation runtime is

point at LT S#k and compute x(t + h)
*/
7

8
9
10

12
13
14
15

kmTbs + K(TH + Te ) + Tserial ,

[x(t + h), Vm , Hm , v] =
MATEX Arnoldi(L, U, X2 ,
h, t, x(t), ǫ, Pk (t, h), Fk (t, h));
alts = t;
end
else
/* Obtain x(t + h) at Snapshot with
computed Krylov subspace

11

(28)

(29)

where K(TH + Te ) contains the portion of computing
Snapshot in DR-MATEX mode. The speedup of DR-MATEX
over single MATEX is
Speedup =
*/

KmTbs + K(TH + Te ) + Tserial
.
kmTbs + K(TH + Te ) + Tserial

(30)

For R-MATEX, we have small m. Besides, Tbs is relatively
larger than (TH + Te ) in our targeted problem. Therefore, the
most dominating part is the KmTbs in Eq. (28). We can always
decompose input source transitions, and make k smaller than
K.
In contrast, suppose the traditional method with fixed step
size has N steps for the entire simulation, the runtime is

ha = t + h − alts ;
x(t + h) = kvkVm eha Hm e1 − Pk (t, h);
end
t = t + h;
end

N Tbs + Tserial .
9

time stepping. Actually, multiple factorized matrices can be
deployed [55], [56]. We can choose one of them during the
stepping. The problem is the memory and runtime overhead
for the multiple matrix factorizations. Another point is if large
time step h′ is chosen, the standard low order scheme cannot
maintain the accuracy.
Experiment is conducted on a single computing node. In
Table II, we record the total simulation runtime Total(s),
which includes the processes of DC and transient simulation,
but excludes the non-numerical computation before DC, e.g.,
netlist parsing and matrix stamping. We also record the part
of transient simulation Tran(s), excluding DC analysis and
LU decompositions. The speedup of I-MATEX is not as
large as R-MATEX, because I-MATEX with a large spectrum of A generates large dimension m of Krylov subspace.
Meanwhile, the step size is not large enough to let it fully
harvest the gain from time marching with stepping. In contrast,
R-MATEX needs small dimension numbers m of rational
Krylov subspace, which ranges from 2 to 8 in those cases.
Therefore, they can benefit from large time stepping, shown as
SPDPrtr . For ibmpg4t, R-MATEX achieves maximum speedup
resulted from the relatively small number of breakpoints in that
benchmark, which is around 44 points, while the majority of
others have over 140 points.
In Table II, our single mode R-MATEX achieves the average
speedup 5× over TR-FTS. Note the average speedup number
of single mode R-MATEX over TR-FTS for the original IBM
benchmark (ibmpg1t∼ibmpg6t) is less than the speedup of
the new test cases (ibmpg1t new∼ibmpg6t new). As we mentioned before, ibmpg1t new∼ibmpg6t new have harsher input
constraints, making the available step size only 1ps. Therefore,
the adaptive stepping by R-MATEX is more beneficial to
the runtime performance in ibmpg1t new∼ibmpg6t new than
ibmpg1t∼ibmpg6t.

TABLE I
S PECIFICATIONS OF IBM POWER GRID BENCHMARKS .
Design
ibmpg1t
ibmpg2t
ibmpg3t
ibmpg4t
ibmpg5t
ibmpg6t

#R
41K
245K
1.6M
1.8M
1.6M
2.4M

#C
11K
37K
201K
266K
473K
761K

#L
277
330
955
962
277
381

#I
11K
37K
201K
266K
473K
761K

#V
14K
330
955
962
539K
836K

#Nodes
54K
165K
1.0M
1.2M
2.1M
3.2M

Then, the speedup of distributed DR-MATEX over the traditional method is
N Tbs + Tserial
.
(31)
Speedup′ =
kmTbs + K(TH + Te ) + Tserial
Note that, when the minimum distance among input source
breakpoints decreases, large time span or many cycles is
required to simulate PDNs, the schemes with such uniform
step size would degrade runtime performance furthermore due
to the increase of N . In contrast, in MATEX PDN solver,
K is not so sensitive to such constraints. Besides, k can be
maintained in a small number based on the decomposition
strategy. Therefore, the speedups of our proposed methods
tend to be larger when the simulation requirements become
harsher.
V. E XPERIMENTAL R ESULTS
We implement all the algorithms in MATLAB R2014b7
and use UMFPACK package for LU factorization. First, we
compare I-MATEX, R-MATEX and TR in order to show our
runtime improvements in single machine framework in Table
II. Second, we show our distributed framework DR-MATEX
achieves large speedups in Table III. The experiments are
conducted on the server with Intel(R) Xeon (R) E5-2640 v3
2.60GHz processor and 125GB memory.

B. Performance of DR-MATEX in Sec. IV-B

A. Performance of I-MATEX and R-MATEX in Sec. IV-A

We test our distributed DR-MATEX in the following experiments with the same IBM power grid benchmarks. These
cases have many input transitions (GT S) that limit step sizes
of R-MATEX. We divide the region before the computation of
simulation. We decompose the input sources by the approach
discussed in Sec. IV-B3 and obtain much fewer transitions
of LT S for computing nodes. The original input source
numbers are over ten thousand in the benchmarks. However,
based on “bump” feature (as shown in Fig. 10), we obtain
a fairly small numbers for each computing node, which is
shown as Group # in Table III. (Now, the fact that hundred
machines to process in parallel is quite normal [38], [57] in
the industry.) We pre-compute GT S and LT S groups and
assign sub-tasks to corresponding nodes8 . MATEX scheduler
is only responsible for simple superposition calculation at the
end of simulation. Since the slave nodes are in charge of all
the computing procedures (Fig. 11) for the computation of

We compare our proposed I-MATEX and R-MATEX against
the popular TR-FTS on the IBM power grid benchmarks [22].
Among the current sources, the smallest interval between two
breakpoints is hupper = 10ps, which puts the upper limit
of the TR’s step size. All of these cases have very large
numbers of input current sources. Table I shows the details of
each benchmark circuit of which size ranges from 54K up to
3.2M. The simulation time is 10ns. From ibmpg1t to ibmpg6t,
TR uses fixed step size in 10ps. We also change the IBM
power grid benchmark to make the smallest distance among
breakpoints 1ps by interleaving input sources’ breakpoints
(similar as Fig. 1). Therefore, the fixed step size method
can only use at most 1ps. The names of those benchmarks
are ibmpg1t new, ibmpg2t new, ibmpg3t new, ibmpg4t new,
ibmpg5t new and ibmpg6t new.
After DC analysis in TR-FTS, we LU factorize matrix
once for the later transient simulation, which only contains

8 Based on the feature of input sources available, the preprocessing is very
efficient, which takes linear time complexity to obtain GTS, LTS and separates
the sources into different groups.

7 Measurements reported are on MATLAB implementations. They are
subject to limitations and are not directly comparable to C++ implementations
reported in literature such as [44].

10

TABLE II
P ERFORMANCE COMPARISONS ( SINGLE COMPUTING NODE ): TR-FTS, I-MATEX, AND R-MATEX. DC(s): RUNTIME OF DC ANALYSIS ( SECONDS ); mI :
T HE MAXIMUM m OF K RYLOV SUBSPACE IN I-MATEX. Tran(s): RUNTIME OF TRANSIENT SIMULATION AFTER DC ( SECONDS ), EXCLUDING THE
MATRIX FACTORIZATION RUNTIME ; Total(s): R UNTIME OF OVERALL TRANSIENT SIMULATION ( SECONDS ); Df(uV): M AXIMUM AND AVERAGE VOLTAGE
DIFFERENCES COMPARED TO PROVIDED SOLUTIONS ( U V); mR : T HE MAXIMUM m OF K RYLOV SUBSPACE IN R-MATEX SPDPr
tr : S PEEDUP OF
R-MATEX OVER TR-FTS WITH RESPECT TO Tran(s); SPDPri : S PEEDUP OF R-MATEX OVER I-MATEX WITH RESPECT TO Tran(s).
Design
ibmpg1t
ibmpg2t
ibmpg3t
ibmpg4t
ibmpg5t
ibmpg6t
ibmpg1t new
ibmpg2t new
ibmpg3t new
ibmpg4t new
ibmpg5t new
ibmpg6t new
Average

DC(s)
0.2
0.8
16.4
13.5
9.0
15.3
0.2
0.9
16.3
18.3
10.5
13.1
—

TR-FTS
Tran(s)
Total(s)
5.7
6.00
40.0
41.9
263.2
295.0
460.8
501.9
476.6
498.0
716.0
749.1
51.3
51.7
431.4
433.5
3716.5
3749.0
5044.6
5085.3
5065.9
5110.1
7015.3
7059.7
—
—

mI
30
28
29
29
30
25
30
28
29
29
30
25
—

Tran(s)
28.8
130.0
1102.5
433.8
1934.4
2698.9
27.2
114.9
1219.3
753.5
2494.0
3647.9
—

I-MATEX
Total(s)
28.9
130.9
1115.1
458.2
1944.5
2713.7
27.4
115.7
1232.6
776.4
2504.7
3663.1
—

Df(uV)
58\9.8
92\10.5
95\20.4
101\39.3
29\5.6
39\8.6
58\9.8
49\10.5
95\20.4
101\39.3
30\5.6
39\8.6
65\15.7

their own transient simulation tasks, and have no communications with others, our framework falls into the category
of Embarrassingly Parallelism model. We can easily emulate
the multiple-node environment. We simulate each group using
the command “matlab -singleCompThread” in our server. We
record the runtime numbers for each process (slave nodes) and
report the maximum runtime as the total runtime “Total(s)”
of DR-MATEX in Table III. We also record “pure transient
simulation” as “Tran(s)”, which is the maximum runtime of
the counterparts among all computing nodes.
For TR-FTS, we use h = 10ps, so there are 1,000 pairs
of forward and backward substitutions during the process
of pure transient simulation for ibmpg1t∼ibmpg6t; We use
h = 1ps for ibmpg1t new∼ibmpg6t new. Therefore, we
have 10,000 pairs of forward and backward substitutions for
stepping. In DR-MATEX, the circuit solver uses R-MATEX
with γ = 10−10 , which is set to sit among the order of varied
time steps during the simulation (since Sec. III-E discusses the
insensitivity of γ around the step size of interests). TR-FTS is
not distributed because it has no gain by dividing the current
source as we do for the DR-MATEX. TR-FTS cannot avoid the
repeated pairs of forward and backward substitutions. Besides,
adaptive stepping for TR-FTS only degrades the performance,
since the process requires extra matrix factorizations.
In Table III, our distributed mode gains up to 98× for
the pure transient computing. The average peak dimension
m of rational Krylov subspace is 7. The memory overhead
ratio for each node (around 1.6× over TR-FTS on average) is
slightly larger, which is worthwhile with respect to the large
runtime improvement. With the huge reduction of runtime for
Krylov subspace generations, the serial parts, including LU
and DC, play more dominant roles in DR-MATEX, which
can be further improved using advance matrix solvers, such
as [58].
VI. C ONCLUSIONS

AND

mR
5
5
5
5
5
5
5
5
5
6
5
6
—

R-MATEX
Tran(s)
Total(s)
10.1
10.3
35.6
37.4
275.5
301.0
200.5
239.1
383.1
401.9
773.5
800.5
11.7
12.1
43.3
44.9
481.7
508.2
350.9
387.2
746.2
766.4
895.1
923.1
—
—

Df(uV)
45\6.8
45\6.8
95\18.5
99\34.2
29\4.4
33\5.6
53\6.9
33\5.6
95\18.9
99\34.2
30\4.4
33\7.3
57\12.8

Speedups
SPDPrtr
SPDPri
0.6×
2.9×
1.1×
3.7×
1.0×
4.0×
2.3×
2.2×
1.2×
5.0×
0.9×
3.5×
4.4×
2.3×
10.0×
2.7×
7.7×
2.5×
14.4×
2.1×
6.8×
3.3×
7.8×
4.1×
5×
3×

tial integration scheme. We visualize the error distributions to
show the advantages of using rational (R-MATEX) and invert
(I-MATEX) Krylov subspace methods for matrix exponential
and vector product (MEVP) over standard Krylov subspace
method (MEXP). For the PDN simulation, our time integration
scheme can perform adaptive time stepping without repeating
matrix factorizations, which cannot be achieved by traditional
methods using implicit numerical integration with fixed timestep scheme. Compared to the commonly adopted framework
TR with fixed time step (TR-FTS), our single mode framework
(R-MATEX) gains runtime speedup up to around 15×. We also
show that the distributed MATEX framework (DR-MATEX)
leverages the superposition property of linear system and
decomposes the task based on the feature of input sources,
so that we reduce chances of Krylov subspace generations
for each node. We achieve runtime improvement up to 98×
speedup.
We show the exponential integration with Krylov subspace
methods maintains high order accuracy and flexible time
stepping ability. The exponential integration framework was
actually mentioned by the very early work in circuit simulation
algorithms [29], but it had not attracted too much attention due
to the high computational complexities of matrix exponential
during that time. Nowadays, the progress of Krylov subspace
methods provides efficient way to compute matrix exponential
and vector product, so that we can utilize certain features of
exponential integration, which are hardly obtained by traditional time integration schemes. Exponential integration can
also serve as stable explicit schemes [59], [60] for general
dynamical systems. It is a promising framework for the future
circuit simulation algorithms and software. The opportunities
of parallel and distributed computing with the cutting-edge
multi-core and many-core hardware are also worth exploring
for the further parallelism and runtime improvement.
ACKNOWLEDGMENT

F UTURE D IRECTIONS

In this work, we propose an efficient framework MATEX for
accurate PDN time-domain simulation based on the exponen-

We thank Prof. Scott B. Baden, Prof. Mike Botchev, Dr.
Quan Chen, Dr. Peng Du, Dr. Martin Gander, Prof. Nicholas
11

TABLE III
T HE PERFORMANCE OF DR-MATEX (D ISTRIBUTED R-MATEX). Group #: G ROUP NUMBER OF THE TESTCASES . T HIS NUMBER REPRESENTS THE
TOTAL NUMBER OF SIMULATION SUB - TASKS FOR THE DESIGN ; Tran(s): R UNTIME OF TRANSIENT SIMULATION AFTER DC ( SECONDS ); Total(s):
RUNTIME OF OVERALL TRANSIENT SIMULATION ( SECONDS ); Max. Df.(V) AND Avg. Df.(V): MAXIMUM AND AVERAGE DIFFERENCES COMPARED TO THE
SOLUTIONS OF ALL OUTPUT NODES PROVIDED BY IBM POWER GRID BENCHMARKS . SPDPtr : S PEEDUP OVER TR-FTS’ S Tran(s) IN TABLE II; SPDPr :
S PEEDUP OVER R-MATEX’ S Tran(s) IN TABLE II; Peak m: THE PEAK DIMENSION USED IN DR-MATEX FOR MEVP; Mem. Ratio over TR-FTS: T HE
PEAK MEMORY COMPARISON BETWEEN THE MAXIMUM MEMORY CONSUMPTION OF DR-MATEX OVER TR-FTS IN TABLE II.
Design
ibmpg1t
ibmpg2t
ibmpg3t
ibmpg4t
ibmpg5t
ibmpg6t
ibmpg1t new
ibmpg2t new
ibmpg3t new
ibmpg4t new
ibmpg5t new
ibmpg6t new
Average

Group #
100
100
100
15
100
100
100
100
100
15
100
100
—

Tran(s)
1.4
8.9
91.7
52.3
148.4
189.9
2.4
5.6
103.0
51.5
185.6
274.8
—

DR-MATEX
Total(s)
Max
1.9
11.4
129.9
112.2
178.9
234.2
2.8
7.0
140.9
108.4
227.8
317.7
—

Df.(V)
5.3e-5
4.6e-5
9.6e-5
9.9e-5
9.0e-5
3.4e-5
5.3e-5
4.6e-5
9.8e-5
9.9e-5
9.9e-5
3.4e-5
7.1e-5

Higham, Prof. Marlis Hochbruck, Junkai Jiang, Dr. John Loffeld, Dr. Jingwei Lu, Mulong Luo, Prof. Alexander Ostermann,
Prof. Mayya Tokman and Chicheng Zhang, for the informative
discussions. Hao Zhuang thanks the supports from UCSD’s
Powell Fellowship and Qualcomm FMA Fellowship. Last but
not the least, we thank for all the insightful comments from
the reviewers.

Avg Df.(V)
8.6e-6
8.6e-6
19.7e-6
27.9e-6
1.1e-6
7.2e-6
8.6e-6
8.6e-6
19.9e-6
27.6e-6
2.2e-6
7.1e-6
12.3e-6

Speedups
SPDPtr
SPDPr
4.0×
7.1×
4.5×
4.0×
2.9×
4.4×
8.8×
3.8×
3.2×
2.6×
3.8×
4.1×
21.8×
5.0×
61.6×
6.2×
25.6×
3.3×
98.0×
6.8×
27.3×
4.0×
25.5×
3.3×
26×
5×

Peak
m
6
7
6
8
7
7
6
7
7
8
7
7
6.7

Mem. Ratio
over TR-FTS
1.9
1.9
1.5
1.4
1.5
1.5
1.9
1.9
1.5
1.4
1.5
1.5
1.6

[13] Y. Zhang and C. Chu, “GDRouter: Interleaved global routing and
detailed routing for ultimate routability,” in Proc. IEEE/ACM Design
Autom. Conf., pp. 597–602, 2012.
[14] A. B. Kahng, S. Kang, H. Lee, I. L. Markov, and P. Thapar, “Highperformance gate sizing with a signoff timer,” in Proc. IEEE/ACM Int.
Conf. Comput.-Aided Design, pp. 450–457, 2013.
[15] C. Zhuo, G. Wilke, R. Chakraborty, A. Aydiner, S. Chakravarty, and W.K. Shih, “A silicon-validated methodology for power delivery modeling
and simulation,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided Design,
pp. 255–262, 2012.
[16] Z. Zeng, X. Ye, Z. Feng, and P. Li, “Tradeoff analysis and optimization
of power delivery networks with on-chip voltage regulation,” in Proc.
IEEE/ACM Design Autom. Conf., pp. 831–836, 2010.
[17] H. Zhuang, J. Lu, K. Samadi, Y. Du, and C. K. Cheng, “Performancedriven placement for design of rotation and right arithmetic shifters in
monolithic 3D ICs,” in Proc. IEEE Intl. Conf. Communications, Circuits
and Systems, vol. 2, pp. 509–513, 2013.
[18] S. K. Samal, K. Samadi, P. Kamal, Y. Du, and S. K. Lim, “Full chip
impact study of power delivery network designs in monolithic 3D ICs,”
in Proc. IEEE/ACM Int. Conf. Comput.-Aided Design, pp. 565–572,
2014.
[19] D. Xiang and K. Shen, “A thermal-driven test application scheme for
pre-bond and post-bond testing of three-dimensional ICs,” ACM Journal
on Emerging Technologies of Computing Systems, vol. 10, no. 2, p. Art.
no. 18, 2014.
[20] H. Zhuang, X. Wang, I. Kang, J.-H. Lin, and C. K. Cheng, “Dynamic
analysis of power delivery network with nonlinear components using
matrix exponential method,” in IEEE Intl. Symp. EMC&SI, 2015.
[21] H. Zhuang, W. Yu, I. Kang, X. Wang, and C. K. Cheng, “An algorithmic
framework for efficient large-scale circuit simulation using exponential
integrators,” in Proc. IEEE/ACM Design Autom. Conf., 2015.
[22] S. R. Nassif, “Power grid analysis benchmarks,” in Proc. IEEE/ACM
Asia South Pac. Design Autom. Conf., pp. 376–381, 2008.
[23] Z. Ye, Z. Zhu, and J. R. Phillips, “Sparse implicit projection (SIP)
for reduction of general many-terminal networks,” in Proc. Intl. Conf.
Computer-Aided Design, pp. 736–741, 2008.
[24] Z. Li, R. Balasubramanian, F. Liu, and S. Nassif, “2012 tau power grid
simulation contest: benchmark suite and results,” in Proc. IEEE/ACM
Int. Conf. Comput.-Aided Design, pp. 643–646, 2012.
[25] C. Zhuo, H. Gan, and W.-K. Shih, “Early-stage power grid design:
Extraction, modeling and optimization,” in Proc. IEEE/ACM Design
Autom. Conf., pp. 1–6, 2014.
[26] H. Zhuang, W. Yu, G. Hu, Z. Liu, and Z. Ye, “Fast floating random
walk algorithm for multi-dielectric capacitance extraction with numerical
characterization of Green’s functions,” in Proc. IEEE/ACM Asia South
Pac. Design Autom. Conf., pp. 377–382, 2012.
[27] C. Zhang and W. Yu, “Efficient space management techniques for largescale interconnect capacitance extraction with floating random walks,”
IEEE Trans. Comput.-Aided Design Integr. Circuits Syst, vol. 32, no. 10,
pp. 1633–1637, 2013.

R EFERENCES
[1] D. Kouroussis and F. N. Najm, “A static pattern-independent technique
for power grid voltage integrity verification,” in Proc. IEEE/ACM Design
Autom. Conf., pp. 99–104, 2003.
[2] S. R. Nassif and J. N. Kozhaya, “Fast power grid simulation,” in Proc.
IEEE/ACM Design Autom. Conf., pp. 156–161, 2000.
[3] M. S. Gupta, J. L. Oatley, R. Joseph, G.-Y. Wei, and D. M. Brooks,
“Understanding voltage variations in chip multiprocessors using a distributed power-delivery network,” in Proc. DATE, pp. 1–6, 2007.
[4] S. Lin, M. Nagata, K. Shimazake, K. Satoh, M. Sumita, H. Tsujikawa,
and A. T. Yang, “Full-chip vectorless dynamic power integrity analysis
and verification against 100uv/100ps-resolution measurement,” in Proc.
IEEE CICC, pp. 509–512, 2004.
[5] S. Lin and N. Chang, “Challenges in power-ground integrity,” in Proc.
IEEE/ACM Int. Conf. Comput.-Aided Design, pp. 651–654, 2001.
[6] R. Zhang, B. H. Meyer, W. Huang, K. Skadron, and M. R. Stan, “Some
limits of power delivery in the multicore era,” WEED, 2012.
[7] R. Zhang, K. Wang, B. H. Meyer, M. R. Stan, and K. Skadron,
“Architecture implications of pads as a scarce resource,” in Proc.
IEEE/ACM Intl Symp. Computer Architecture, pp. 373–384, 2014.
[8] K. Wang, B. H. Meyer, R. Zhang, K. Skadron, and M. R. Stan,
“Walking pads: Fast power-supply pad-placement optimization.,” in
Proc. IEEE/ACM Asia South Pac. Design Autom. Conf., vol. 20, p. 4,
2014.
[9] J. Lu, P. Chen, C.-C. Chang, L. Sha, D. Huang, C.-C. Teng, and C.K. Cheng, “ePlace: Electrostatics based placement using Nesterov’s
method,” in Proc. IEEE/ACM Design Autom. Conf., 2014.
[10] M. Pan, N. Viswanathan, and C. Chu, “An efficient and effective detailed
placement algorithm,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided
Design, pp. 48–55, 2005.
[11] J. Lu, H. Zhuang, P. Chen, H. Chang, C.-C. Chang, Y.-C. Wong,
L. Sha, D. Huang, Y. Luo, C.-C. Teng, and C. K. Cheng, “ePlaceMS: Electrostatics based placement for mixed-size circuits,” IEEE Trans.
Comput.-Aided Design Integr. Circuits Syst., vol. 34, no. 5, pp. 685–698,
2015.
[12] L. Xiao, Z. Xiao, Z. Qian, Y. Jiang, T. Huang, H. Tian, and E. F. Y.
Young, “Local clock skew minimization using blockage-aware mixed
tree-mesh clock network,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided
Design, pp. 458–462, 2010.

12

[28] W. Yu, H. Zhuang, C. Zhang, G. Hu, and Z. Liu, “RWCap: A floating
random walk solver for 3-D capacitance extraction of very-large-scale
integration interconnects,” IEEE Trans. Comput.-Aided Design Integr.
Circuits Syst, vol. 32, no. 3, pp. 353–366, 2013.
[29] L. O. Chua and P.-M. Lin, Computer Aided Analysis of Electric Circuits:
Algorithms and Computational Techniques. Prentice-Hall, 1975.
[30] Y. Saad, Iteravite Methods for Sparse Linear Systems. SIAM, 2003.
[31] T. A. Davis, Direct Method for Sparse Linear Systems. SIAM, 2006.
[32] T. Yu and M. D.-F. Wong, “PGT SOLVER: An efficient solver for power
grid transient analysis,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided
Design, pp. 647–652, 2012.
[33] J. Yang, Z. Li, Y. Cai, and Q. Zhou, “Powerrush: Efficient transient
simulation for power grid analysis,” in Proc. IEEE/ACM Int. Conf.
Comput.-Aided Design, pp. 653–659, 2012.
[34] X. Xiong and J. Wang, “Parallel forward and back substitution for
efficient power grid simulation,” in Proc. IEEE/ACM Int. Conf. Comput.Aided Design, pp. 660–663, 2012.
[35] L. Nagel, SPICE2: A computer program to simulate semiconductor
circuits. Ph.D. dissertation, 1975.
[36] S.-H. Weng, Q. Chen, and C. K. Cheng, “Time-domain analysis of
large-scale circuits by matrix exponential method with adaptive control,”
IEEE Trans. Comput.-Aided Design Integr. Circuits Syst, vol. 31, no. 8,
pp. 1180–1193, 2012.
[37] Y. Saad, “Analysis of some krylov subspace approximations to the matrix
exponential operator,” SIAM J. Numer. Anal., vol. 29, no. 1, pp. 209–228,
1992.
[38] M. Armbrust, A. Fox, R. Griffith, A. D. Joseph, R. Katz, A. Konwinski,
G. Lee, D. Patterson, A. Rabkin, I. Stoica, and M. Zaharia, “A view of
cloud computing,” Communications of the ACM, vol. 53, no. 4, pp. 50–
58, 2010.
[39] J. Dean and S. Ghemawat, “Mapreduce: simplified data processing on
large clusters,” CACM, vol. 51, no. 1, pp. 107–113, 2008.
[40] D. Xiang, Y. Zhang, and Y. Pan, “Practical deadlock-free fault-tolerant
routing in meshes based on the planar network fault model,” IEEE Trans.
Computers, vol. 58, no. 5, pp. 620–633, 2009.
[41] M. Zaharia, M. Chowdhury, M. J. Franklin, S. Shenker, and I. Stoica,
“Spark: cluster computing with working sets,” in HotCloud, pp. 10–10,
2010.
[42] Q. He, W. Au, A. Korobkov, and S. Venkateswaran, “Parallel power
grid analysis using distributed direct linear solver,” in IEEE Intl. Symp.
EMC, pp. 866–871, 2014.
[43] N. Gupte and J. Wang, “Secure power grid simulation on cloud,”
IEEE Trans. Comput.-Aided Design Integr. Circuits Syst, vol. 34, no. 3,
pp. 422–432, 2015.
[44] J. Wang and X. Xiong, “Scalable power grid transient analysis via
MOR-assisted time-domain simulations,” in Proc. IEEE/ACM Int. Conf.
Comput.-Aided Design, 2013.
[45] X.-X. Liu, H. Wang, and S. X.-D. Tan, “Parallel power grid analysis
using preconditioned gmres solver on CPU-GPU platforms,” in Proc.
IEEE/ACM Int. Conf. Comput.-Aided Design, pp. 561–568, IEEE, 2013.
[46] Z. Feng and P. Li, “Multigrid on GPU: tackling power grid analysis on
parallel simt platforms,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided
Design, pp. 647–654, 2008.
[47] H. Zhuang, S.-H. Weng, J.-H. Lin, and C. K. Cheng, “MATEX: A
distributed framework of transient simulation of power distribution
networks,” in Proc. IEEE/ACM Design Autom. Conf., 2014.
[48] M. J. Gander and S. Guttel, “PARAEXP: A parallel integrator for linear
initial-value problems,” SIAM Journal on Scientific Computing, vol. 35,
no. 2, pp. C123–C142, 2013.
[49] S.-H. Weng, Q. Chen, N. Wong, and C. K. Cheng, “Circuit simulation
via matrix exponential method for stiffness handling and parallel processing,” in Proc. IEEE/ACM Int. Conf. Comput.-Aided Design, pp. 407–
414, 2012.
[50] J. van den Eshof and M. Hochbruck, “Preconditioning Lanczos approximations to the matrix exponential,” SIAM J. Sci. Comput., vol. 27, no. 4,
pp. 1438–1457, 2006.
[51] M. A. Botchev, “A short guide to exponential krylov subspace time
integration for maxwell’s equations,” Dept. of Applied Mathematics,
Univ. of Twente, 2012.
[52] H. Zhuang, S.-H. Weng, and C. K. Cheng, “Power grid simulation using
matrix exponential method with rational krylov subspaces,” in Proc.
IEEE Intl. Conf. ASIC, 2013.
[53] Q. Chen, S.-H. Weng, and C. K. Cheng, “A practical regularization
technique for modified nodal analysis in large-scale time-domain circuit
simulation,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst,
vol. 31, no. 7, pp. 1031–1040, 2012.

[54] J. Wilkinson, “Kronecker’s canonical form and the QZ algorithm,”
Linear Algebra and its Applications, vol. 28, pp. 285–303, 1979.
[55] P. Li, “Parallel circuit simulation: A historical perspective and recent developments,” Foundations and Trends in Electronic Design Automation,
vol. 5, no. 4, pp. 211–318, 2012.
[56] X. Ye, M. Zhao, R. Panda, P. Li, and J. Hu, “Accelerating clock
mesh simulation using matrix-level macromodels and dynamic time step
rounding,” in ISQED, pp. 627–632, 2008.
[57] B. Hindman, A. Konwinski, M. Zaharia, A. Ghodsi, A. D. Joseph,
R. Katz, S. Shenker, and I. Stoica, “Mesos: A platform for fine-grained
resource sharing in the data center,” in NSDI, pp. 22–22, 2011.
[58] X. Chen, Y. Wang, and H. Yang, “NICSLU: An adaptive sparse matrix
solver for parallel circuit simulation,” IEEE Trans. Comput.-Aided
Design Integr. Circuits Syst, vol. 32, no. 2, pp. 261–274, 2013.
[59] M. Hochbruck and A. Ostermann, “Exponential integrators,” Acta Numerica, vol. 19, pp. 209–286, 2010.
[60] M. Hochbruck, C. Lubich, and H. Selhofer, “Exponential integrators for
large systems of differential equations,” SIAM J. Sci. Comput., vol. 19,
no. 5, pp. 1552–1574, 1998.

Hao Zhuang (S’11) is a Ph.D. candidate at the
Department of Computer Science and Engineering,
University of California, San Diego, CA, USA (UC
San Diego). He received his C.Phil. degree in Computer Science from UC San Diego in June, 2015.
His current research interests include numerical and
optimization algorithms, with the applications in
design automation, very large-scale integration systems, network and data analysis.
He has industrial experiences in several companies, including Synopsys, Inc., Mountain View, CA,
USA, and Ansys, Inc., San Jose, CA, USA, where he worked on the largescale circuit analysis and dynamic power network simulation algorithms
within several products. He was one of the main software developers of
RWCap, a parallel program for VLSI capacitance extraction using floating
random walk algorithm at Tsinghua University. At UC San Diego, he designed
the numerical algorithms for large-scale circuit simulation using matrix
exponentials, and optimization algorithms for electrostatics based VLSI global
placement.
Mr. Zhuang was the recipient of the Charles Lee Powell Fellowship and
the Qualcomm FMA Fellowship. He is a student member of IEEE, ACM and
SIAM. He serves as a reviewer for IEEE Transactions on Computer Aided
Design (TCAD), and an external reviewer for Design Automation Conference
(DAC) and International Symposium on Physical Design (ISPD).

Wenjian Yu (S’01-M’04-SM’10) received the B.S.
and Ph.D. degrees in computer science from Tsinghua University, Beijing, China, in 1999 and 2003,
respectively.
In 2003, he joined Tsinghua University, where
he is an Associate Professor with the Department
of Computer Science and Technology. He was a
Visiting Scholar with the Department of Computer
Science and Engineering, University of California,
San Diego, La Jolla, CA, USA, twice during the
period from 2005 to 2008. He serves as an associate
editor of IEEE Transactions on Computer Aided Design since 2016. His
current research interests include numerical modeling and simulation techniques for designing complex systems (like integrated circuit, touchscreen,
cyber physical system, etc.), and the matrix computation algorithms for BIG
DATA analysis. He has authored three books and over 130 papers in refereed
journals and conferences.
Dr. Yu was a recipient of the Distinguished Ph.D. Award from Tsinghua
University in 2003 and the Excellent Young Scholar Award from the National
Science Foundation of China in 2014.

13

Shih-Hung Weng is currently with Facebook as
Research Scientist and working on large scale distributed system for payments services. He received
the B.S. and M.S degrees in Computer Science from
National Tsing Hua University, Hsinchu, Taiwan, in
2006 and 2008, and the Ph.D. degree from University of California San Diego, La Jolla, U.S.A in 2013
under the supervision of Dr. Chung-Kuan Cheng.
His thesis focused on parallel circuit simulation
for power-grid noise analysis in very large scale
integration (VLSI) system.

Chung-Kuan Cheng (S’82-M’84-SM’95-F’00) received the B.S. and M.S. degrees in electrical engineering from National Taiwan University, and the
Ph.D. degree in electrical engineering and computer
sciences from University of California, Berkeley in
1984.
From 1984 to 1986 he was a senior CAD engineer
at Advanced Micro Devices Inc. In 1986, he joined
the University of California, San Diego, where he is
a Distinguished Professor in the Computer Science
and Engineering Department, an Adjunct Professor
in the Electrical and Computer Engineering Department. He served as a
principal engineer at Mentor Graphics in 1999. He was an associate editor
of IEEE Transactions on Computer Aided Design for 1994-2003. He is a
recipient of the best paper awards, IEEE Trans. on Computer-Aided Design
in 1997, and in 2002, the NCR excellence in teaching award, School of
Engineering, UCSD in 1991, IBM Faculty Awards in 2004, 2006, and
2007, the Distinguished Faculty Certificate of Achievement, UJIMA Network,
UCSD in 2013. He is appointed as an Honorary Guest Professor of Tsinghua
University 2002-2008, and a Visiting Professor of National Taiwan University
2011, and 2015.
His research interests include medical modeling and analysis, network
optimization and design automation on microelectronic circuits.

Ilgweon Kang (S’12) received the B.S and M.S.
degrees in electrical and electronic engineering from
Yonsei University, Seoul, Korea, in 2006 and 2008,
respectively. He is currently pursuing the Ph.D
degree at University of California at San Diego
(UCSD). He was with the Research and Development Division, SK Hynix, Icheon, Korea, from 2008
to 2012, where he was involved in development of
3D stacked DRAM with (Through Silicon-Via) TSV
technology.
His current research interests include 3D IC design and optimization, CAD, and design for testability

Jeng-Hau Lin (S’15) is a Ph.D student at Department of Computer Science and Engineering, University of California, San Diego, CA, USA. He received
his B.S degree in Electrical Engineering in 2005 and
M.S. degree in Communication Engineering in 2007
from National Taiwan University.
His research interests includes time-frequency
analyses on signal integrity and the consensus of
multi-agent system.

Xiang Zhang (M’13) received the B.Eng. in Information Engineering from Shanghai Jiao Tong
University (SJTU), China and M.S. in Electrical and
Computer Engineering from University of Arizona,
Tucson, AZ, in 2010.
From 2011 to 2014, he was a Senior Hardware Engineer with Qualcomm Technologies, Inc. In 2014,
he joined Apple Inc as an iPhone Hardware System
Engineer. Since 2012, he has been working towards
the Ph.D. degree in Department of Electrical and
Computer Engineering at University of California,
San Diego, CA, USA.
His current research interests include power distribution network design
and optimization, circuit simulation and system-on-chip design.

Ryan Coutts received a B.S. in electrical engineering from the University of California, Riverside
in 2005. He then received his M.S. in electrical
engineering from Stanford University in 2006. After
receiving his M.S. degree, he worked as a signal and
power integrity engineer at NVIDIA Inc. specializing in FSB, DDR and power integrity. He is currently
working at Qualcomm Inc. where he works in low
power design methodologies for mobile covering
the range from power integrity, thermal optimization
and on-chip timing. Mr. Coutts has filed 13 patents
related to his innovation as well and is currently pursuing a Ph.D at the
University of California, San Diego.

14

