{Simulation Algorithms with Exponential Integration for Time-Domain Analysis of Large-Scale Power Delivery Networks
Hao~Zhuang
\IEEEmembership{Student Member
IEEE
Hao~Zhuang
\IEEEmembership{Student Member
IEEE
%\input{new_abst.tex \input{abst_tcad.tex
abst_tcad.tex
We design an algorithmic framework using matrix exponentials for time-domain simulation of power delivery network (PDN). %By the explicit formulation via exponential time integration, 
Our framework can reuse factorized matrices to simulate the large-scale linear PDN system with variable step-sizes. In contrast, current conventional PDN simulation solvers have to use
% implicit numerical integration schemes with 
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 \emph{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 \emph{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 \emph{Embarrassingly Parallel}} model. 

Experimental results show R-MATEX and DR-MATEX can achieve up to around  $14.4\times$ and $98.0\times$  runtime speedups over traditional trapezoidal integration based solver with fixed time-step approach.  
  
\begin{comment}


In this work, we design an efficient and accurate algorithmic framework using matrix exponentials for time-domain simulation of power  delivery network (PDN). Thanks to the  explicit exponential time integration scheme with high order approximation of differential equation system, our framework can reuse factorized matrices for adaptive time stepping without loss of accuracy. The key operation of matrix exponential and vector product (MEVP) is computed by proposed efficient rational Krylov subspace method and helps achieve large stepping.  

With the enhancing capability of time marching and high-order approximation capability, we design R-MATEX, which outperforms the classical PDN simulation method using trapezoidal formulation with fixed step size (TR-FTS). We also  propose a distributed computing framework, DR-MATEX, and highly  accelerate the simulation speedup by reducing Krylov subspace generations caused by frequent breakpoints from the side of current sources. 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 current sources. Then, the subtasks are  processed in parallel at different computing nodes and summed up at the end of simulation to provide the accurate solutions. The experimental results show R-MATEX and DR-MATEX can achieve $11.4\times$ and $68.0\times$ runtime speedups on average over TR-FTS.  





\emph{MATEX} is a transient circuit simulation framework for
 power distribution networks (PDNs), which utilizes \underline{mat}rix \underline{ex}ponential kernel with Krylov subspace approximations to solve differential equations of linear circuit.
 MATEX is also a distributed computing framework.
First, the whole simulation task is divided into subtasks 
based on decompositions of current sources, in order to reduce
the frequency of Krylov subspace generations. 
Then these subtasks are distributed to different computing nodes and 
processed in parallel.
%parallely processed. 
Within each node, after the  matrix factorization at the beginning of simulation, the adaptive time stepping solver is performed without extra matrix re-factorizations.
% and with less Krylov subspace generations. 
MATEX overcomes the stiffness hinder of previous 
matrix exponential-based circuit simulator 
by \emph{rational} Krylov subspace method, which leads to larger step sizes with  smaller dimensions of Krylov subspace bases and highly accelerates the whole computation.
In this work, we also discuss the matrix exponential integrators and how to combined it with low order scheme to further accelerate simulation without loss of accuracy.
\end{comment}

 
%This greatly reduces both memory and computational complexities.
%Since the integration kernel an inherently high polynomial approximation, 
%MATEX can utilize larger time steps to do time marching than
%circuit simulation algorithms using traditional integration kernels, 
%e.g. trapezoidal, backward Euler.
%During the simulation, MATEX has the capability of deploying adaptive time stepping 
%and still reusing factorized matrices at the very beginning, 
% even with different time step sizes,
%MATEX outperforms both traditional fixed and adaptive time stepping methods,
%Besides, by utilizing superpositions and
%distributed parallel processing, 
%the scalablity of MATEX reduces runtimes furthermore, 
%e.g., achieving around 13X over trapezoidal method with fixed time step in transient computing on IBM power grid benchmarks.


ack.tex
%We acknowledge the support from NSF-CCF 1017864, NSFC 61422402 and NSFC 61274033.  
We thank Prof. Scott B. Baden,  Prof. Mike Botchev, Dr. Quan Chen, Dr. Peng Du,  Dr. Martin Gander, Prof. Nicholas 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.
%Dr. Steve McCormick, Dr. Shen Lin, Dr. Norman Chang,
  
adaptive_step_control.tex

The new proposed MEXP for the power grid simulation can significantly benefit from the adaptive step size control because 
the rational Krylov subspace still preserves the capability of adaptivity, which has been demonstrated the effectiveness 
in linear circuits~\cite{Weng12_TCAD}. And, with the relaxation of the stiffness constraint, MEXP can now utilize even larger 
step size, e.g., 1\emph{ns}, to perform transient analysis. Such large step size is particularly useful for the power grid 
simulation where the voltage noise usually includes the low frequency behavior of both package and board levels. 

The idea of our adaptive step scheme is to step forward as much as possible so that MEXP can finish the simulation 
quickly. With the insight from Eqn. (\ref{eq:scale}), MEXP is able to adjust $\alpha$ to calculate result of arbitrary 
step size with only performing one Arnoldi process. However, because MEXP has to guarantee slopes of input sources must 
be constant within the step size, the maximum allowed step size is limited. Therefore, MEXP will first compute the rational 
Krylov subspace using $\mbf{x}(t)$, and then results within the maximum allowed step size, which is determined by input 
sources, can be calculated through the step size scaling. 

Algorithm~\ref{algo:adaptive_step} shows MEXP with adaptive step control scheme. In order to comply with the required 
accuracy during the simulation, the allowed error $err_k$ at certain time instant $t_k$ is defined as 
\begin{eqnarray}
    err_k \le \frac{E_{Tol}}{T}\times h_k, \nonumber
\end{eqnarray}
where $E_{Tol}$ is the error tolerance in the whole simulation process, $T$ is the simulation time span, $h_k$ is the 
step size at time $t_k$, and $err_k$ is the posterior error of MEXP from Eqn. (\ref{eq:post_err_est}). Hence, when 
we construct the rational Krylov subspace at the current time, we will increase $m$ until the $err_k$ satisfies the 
error tolerance. 

\begin{algorithm} 
    \label{algo:adaptive_step}
    \caption{MEXP with Adaptive Step Control}
    \KwIn{$\mbf{C}$, $\mbf{G}$, $\mbf{B}$, $\mbf{u}(t)$, $\tau$, error tolerance $E_{Tol}$ and simulation time $T$}
    \KwOut{$\mbf{x}(t)$}
    {
        $t = 0$; 
	$\mbf{x}(0)$ = DC\_analysis\;
        %Choose $\tau$ for rational Krylov subspace method\;
        $[\mbf{L}_{sub},\mbf{U}_{sub}] = $LU($\mbf{C} +\gamma \mbf{G}$)\;
        \While{$t \leq T$}
        {	
	    %Update $v_{k+1}$ with $\mbf{L}, \mbf{U}, \mbf{W}$ by (\ref{eq:blu_f})\;
	    %and update $\mathbf{L},\mathbf{U}$ by (\ref{eq:blu})\;
	    Compute maximum allowed step size $h$ from $\mbf{u}(t)$\;
	    Construct $\mbf{H_m}$, $\mbf{V_m}$ by Arnoldi process until $err_k \le \frac{E_{Tol}}{T}h$\;

	    $\alpha = \frac{h}{\gamma}$ \;
	    Compute $\mbf{x}(t+h)$ by (\ref{eq:scale})\;

	    $t = t+h$; 
	}
    }
\end{algorithm}
The complexity of MEXP with adaptive step size is mainly determined by the total number of required backward/forward
substitutions during the simulation process. This is because the dimension of $\mbf{H_m}$ is small, e.g., around 
$5 \sim 20$, and the computation overhead of the step size scaling is negligible. The number of total substitution 
operations is 
\begin{eqnarray}
    \sum_{i=0}^N m_{i}, \nonumber
\end{eqnarray}
where $N$ is total time steps, and $m_i$ is required dimension of rational Krylov subspace of each time step. Compared to 
the trapezoidal method where the number of substitution operations depends on the fixed step size, MEXP could use less 
substitution operations as long as the maximum allowed step size is much larger than the fixed step size. Our experiments 
in the following section demonstrates it is usually the case for the power grid simulation.

bio.tex
%\newpage
\begin{IEEEbiography}[{\includegraphics[width=1in,clip,keepaspectratio]{./bios/zhuang_bio_photo.eps}}]{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 large-scale 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). 
\end{IEEEbiography}
%\vspace{-1.2cm} 
\begin{IEEEbiography}[{\includegraphics[width=1in,clip,keepaspectratio]{./bios/WenjianYU.eps}}]
{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.
\end{IEEEbiography} 
\begin{IEEEbiography}[{\includegraphics[width=1in,clip,keepaspectratio]{./bios/ShihHungWeng.eps}}]{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. 
\end{IEEEbiography}
%\vspace{-1.2cm}
\begin{IEEEbiography}[{\includegraphics[width=1in,clip,keepaspectratio]{./bios/ilgweon_Kang.eps}}]{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
\end{IEEEbiography} 
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{./bios/JengHauLin02.eps}}]{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.  
\end{IEEEbiography}
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{./bios/xiang_photo.eps}}]{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.
\end{IEEEbiography} 
\begin{IEEEbiography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{./bios/Ryan.eps}}]{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.
\end{IEEEbiography} 
\begin{IEEEbiography}[{\includegraphics[height=1.1in,clip,keepaspectratio]{./bios/ck2010.eps}}]{Chung-Kuan Cheng} (S'82-M'84-SM'95-F'00)  
\begin{comment}
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.
\end{comment}
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.
\end{IEEEbiography}  
%\begin{IEEEbiography}[]{}
%\end{IEEEbiography}
%\section{Appendix}
%\label{sec:append}
% \input{append}
% that's all folks

copyright.tex
\noindent 
Copyright (c) 2015 IEEE. Personal use of this material is permitted. However, permission to use
this material for any other purposes must be obtained from the IEEE by sending an email to
pubs-permissions@ieee.org.
\newpage 
drmatex.tex
\subsubsection{Motivation}
With the ability of large step time-marching, the input term $\mbf b$ embedded in Eq. (\ref{eqn:discrete_sol}) actually serves a double-edged sword in MEVP computation. There are usually many input currents in PDNs as well as their breakpoints, which narrow the regions for the time stepping of matrix exponential-based method due to the unaligned breakpoints. The region before the next transition $t_s$ may be shortened when there are a lot of independent input sources injected into the linear system. It leads to more chances of generating new Krylov subspace. We want to reduce the number of  Krylov subspace generations to improve the runtime performance. \footnote{%\bf 
The breakpoints also put the constraint on TR-FTS and other low order integration schemes (Backward Euler, etc.). However, the incapability of large step time-marching and theoretical limitations of low order accuracy refrain them from reaching this problem in the first place.}
\vspace{+0.1in}

\subsubsection{Treatment and Methodology} 
\label{ltsgts}
% In order to solve the frequent breakpoints problem by distributed computing, 
In matrix exponential integration framework, the flexible time stepping can choose any time spot until the next input transition spot $t_s$, as long as the 
%matrix exponential and vector 
approximation of $e^{\mbf A}\mbf v$
is accurate enough. The Krylov subspace can be reused 
when $t+h \in [t,t_s]$, only by scaling $\mbf H_m$ with $h$ 
in 
$$\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V_m} e^{h \mbf H_m} \mbf e_1-\mbf P(t,h).$$
This is an important feature that even doing the adaptive time stepping, we can still use the last Krylov subspaces. 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 for  breakpoints  of input sources:
\vspace{+0.08in}
\begin{itemize}
\item 
%\textbf{Definition}: 
\emph{Local Transition Spot} ($LTS$): the set of $TS$ at an input source to the PDN.
\item
%  \textbf{Definition}: 
\emph{{Global Transition Spot}} ($GTS$): the union of $LTS$ among all the input sources to the PDN. 
\item  
%\textbf{Definition}: 
\emph{{Snapshot}}:  a set $GTS \setminus LTS$ at {an input source}. 
\end{itemize}
\vspace{+0.08in}

If we simulate the PDN with respective to all the inputs,
$GTS$ are the places where generations of Krylov subspace 
cannot be avoided. 
For example, there are three input sources in a PDN (Fig. \ref{fig:rc_ckt}). 
The input waveforms are shown in Fig. \ref{fig:input}.
The first line is $GTS$, which is contributed by the union of $LTS$ 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 $LTS$ and keeps track of $Snapshot$ for the later usage of summation via superposition. {\color{red} Between two LTS points $t$ and $t+h$,  the \emph{Snapshot} points $t+h_1 < t+h_2 < \cdots < t+h_l  \in (t, t+h]$ can reuse the Krylov subspace generated at $t$}, which is mentioned in Sec. \ref{sec:discuss}. For each node, {the chances of generation of Krylov subspaces are reduced} and the time periods of reusing these subspaces are enlarged locally, which bring huge computational benefit when processing these subtasks in parallel. {\color{red} Besides, when subtasks are assigned, there are no communication among those computing nodes, leading to so-called ``Embarrassingly Parallel" computation model.}
 
Above, we divide the simulation task by input sources. We can, more aggressively, decompose the  task according to the ``bump'' shapes within such input pulse sources.  We group the input sources which have the same (t\_delay, t\_rise, t\_fall, t\_width) into one set. 
For example, the input source \#1 of Fig. \ref{fig:input} is divided to  \#1.1 and \#1.2 in \ref{fig:ckt_ts_decomp}. The input source \#2 in  Fig. \ref{fig:ckt_ts_decomp} is divided to  \#2.1 and \#2.2 in Fig. \ref{fig:ckt_ts_decomp}.  Therefore, there are four groups in Fig. \ref{fig:ckt_ts_decomp}, Group 1 contains $LTS\#1.1$, Group 2 contains $LTS\#2.1$, Group 3 contains $LTS \#2.2$, and Group 4 contains $LTS \#1.2$ and $\#3$. Our proposed framework MATEX is shown in Fig. \ref{fig:drmatex}. After pre-computing $GTS$ and decomposing $LTS$ based on ``bump'' shape (Fig. \ref{fig:ckt_ts_decomp}), we group them and form $LTS$ $\#1 \sim \#K$.\footnote{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}  

MATEX scheduler sends out $GTS$ and $LTS$ to different MATEX slave node.
Then the simulations are  performed in parallel.
There are still no communications among nodes before the ``write back'' (in Fig. \ref{fig:drmatex}), by when all nodes complete their transient simulations.
Within each slave node, ``circuit solver'' (Alg. \ref{algo:ckt_solver}) 
computes transient response with varied time steps. 
Solutions are obtained without re-factorizing matrix during the transient computing. 
After finishing all simulations from slave nodes, they writes back the results
and informs the MATEX scheduler.
\begin{figure*}[ht]
\centering{
\includegraphics[width=0.8\textwidth]%{./figs/hierachical_ckt_sim.eps}
{./figs/hierachical_ckt_sim_new.eps}
%\vspace{-0.2in}
\centering{
\caption{DR-MATEX: the distributed MATEX framework using R-MATEX circuit solver}
}
\label{fig:drmatex}
}
\end{figure*}
%\textbf{Why does MATEX can leverage such distributed computing framework, but trditional methods cannot?}

\begin{figure}[ht] 
    \centering
    \includegraphics[width=3.0in]{./figs/rc_ckt.eps}
     \caption{Part of a PDN model with input sources from Fig. \ref{fig:input}.}
	\label{fig:rc_ckt}
\end{figure}

\begin{figure}[ht] 
%\scriptsize
    \centering
    \includegraphics[width=3.7in]%{./figs/time_spots.eps}
    {./figs/time_spots_2.eps}
      %  \vspace{-0.15in}
    \caption{
	{ 
	{Illustration of input transitions. 
	GTS: Global Transition Spots; LTS: Local Transition Spots;
	{\bf Snapshots: the crossing positions by dash lines and LTS \#$k$ without solid points.} 
	%that is the snapshots at LTS \#$k$. % which is the set of $GTS \setminus LTS$. 
	}
	}
    }
    \label{fig:input}
\end{figure}
%\vspace{-0.03in}


\begin{figure}[ht]
%\scriptsize
    \centering
    \includegraphics[width=3.65in]%{./figs/ckt_ts_decomp.eps}
    {./figs/source_large.eps}
    %   \vspace{-0.15in}
    \caption{ 
	{Grouping of ``Bump'' shape transitions for sub-task simulation.
	The matrix exponential-based method can  utilize  adaptive stepping in each LTS and reuse the Krylov subspace generated at the most recent solid point. 
	However, traditional methods (TR, BE, etc) still need to do time marching step by step, either by pairs of backward and forward substitutions and proceed with fixed time step, or re-factorize matrix and solve linear system when using adaptive stepping. 
	    (Pulse input information:
	        t\_delay: initial delay time;
	        t\_rise: rise time;
	        t\_width: width of pulse-wise; 
	        t\_fall: fall time; and
	        t\_period: period).
	}
    }
	\label{fig:ckt_ts_decomp}
\end{figure}

\begin{algorithm}
\label{algo:ckt_solver} 
{ 
\caption{DR-MATEX: the distributed NATEX framework using R-MATEX} 
\KwIn{ $LTS$ \#$k$,  $GTS$, $\mbf X_1$, $\mbf X_2$, and $\mbf P_k$, $\mbf F_k$, which contain the corresponding $\mbf b$ for node $k$; error
tolerance $E_{tol}$, and simulation time $T$.}
%\KwOut{ $\mbf V_m, \mbf H_m,  m$ }
\KwOut{Local solution $\mbf x$ along $GTS$ in node 
$k \in [1,\cdots,S]$, where $S$ is the number of nodes}
    {   
%	[$\mbf L, \mbf U$]=LU($\mbf X_1$)\\
%	Normalize the intial vector 
	$t=T_{start}$\;
	%$s=T_{start}$\;
	$\mbf x(t) = \text{Local\_Initial\_Solution}$\;
	%$\mbf x(t) = local\_DC\_analysis$; 
%	{\tcp*{ DC for input at \#$s$}} 
	%$[\mbf L,\mbf U] = \text{LU\_Decompose}( \mbf C+\gamma \mbf G)$\;
	$[\mbf L,\mbf U] = \text{LU\_Decompose}( \mbf X_1)$\;
	\While { $t \leq T_{end}$}
	{
	Compute maximum allowed step size $h$ based on $GTS$\;
	% and	$ \alpha = \frac{2h}{\gamma}$\; 
	 {
    %\tcp{\color{blue} R-MATEX mode: acurate large step marching}
	\If{$t \in $  LTS \#$k$}
	{
	\tcc{Generate the Krylov subspace for the time point of $LTS$ and compute $\mbf x$}
%	\tcp*{Reach LTS and update $\mbf V_m, \mbf H_m, \mbf v, s$}
	$[\mbf x(t+h), \mbf V_m, \mbf H_m,\mbf v] = \text{MATEX\_Arnoldi}(\mbf L,\mbf U,\mbf X_2,$ $ h, t, \mbf x(t), \epsilon, \mbf P_k(t,h), \mbf F_k(t,h))$\;
	%\label{algo2:arnoldi}
%	Update $\mbf x(t +h) $ based on Eq. (\ref{eqn:new_exact})\; 
	$a_{lts} = t$\; 
	%\tcp{put an LTS anchor for current $t$} 
	}
	\Else{
	\tcc{Compute the $\mbf x$ at $Snapshot$ by reusing the latest Krylov subspace}
%	    \For  {$s_i \in (t_i + h,t_{i+1}]$}
	    {
	    %// Take snapshots by reusing $\mbf V_m , \mbf H_m, \mbf v$ \;
	    $h_{a} = t + h - a_{lts}$\;
	    $\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V_m} e^{h_{a} \mbf H_m} \mbf e_1 - \mbf P_k(t,h)$\;
	   % \label{algo2:h}
%	    Update $\mbf x(t +h) $ based on Eq. (\ref{eqn:new_exact})\; 
	%$\lVert r_m(h) \rVert  =  \lVert \mbf v \rVert \left| h_{m+1,m} e_m^T e^{h\mbf H_m} e_1 \right|$
	    }
	}
	$t= t+h$\;
	}
}    }
}
\end{algorithm}

\subsection{Complexity Analysis} 
\label{sec:complexity} 
Suppose on average we have Krylov subspace basis dimension $m$ at each time step along the time span,
one pair of backward and forward substitutions has time complexity $T_{bs}$. 
The total time complexity of serial parts is $T_{serial}$, which includes matrix factorizations, local results collection, etc.
The matrix exponential evaluation using $\mbf H_m$ is $T_{H}$ which costs time complexity $O(m^3)$, plus extra $T_e$ to form $\mbf x$, which costs $O(n m^2)$.
%$T$, $T_s$ and $T_e$ are positive proportional to the nonzero terms in operated matrices.
Given $K$ points of $GTS$, without decomposition of input transitions,
the time complexity is 
$KmT_{bs} + K (T_{H}+ T_e)  + T_{serial}$.
%$O(KmT + K (m^3+ T_e)  + T_s)$, 
After dividing the input transitions and sending to enough computing nodes, 
we have $k$ points of $LTS$ for each node based on feature extraction and grouping 
(e.g., $k=4$ for one ``bump'' shape feature).
%suppose we have enough computing nodes to parallel process all input sets, 
The total computation complexity is 
$kmT_{bs} + K (T_{H}+ T_e)  + T_{serial}$, where  $K (T_{H}+ T_e)$ contains the portion of computing $Snapshot$.
%$O( kmT + K(m^3 + T_e) + T_s)$.
%Since no commnication required between the different processors, the extra complexity form
The speedup of distributed computing over single MATEX, 
\begin{eqnarray}
%\footnotesize 
%\text{Speedup}=\frac{ O(Km^3)+ O(Ktm)+O_{s} }  {O(Km^3)+O(ktm)+O_{s}}
\text{Speedup}=\frac{ KmT_{bs} + K(T_{H}+ T_e)+ T_{serial}}  
{ kmT_{bs}+ K(T_{H}+ T_e)+T_{serial}}.
%\text{Speedup}=\frac{ KmT + K(m^3+ T_e)+ T_{s} }  { kmT+ K(m^3+T_e)+T_{s}}
\end{eqnarray}
In R-MATEX, we have very small $m$. Besides, $T_{bs}$ is larger than $T_{H}+T_e$. 
Therefore, the most dominating part is the $KmT_{bs}$.
We can always decompose input transitions, and make $k$ very small compared to $K$.
Traditional method with fixed step size has $N$ steps for the  entire simulation.
The complexity is $NT_{bs}+T_{serial}$. Then the speedup of distributed MATEX over the ones with fixed step size is,
\begin{eqnarray}
\label{eqn:tr_spd}
%\text{Speedup}=\frac{ O(Km^3)+ O(Ktm)+O_{s} }  {O(Km^3)+O(ktm)+O_{s}}
\text{Speedup}'=\frac{NT_{bs}+T_{serial} }  { kmT_{bs}+ K(T_{H}+T_e)+T_{serial}}.
\end{eqnarray}
Usually, $N$ is much larger $K$ and $km$. 
Uniform step sizes make $N$ increase due to resolution of input transitions, 
to which $K$ is not so sensitive.
As mentioned before, $k$ can be maintained in a small number.
When elongating time span of simulation,
$N$ will increase,
however, $k$ will not change due to its irrelevance to time span 
(may bring more input transition features and increase computing nodes),
and speedups tend to become larger.
Therefore, our MATEX is more robust and has promising theoretical speedups.
\vspace{+0.2in}
exp_results.tex
We implement all the algorithms % our proposed MATEX 
in MATLAB R2014b\footnote{Measurements reported are on MATLAB implementations. They are subject to limitations and are not directly comparable to C++ implementations reported in literature such as \cite{Wang13}.} and use UMFPACK package for LU factorization. 
%First, we compare the different Krylov methods (MEXP, I-MATEX, and R-MATEX) to show rational Krylov subspace is the best for linear PDN circuit simulation in Table. \ref{tab:rc_ckt}. 
%Second, we test our HR-MATEX (MATEX with low order scheme TR),
First, we compare I-MATEX, R-MATEX and TR in order to show our runtime improvements in single machine framework in Table \ref{tab:perf_mexp}. 
%It shows HR-MATEX is a treatment for single computing node to ease the small time-step constraint for MATEX-like PDN circuit simulation. 
Second, we show our distributed framework DR-MATEX achieves large speedups in Table \ref{tab:pg_rmatex}. 
The experiments are conducted on the server with Intel(R) Xeon (R)~E5-2640 v3 2.60GHz processor and 125GB memory.  

%\subsubsection{Simple $\mbf C$ matrix cases}

%$\mbf C$ is invertible and has a small range of capacitance values. 
%In the following parts, we apply our I-MATEX and R-MATEX to futher cases and make comparisions 
%with other methods.

%but it spread  more stable along the total simulation in R-MATEX and $m_p$ is smaller.


%\vspace{-0.02in}
%\subsection{Power Distribution Network Simulation}


\begin{comment}
\subsubsection{Industrial Cases (Impulse Response)} \

We utilize industrial power distribution network cases to show MATEX's adaptivity.
The rationale is that designers rely on the simulation result of impulse response of the 
PDN design in order to characterize the design. 
Many previous works \cite{Du10}\cite{Tian96}\cite{Hu12_TVLSI}
have proposed different PDN analysis based on the impulse 
response. 
The nature of impulse response of the PDN design, 
which contains low-, mid- and high-frequency components, 
can significantly enjoy the adaptive step size in MATEX.  
We would also like to mention that the impulse response 
based analysis is not only for the PDN design, 
but also for worst-case eye opening analysis in the high speed interconnect
\cite{Analui05}\cite{Shi08}. 

The industrial PDN design are given, including
on-chip power grid, package and board. 
The power grid consists of four metal 
layers: M1, M3, M6 and Redistribution Layer (RDL). 
The physical parameters of each metal layer 
is listed in Table~\ref{tab:spec_metal_layers}. 
The package is modeled as an RL series at each C4 bump, 
and the board is modeled as a lumped RLC network. The
specification of each PDN design is listed in Table~\ref{tab:spec_pdn_designs} 
where the size of each design 
ranges from $45.7$K to $7.40$M. 
\begin{table}[htp]
  \caption{Widths and pitches of metal layers in the PDN design($\mu m$).}
  \centering
  \begin{tabular}{cccccccc}
    \multicolumn{2}{c}{M1} &\multicolumn{2}{c}{M3} 
    &\multicolumn{2}{c}{M6} &\multicolumn{2}{c}{RDL}\\
    \hline pitch &width  &pitch &width  &pitch &width  &pitch &width \\
    \hline 2.5   &0.2  &8.5  &0.25  &30 &4 &400 &30\\
    \hline 
  \end{tabular}
  \label{tab:spec_metal_layers}
\end{table}

\begin{table}[htp]
    \caption{Specifications of PDN Designs}
    \centering
    \begin{tabular}{|c||r|r|r|r|r|}
    \hline
	Design  & Area ($mm^2$)	& \# R	    & \# C	& \# L	    & \# Nodes  \\ \hline
	D1	& $0.35^2$	& 23221	    & 15193	& 15193	    & 45.7K	\\ \hline
	D2	& $1.40^2$  	& 348582    & 228952	& 228952    & 688K	\\ \hline
	D3	& $2.80^2$  	& 1468863   & 965540	& 965540    & 2.90M	\\ \hline
	D4	& $5.00^2$  	& 3748974   & 2467400	& 2464819   & 7.40M	\\ 
    \hline
    \end{tabular}
    \label{tab:spec_pdn_designs}
\end{table}
The impulse response can be derived from the simulation result 
of a step input from $0$V to $1$V with a small transition time. 
Hence, we inject a step input to each PDN design and 
compare the performance of adaptive time-stepping 
MATEXs (I-MATEX, S-MATEX), and TR (fixed time step $h=10ps$), adaptive time-steping TR. 
The transition time of 
the step input and the simulation time span is $10ps$ and $1\mu s$ 
for observing both high- and low-frequency behavior. 
For such large size application, in the respective of memory counsumption, 
we limit the Krylov subspace order $m$ within $[6,~30]$. During the simulation. 
We find such region can provide enough accuracy for matrix exponential approximation by I-MATEX and S-MATEX.
Table~\ref{tab:results_pdn_cases} shows the simulation runtime of MATEXs and TRs 
%where the fixed step size is set as $10ps$ 
to comply with the transition time. 
In the table, ``DC'', ``LU'' and ``Total'' 
indicate the runtime for DC analysis, LU 
factorization and the overall simulation, respectively. 
DC analysis is also via the LU factorization. 
For large cases, we can also adopt 
other techniques \cite{Xiong12}\cite{Yang12}\cite{Yu12} 
to improve the performance of DC analysis for both methods. 
``Spdp1'' shows the speedup of R-MATEX over fixed time step TR ($h=10ps$). 
``Spdp2'' shows the speedup of R-MATEX over TR (Adaptive). 
``Spdp3'' shows the speedup of R-MATEX over I-MATEX. 

\begin{table*}[ht]
    \centering
    \caption{Simulation runtime of PDN designs}
    \begin{tabular}{|c||r|| r|r|| r|| r| r || r|r|r| r| r|}
    \hline
	\multirow{2}{*}{Design} & \multirow{2}{*}{DC(s)} 
	&\multicolumn{2}{c||}{TR ($h=10ps$)} 
	&\multicolumn{1}{c||}{TR (Adaptive)} 
	& \multicolumn{2}{c||}{I-MATEX }  
	& \multicolumn{5}{c|}{R-MATEX ($\gamma=10^{-10}$) }  
	\\ 
	\cline{3-12} 
				 &	  &LU(s)   &Total   
				 & \multicolumn{1}{c||}{Total}	
				 &LU(s)	 &Total		  
				 &LU(s)	 &Total		& Spdp1 
				 & Spdp2
				 & Spdp3
				 \\ \hline
        D1			& 0.11	  & 0.12   & 392.46s	
				& 52.44s %44.80s         % TR A
				& 0.13    & 33.67s    % I-MATEX 
				& 0.13    & 23.76s	
				& 16.5X 
				& 2.2X 
				& 1.4X 
				\\ \hline
        D2 			& 2.72 
				& 3.41   & 2.80h	
				& 0.27h        
				& 3.78   & 0.21h   
				& 3.98	 & 0.19h
				& 14.7X
				& 1.4X 
				& 1.1X \\ \hline
        D3 			& 17.86   
				& 23.22 &  12.59h	
				& 1.42h %0.26h        
				& 24.88 & 1.06h   
				& 23.80	 &  0.68h	
				& 18.5X 
				& 2.1X
				& 1.5X 
				\\ \hline
        D4 			& 62.11 
				& 76.07  & 37.24h 
				& 4.27h 
				& 88.82  & 3.26h 
				& 77.96	 & 2.05h 
				& 18.2X
				& 2.1X 
				& 1.5X \\ 
    \hline
    \end{tabular}
    \label{tab:results_pdn_cases}
\end{table*}


\begin{figure}[ht]
    \centering
    \includegraphics[width=3.7in]{./figs/result_pdn_case1.eps}
    \caption{Result of D1}
    \label{fig:result_pdn_case1}
\end{figure}
\end{comment}
 

\begin{comment}
\begin{table*}[t]
 \centering
    \caption{TR with adaptive stepping (TR(adpt)) vs. I-MATEX vs. R-MATEX.
    ``Total(s)'' is the total runtime of test cases.
    ``DC(s)'' records the time to obtain initial condition.
    Spdp$^1$: Speedup of I-MATEX over TR(adpt); 
    Spdp$^2$: Speedup of R-MATEX over TR(adpt);  
    Spdp$^3$: Speedup of R-MATEX over I-MATEX.  
   }
    \begin{tabular}{|c|r|r||r|r|| r| r|  r| }
    \hline
    \multirow{2}{*}{Design} & 
    \multirow{2}{*}{DC(s)} &
    \multicolumn{1}{c||}{ TR(adpt)} & 
    \multicolumn{2}{c||}{I-MATEX} & 
    \multicolumn{3}{c|}{R-MATEX } 
    \\ 
    \cline{3-8} 
    &
			     &Total(s) % s I-M 
			     &Total(s) 
			     &SPDP$^1$
			     &Total(s) 
			     &SPDP$^2$
			     &SPDP$^3$
			     \\ \hline  
	ibmpg1t		     & 0.13     
			      & 29.48
			      & 27.23
			     & 1.3X
			     & 4.93	  
			     & 6.0X 
			     & 5.5X
			     \\ \hline 
    	ibmpg2t 	     & 0.80  
			    & 179.74 
			   &  124.44
			   & 1.4X 
			   & 25.90  
			   & 6.9X 
			   & 4.8X
			    \\ \hline 
    	ibmpg3t 	     & 13.83 
			    & 2792.96 
			   & 1246.69  %1988.98
			   & 2.2X
			   & 244.56  
			   & 11.4X
			   & 5.1X
			   \\ \hline
    	ibmpg4t 	    & 16.69
			    & 1773.14 
			    & 484.83 %625.12
			   & 3.7X 
			    & 140.42 %231.06 
			    & 12.6X
			    & 3.5X
			    \\ \hline
    	ibmpg5t 	    & 8.16 
			   & 2359.11
			    & 1821.60 
			    & 1.3X 
			    & 337.74 
			    & 7.0X
			    & 5.4X 
			    \\ \hline 
    	ibmpg6t 	     & 11.17
		    	& 3184.59
			    & 2784.46 %4059.30
			   & 1.1X
			    & 482.42  
			    &   6.6X 
			   &5.8X 
			    \\ 
    \hline
    \end{tabular}
    \label{tab:perf_mexp}
   % \begin{tablenotes}\scriptsize
%	\item [1] Speedup of I-MATEX (mn) over TR ($h=10ps$) 
%	\item [2] Speedup of R-MATEX (mn)over I-MATEX (mn) 
%	\item [3] Speedup of R-MATEX (mn) over TR ($h=10ps$)
%	\item [4] Speedup of R-MATEX (mn) over TR (adpt)      
%	\vspace{-0.5in}
%    \end{tablenotes}
    %\end{threeparttable}
\end{table*}
\end{comment}

%\subsection{Distributed MATEX Performance}

\subsection{Performance of I-MATEX and R-MATEX in Sec. \ref{sec:frame_motivation}}
\label{sec:resultI}
We compare our proposed I-MATEX and R-MATEX against the popular TR-FTS on the IBM power grid benchmarks \cite{Nassif08_Power}. Among the current sources, the smallest interval between two breakpoints is $h_{upper}=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~\ref{tab:spec_ibmpg} shows the details of each benchmark circuit of which size ranges from $54$K up to $3.2$M. 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. \ref{fig:interleave}). 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.  

\begin{table}[t]
\centering
\caption{Specifications of IBM power grid benchmarks.}
\begin{tabular}{|c|r|r|r|r|r|r|r|} 
\hline
Design    & \#R     & \#C      & \#L  & \#I    & \#V  & \#Nodes  \\ \hline
ibmpg1t   & 41K   & 11K    & 277  & 11K  & 14K  & 54K  \\ \hline
ibmpg2t   & 245K  & 37K    & 330  & 37K  & 330    & 165K  \\  \hline
ibmpg3t   & 1.6M & 201K   & 955  & 201K  & 955    & 1.0M \\ \hline
ibmpg4t   & 1.8M & 266K   & 962  & 266K  & 962    & 1.2M \\ \hline
ibmpg5t   & 1.6M & 473K   & 277  & 473K  & 539K   & 2.1M \\ \hline
	%092148 
ibmpg6t   & 2.4M & 761K   & 381  & 761K &836K  & 3.2M \\ \hline
\end{tabular}
\label{tab:spec_ibmpg}
\end{table} 

After DC analysis in TR-FTS, we LU factorize matrix once for the  later transient simulation, which only contains time stepping. % as \cite{Xiong12}. 
Actually, multiple factorized matrices can be deployed \cite{Li11, Ye08_ISQED}. 
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 \ref{tab:perf_mexp}, we record the total simulation runtime {\bf{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 {\bf{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 $\mbf 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 SPDP$^{r}_{tr}$. 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 \ref{tab:perf_mexp}, our single mode R-MATEX achieves the average speedup {\bf $5\times$} over TR-FTS. Note the average speedup number of single mode R-MATEX over TR-FTS for the original IBM benchmark (ibmpg1t$\sim$ibmpg6t) is less than the speedup of the new test cases (ibmpg1t\_new$\sim$ibmpg6t\_new). As we mentioned before, ibmpg1t\_new$\sim$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$\sim$ibmpg6t\_new than ibmpg1t$\sim$ibmpg6t. 
%, the later of which use larger step size $10ps$ in TR-FTS.}

\begin{table*}[t]
	\centering
	 	\caption{%\color{red}
		 {Performance comparisons (single computing node): TR-FTS, I-MATEX, and R-MATEX. 
		 {\bf DC(s)}: Runtime of DC analysis (seconds);
		 {{\bf $m_I$}: The maximum $m$ of Krylov subspace in I-MATEX.}
		 {\bf Tran(s)}: Runtime of transient simulation after DC (seconds), excluding the matrix factorization runtime; 
		 {\bf Total(s)}: Runtime of overall transient simulation (seconds); \
		  {{\bf Df(uV)}: Maximum and average voltage differences compared to provided solutions (uV); } 
		 { {\bf  $m_R$}: The maximum $m$ of Krylov subspace in R-MATEX}
		  %((seconds)  
		 {\bf SPDP$^{r}_{tr}$}: Speedup of R-MATEX over TR-FTS with respect to {\bf Tran(s)};
		  {\bf SPDP$^{r}_{i}$}:  Speedup of R-MATEX over I-MATEX with respect to {\bf Tran(s)}. 
		% {\bf SPDP$^{hr}_{tr}$}:  Speedup of HR-MATEX over TR-FTS.
	}
	}
	\begin{threeparttable}[b]
		\begin{tabular}{|c|r|r|r||r|r|r|r||r|r|r| r||r|r| 
				}
			\hline
			\multirow{2}{*}{Design} & 
			\multirow{2}{*}{DC(s)} &
			\multicolumn{2}{c||}{TR-FTS} & 
			\multicolumn{4}{c||}{I-MATEX} & 
			\multicolumn{4}{c||}{R-MATEX} & 
			\multicolumn{2}{c|}{Speedups} 
		%	}
			\\ 
			\cline{3-14} 
			&	
			& Tran(s)	&Total(s)	
			% &LU(s)
			& $m_I$	
			& Tran(s)
			& Total(s) 
			& Df(uV)
			& $m_R$	
			& Tran(s)
			& Total(s) 
			& Df(uV)
			% s I-M 
		%	&Tran(s)
		%	& Total(s)		
			& SPDP$^{r}_{tr}$ %\footnote{ over TR }%  (/TR.)   
			& SPDP$^{r}_{i}$ %\footnote{ over TR }% (/TR.)	
%			& SPDP$^{r}_{tr}$ %\footnote{ over TR (adpt)}%  (/TR.) 
			% m I-M	
		%	&LU(s)	&Total(s) % s R-M	
			%&Spdp2\footnote{ over TR }
			%	(/TR.)
%			&Total(s) % m R-M 
			%&SPDP\footnote{ over I-MATEX(mn) } %	(/I- (mn))
			\\ \hline  \hline 
			ibmpg1t		     
			& 0.2  & 5.7 	& 6.00 
			& 30  	& 28.8 	& 28.9	& $58 \backslash 9.8$  
			& 5 	& 10.1 & 10.3  & $45 \backslash 6.8 $ & $0.6\times$ 
			& $2.9\times$
			\\ \hline 
			ibmpg2t    
			& 0.8  & 40.0 & 41.9 	
			& 28 & 130.0 &  130.9	& $92  \backslash 10.5$
			& 5  & 35.6 & 37.4 & $45 \backslash 6.8$ 
			&  $1.1\times$	&  $ 3.7\times$ 
			%	&  $1.9\times$	 
			\\ \hline 
			ibmpg3t 	     
			&  16.4  &  263.2  	& 295.0  
			& 29 & 1102.5  & 1115.1 	& $95 \backslash 20.4$ 
			&  5 & 	275.5  & 301.0 	&  $95 \backslash 18.5$
			& 	$1.0\times$ 	&	$4.0 \times$		% & $ 1.6\times$  
			\\ \hline
			ibmpg4t    
			& 13.5 & 460.8 	& 501.9
			&  29 
			 & 433.8	& 458.2 &  $101  \backslash 39.3$
			&  5& 200.5  & 239.1   & $99 \backslash 34.2$ 
			 & $ 2.3\times$
			& 	 $ 2.2\times$ %	& $ 6.4\times$
			\\ \hline
			ibmpg5t   
			&  9.0 	&  476.6 	&  	498.0 
			&  30 &  1934.4	& 1944.5  	&  $29  \backslash 5.6$		
			&  5 & 383.1   & 401.9 &  $29  \backslash 4.4$    
			&  $ 1.2\times$
			& $ 5.0\times$	% &$ 2.1\times$
			\\ \hline 
			ibmpg6t    
			&  15.3 	&  716.0 	& 749.1  	
			& 25  &  2698.9	& 2713.7  	&   $39\backslash 8.6$
			& 5 & 773.5 & 800.5 	&   $33\backslash 5.6$	
%			& 323.20 & 357.10
	&  $0.9\times$ 
			& $3.5\times$ %	&  $2.2\times$
			\\ 
			\hline
		 % \hline 
			ibmpg1t\_new	   
			&  0.2	&  51.3 & 51.7
			& 30 &  27.2  & 27.4 &  $58  \backslash 9.8$	  
			& 5 & 11.7   & 12.1  & $53 \backslash 6.9 $	 
			& $4.4\times$ 	& $2.3\times$% & $13.4\times$	 
			\\ \hline 
			ibmpg2t\_new 	      
			&  0.9   &  431.4 	&  	433.5
			&  28 	& 114.9  & 115.7 	&  $49 \backslash 10.5$	   
			& 5 	& 43.3 	& 44.9  &  $33\backslash 5.6$	  
			& $10.0\times$ 		& $2.7\times$	 %	&  $20.5\times$
			\\ \hline 
			ibmpg3t\_new	        
			&  16.3 &  3716.5  	&  	3749.0
			&  29	&  1219.3 & 1232.6  	&  $95 \backslash 20.4$	
			&  5	& 481.7   & 508.2 & $95 \backslash 18.9$	    
			&  $7.7\times$ 
			&  $2.5\times$	 
			\\ \hline
			ibmpg4t\_new 	     
		 	& 18.3 	& 5044.6 & 5085.3 
		 	&  29   
		 	& 753.5 	& 776.4 	&  $101 \backslash 39.3$		 
			& 6	& 350.9 & 387.2 & $99 \backslash 34.2$		
			&  ${14.4}\times$ & 	 $2.1 \times$ 
			\\ \hline
			ibmpg5t\_new 	     
			& 10.5  & 5065.9 & 5110.1
			& 30  	& 2494.0	&  2504.7
			& $ 30 \backslash 5.6 $ 
			& 5  & 746.2	&  766.4 & $ 30 \backslash 4.4 $   &  $6.8\times$			
			&  $3.3\times$	  %	&  $20.8\times$
			\\ \hline 
			ibmpg6t\_new 	       
			&  13.1	&  7015.3 	&  	7059.7
			& 25 & 3647.9 & 3663.1 &  $39 \backslash 8.6 $ 
			& 6	& 895.1	&  923.1 &  $33\backslash 7.3 $ 
            &  $7.8\times$  &	$4.1\times$	 %	&  $19.4\times$
			\\ 
			\hline	
			\hline
			 Average   & --- & ---
			  & --- & --- & --- & --- 
			  & $65\backslash 15.7$ &  --- & ---
			 &  --- & $57\backslash 12.8$ 
		 & $5\times$
			 &$3\times$ % & $13.4\times$
			 \\
			 \hline 		
		\end{tabular} 
	\end{threeparttable}
	\label{tab:perf_mexp}
\end{table*}



\subsection{Performance of DR-MATEX in Sec. \ref{sec:dist_frame}}  
 
We test our distributed DR-MATEX in the following experiments with the same IBM power grid benchmarks. These cases have many input transitions ($GTS$) that limit step sizes of R-MATEX. 
We divide the region before the computation of simulation. 
%Then stepping can actually benefits from  the large time step. 
We decompose the input sources by the approach discussed in Sec. \ref{bump} and obtain much fewer transitions of $LTS$ for computing nodes.
%In IBM power grid benchmarks, 
The  original input source  numbers are over ten thousand in the benchmarks. 
However, based on ``bump'' feature (as shown in Fig. \ref{fig:ckt_ts_decomp}), we obtain a fairly small numbers for each computing node, which is shown as \emph{Group \#} in Table  \ref{tab:pg_rmatex}. 
(Now, the fact that hundred machines to process in parallel is quite normal \cite{mesos, cloud} in the industry.)
%For scalable experiment and simulations, we utilize the equivalent idea 
%of MATLABMPI implementations \cite{matlabmpi}. 
We pre-compute $GTS$ and $LTS$ groups and assign sub-tasks to corresponding nodes\footnote{ 
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.}.
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. \ref{fig:drmatex}) for the computation of their own transient simulation tasks, and have no communications with  others, our framework falls into the category of \emph{Embarrassingly Parallelism} model. We can easily emulate the multiple-node environment.  
%using our limited number of servers.  
% In this sub-section, we validate our proposed framework without considering too much implementation issues.   
%for this distributed computing. 
%We assign one MATLAB instance at each node with single thread computing mode by using ``matlab -singleCompThread''.  
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 \ref{tab:pg_rmatex}. 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$\sim$ibmpg6t; 
We use $h=1ps$ for ibmpg1t\_new$\sim$ibmpg6t\_new. 
Therefore, we have 10,000 pairs of forward and  backward substitutions for stepping.
%and does not benefit from input transition decomposition.
%Similar TR serves as a comparison baseline, e.g. \ref{wang13}. 
In DR-MATEX, the circuit solver uses R-MATEX with {\bf $\gamma = 10^{-10}$}, which is set to sit among the order of varied time steps during the simulation (since Sec. \ref{sec:compare_krylov} discusses the insensitivity of $\gamma$ 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.

% The experiment is carried on Linux servers with Intel(R) Xeon (R)~E5-2640 v3 2.60GHz processor and over 125GB memory.
In Table \ref{tab:pg_rmatex}, our distributed mode gains {{up to  $98\times$}}   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 \times$  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 \cite{nicslu}.  


\begin{table*}[t] 
	\caption{
		{ 
			{The performance of DR-MATEX (Distributed R-MATEX)}. {\bf Group \#}: Group number of the testcases. 
%			Divided by the first breakpoints. 
			This number represents the total number of simulation sub-tasks for the design;
			{\bf Tran(s)}: Runtime of transient simulation after DC (seconds); 
			{\bf Total(s)}: Runtime of overall transient simulation (seconds); 
			%vs. TR ($h=10ps$); 
			{\bf Max. Df.(V)}  and {\bf Avg. Df.(V)}: maximum and average differences compared to the solutions of all output nodes provided by IBM power grid benchmarks.  
			%by IBM Power Grid Benchmarks solutions  and our internal reference solutions for the modified cases;
		%	{\bf SPDP$_{i}$}: Speedup over R-MATEX in Table \ref{tab:perf_mexp}.
			{\bf SPDP$_{tr}$}: Speedup over TR-FTS's {\bf Tran(s)} in Table \ref{tab:perf_mexp};
			{\bf SPDP$_r$}: Speedup over R-MATEX's {\bf Tran(s)} in Table \ref{tab:perf_mexp};
			{\bf Peak $m$}: the peak dimension used in DR-MATEX  for MEVP; 
			{\bf Mem. Ratio over TR-FTS}: The peak memory comparison between the maximum memory consumption of DR-MATEX  over TR-FTS in Table \ref{tab:perf_mexp}.
			%Spdp$^4$:$ {t_{1000}/t_{rmatex}}$ speedups of transient stepping of MATEX over TR;
			%Spdp$^5$: $ {t_{t\_total}/t_{r\_total}}$ total runtime speedups of MATEX over TR. 
		}
	} 
	\centering
	\begin{tabular}{|c||r|r|r|r|r||r|r||c|c|}
		\hline
		\multirow{2}{*}{Design}&
		\multicolumn{5}{c||}{DR-MATEX}&
		\multicolumn{2}{c||}{Speedups}
		& Peak   &
		{Mem. Ratio} 
		\\
		\cline{2-8}
		& Group \#
		& Tran(s)
		& Total(s)
		& Max Df.(V)
		& Avg Df.(V)
		& SPDP$_{tr}$
	%	& SPDP$_{hr}$
		& SPDP$_{r}$
		& $m$
		& over TR-FTS
		%& Total(s)
		\\
		\hline \hline 
		ibmpg1t & 100  &1.4  &  1.9 
		&  5.3e-5   & 8.6e-6   
		& $4.0\times$ 
		&  $7.1\times$      & 6 & 1.9
		\\
		\hline
		ibmpg2t & 100 & 8.9
		& 11.4 & 
		4.6e-5   & 8.6e-6    
		&  $4.5\times$  
	%	&  $11.7\times$   
		&   $4.0\times$  
		& 7 & 1.9 
		\\
		\hline
		ibmpg3t & 100 & 91.7  & 129.9
		 %  69.71 & 124.41    
		&  9.6e-5 & 19.7e-6   &  $2.9\times$ 
	%	&  $9.5\times$  
		&   $4.4\times$  & 6  &  1.5 
		\\
		\hline
		ibmpg4t & 15  & 52.3 & 112.2	&  9.9e-5  & 27.9e-6    &  $8.8\times$   &    $3.8\times$  &  8 & 1.4  
		\\
		\hline   
		ibmpg5t & 100  	& 148.4 & 178.9  		&  9.0e-5  &  1.1e-6     
		&    $3.2\times$ 
		%&   $6.4\times$  
		&   $2.6\times$   &  7 % 1.4% 6.3/4.6
		& 1.5
		\\
		\hline
		ibmpg6t  & 100    
	& 189.9  & 234.2	&  3.4e-5  &  7.2e-6   
		&  $3.8\times$    
		& $4.1\times$   & 7 
		& 1.5 % 8.5/6.1~3.613.3
		\\
		\hline  
		ibmpg1t\_new  &  {\color{black} 100} 
		&  2.4   & 2.8 
		&  5.3e-5   & 8.6e-6     
		&  $21.8\times$  
		% & 3.8$\times$   
		& $5.0\times$  & 6  & 1.9
		\\
		\hline
		ibmpg2t\_new &{\color{black} 100}  
		% & 2.27 & 5.59
		& 5.6 & 7.0
		%	& 2.27 & 4.39   
		& 4.6e-5 & 8.6e-6  
		& $61.6\times$   
		% & 9.3$\times$  
		 & $6.2\times$ &  7 & 1.9
		\\
		\hline
		ibmpg3t\_new  &  {\color{black} 100} & 	103.0 & 140.9 
		%  69.71 & 124.41
		% 47.61 & 113.87 % 136.96 & 186.28
		&  9.8e-5 & 19.9e-6 
		&  $25.6\times$  
		% &  $ 6.4\times$  
		 &  $3.3\times$  
		%  &  $4.2\times$  &  $3.7\times$   &  $183.6\times$  
		& 7  & 1.5 
		\\
		\hline
		ibmpg4t\_new & {\color{black} 15}   
	& 51.5 & 108.4
	%		& 40.11 & 82.01  
		&  9.9e-5  & 27.6e-6 
		& {$98.0\times$}  
		 % & $2.5\times$ 
		  &  $6.8\times$  & 8   & 1.4 %  7.5/5.7 
		\\
		\hline
		%ibmpg5t & 331.67  & 350.33  &  100 & 35.67 & 54.21 &  0.7E-4 &  1.1E-5 &  9.3X  & 6.5X  
		ibmpg5t\_new  &  {\color{black} 100}  
		& 185.6 & 227.8
%			& 45.56  & 69.20 
	 & 9.9e-5  &  2.2e-6     
		& $27.3\times$ 
		% & $5.3\times$ 
		 & $4.0\times$ 
		& 7 &  1.5
		\\
		\hline
		ibmpg6t\_new  & {\color{black} 100}   
	&  274.8  
&  317.7
%		& 62.92 	& 98.57
		&  3.4e-5  &  7.1e-6
         &  $25.5\times$  %&  $5.7\times$
		 &  $3.3\times$   &  7  & 1.5   
		\\
		\hline 	\hline
		 {Average}
		& --- & --- & --- &  {7.1e-5} &  {12.3e-6}   
		& { $26\times$}
		 % & $6.4\times$ 
		 &  { $5 \times$} &  {6.7} 
		& 
		 {1.6}
		\\
		\hline 
	\end{tabular}
	\label{tab:pg_rmatex}
\end{table*}



%computing runtimes for the transient part exculding LU.
%For R-MATEX, we use the maximum values among this 
%to define the $t_{adpt}$ in Table. \ref{tab:pg_rmatex}.
%MATEX with grouping is parallelizable by simulating each group in 
%multicore or cluster system.  
%The average number of pairs of backward and forward substitutions for Krylov subspace generation is around 60 ($km$ in Eq. (\ref{eqn:tr_spd})),
%while TR($h=10ps$) has 1000 pairs ($N$ in Eq. \ref{eqn:tr_spd}) on each cases.
%The reduction of these substitutions brings large speedups in the transient computing (the part of $K(T_H+T_e)$ in Eq. (\ref{eqn:tr_spd}) also costs portions of computational efforts). 
%and drags the number of speedups. 
%With the huge reduction of simulation time on this part, the serial LU and DC parts in MATEX  expose dominating roles, 



 

\begin{comment}
\begin{figure}[ht]
    \centering
    \includegraphics[width=3.4in]{./figs/result_ibmpg6t.eps}
    \caption{Result of IMBPG6T}
    \label{fig:result_ibmpg6t}
\end{figure}
\end{comment}


%which can be further improved using advance methods.
%~\cite{nicslu},
%which is another important study area.
%There are other way to decompose the input transitions, we use one simple
%strategy based on ``hump'' shapes, and achieve such large speedups.  

%In order to demonstrate the idea, we simply device a naive grouping scheme, which will take identical input sources together as a group. 
\begin{comment}
We assign each group to one computing nodes. 
In the experiment, we notice the message passing between 
our workstations are extreme small.
Time anchors are generated as a wall-lock, 
that visible to each node, so there are not other
communication overhead.
For example, for the message receiving only takes 0.0007 seconds, which is negligible. 
Table~\ref{tab:perf_mexp_group} shows the performance of MATEX with grouping. 
In this table, ``Groups'' is the number of groups for each 
designs, ``Min. Total`` and ``Max. Total'' is the maximum runtime among different nodes, 
%and ``Ideal Spdp'' records the ideal speedup with parallelization. 
Design ibmpg4t has only 5 groups because the alignment is better. 
The performance of MEXP has a significant improvement after grouping. 
The average ideal speedup over TR is 6 times. 
MATEX shows another parallelism in the power grid simulation from the perspective of 
input sources. 
\end{comment}

\begin{comment}
\begin{table*}[ht]
    \centering
    \caption{Performance of R-MATEX w/ Grouping}
    \begin{tabular}{|c||r |r|r|r|r|r|}
    \hline
    \multirow{2}{*}{Design} & \multirow{2}{*}{TR Trans(s)} & \multicolumn{5}{c|} {R-MATEX w/ grouping}		      \\ \cline{3-7} 
			    &			     & Groups	    & Max. Trans (s)	& Min. Trans (s)  &Spdp1  &Spdp2\\ 
    \hline ibmpg1t 	    & 5.94		     & 25	   & 0.84		& 0.48		  & 7.1X  & 12.4X \\
    \hline ibmpg2t 	    & 26.98                   & 25	   & 3.93		& 2.38		  & 6.7X  & 11.3X \\
    \hline ibmpg3t 	    & 245.92                  & 25	   & 38.09     	&21.96		  & 6.5X & 11.2X \\
    \hline ibmpg4t 	    & 329.36                  & 5	   & 33.55		&30.12		  & 9.8X  & 10.9X \\
    \hline ibmpg5t 	    & 331.67                  & 25	   & 58.54 		& 33.59		  & 5.7X & 9.9X \\
    \hline ibmpg6t 	    & 542.04                  & 25	   & 81.40 		&51.68		  & 6.7X & 28.0X \\
    \hline
    \end{tabular}
    \label{tab:perf_mexp_group}
\end{table*}


\subsection{Complexity}

The performance of MATEX and TR are dominated by the number of backward/forward substitutions. This is because the time span 
is relatively long in a power grid, and therefore, the expensive LU fatorization is amortized during the simulation. In MATEX, 
the number of backward/forward substitutions at each time step is actually $m$, which is the number of vectors forming rational 
Krylov subspace. Since $m$ is small, e.g., $5\sim 6$ in IBM benchmark suite, and MATEX uses large step size than TR, MATEX has 
less numbers of substitution operations in general. Table~\ref{tab:subs} presents number of backward/forward substitions for 
each method. ``Sub.'' is the total number of substituion operations along the simulation. MATEX has less operations but close to
TR because the ill alignment of input sources leads to many tiny step sizes whereas MATEX with grouping has much less operations 
due to well aligned inputs. These results also validate Tables~\ref{tab:perf_mexp} and \ref{tab:perf_mexp_group} where MATEX would 
has less speedup due to ill alignment issue, and the grouping scheme could improve the performance. 

\begin{table}[h]
\centering
    \caption{Number of backward/forward substitutions}
    \begin{tabular} {|c|c|c|c|}
    \hline
    \multirow{2}{*}{Design} &TR	 &MATEX		  &MATEX w/ grouping \\ \cline{2-4} 
			    & \# Sub.	 & \# Sub.        &  \# Max Sub.   \\ \hline
    ibmpg1t   		    & 1000	 &  830		  &  120	    \\ \hline 
    ibmpg2t   		    & 1000  	 &  830   	  &  120 	    \\ \hline
    ibmpg3t   		    & 1000  	 &  874   	  &  120 	    \\ \hline 
    ibmpg4t   		    & 1000  	 &  349   	  &  143 	    \\ \hline  
    ibmpg5t   		    & 1000  	 &  982   	  &  139 	    \\ \hline
    ibmpg6t   		    & 1000  	 &  932   	  &  134 	    \\ 
    \hline 
    \end{tabular}
    \label{tab:subs}
\end{table}
\end{comment}

\begin{comment}

\subsection{Scalable IBM Power Grid Suite with Input Grouping}

Following the previous subsection, we not the input is the limitation for 
adaptive steping in our MATEX framework. 
We utlize input grouping~\ref{sec:group}.
Ideal one input sources per group.

\subsubsection{2 groups}

\subsubsection{4 groups}

\subsubsection{8 groups}

\subsubsection{16 groups}

\subsubsection{32 groups}

\subsubsection{64 groups}




Secondly, we consider a single input to run ibmpg1t-ibmpg6t 
to further illustrate our method. We see there are huge speedups than using trapezoidal method.
\end{comment}
%\subsubsection{Full Solver}
%This is the ideal speedups when computing resources are available.

%\subsubsection{Model Order Reduction Preprocessing}
%In this part, we compare the MOR method and put them into trapezoidal method
%and our MATEX framework, in order to compare our method accurate.



\begin{comment}
\begin{table}[h]
    \caption{Comparisons among MEXP, I-MATEX and R-MATEX with RC cases}
    \centering
    \begin{tabular}{|c||c|c|c|c|}
    \hline
	Method  &   $m_{p} $  & Err(\%)
%	&  Time(s)
	& Speedup 
	& Stiffness\\
	%& {$\frac{Re{\lambda_{smallest}(\mbf A)}}{Re({\lambda_{largest}(A)}}$}  \\ 
	\hline
	MEXP	  	&   229  	
	&  0.510
%	&  601.87  
	& -- &
	\multirow{3}{*}{$2.1 \times 10^{16}$ } 
	\\
	\cline{1-4}
	I-MATEX	 	&14       &  0.004
%	& 0.23
	& 2616X   &   \\
	\cline{1-4}
	R-MATEX	   & 12   &  0.004  
%	& 0.22   
	&  2735X   & \\ 
	\hline
	\hline
	MEXP	  	&   224  	&  0.004  
	%&  134.32
	& --  
	& 	    \multirow{3}{*}{$2.1 \times 10^{12}$ } \\
	\cline{1-4}
	I-MATEX		& 14       &0.004   
	%&   0.23
	&  583X  &   \\
	\cline{1-4}
	R-MATEX	&  12   &  0.004  
%	& 0.20 
	&  611X  & \\ 
	\hline
	\hline
	MEXP	   	&   223  & 0.004  
%	& 52.87  
	&  -- & 	    \multirow{3}{*}{$2.1 \times 10^{8}$ } 
	\\
	\cline{1-4}
	I-MATEX  	& 14      & 0.004   
%	&  0.23 
	& 229X  &   \\
	\cline{1-4}
	R-MATEX	  & 12   &  0.004 
%	& 		0.21  
	&  252X  & \\ 
	%\cline{1-6}
	\hline
    \end{tabular}
    \label{tab:rc_cases}
\end{table}
\end{comment}


\begin{comment}
\begin{table*}[htp]
%\scriptsize
    \centering
    \caption{Performance comparison TR, I-MATEX and R-MATEX of single-node computing (sn)/multi-node simulation (mn) versions,
    Spdp$^1$: Speedup of I-MATEX (mn) over TR ($h=10ps$); 
    Spdp$^2$: Speedup of R-MATEX (mn)over I-MATEX (mn); 
    Spdp$^3$: Speedup of R-MATEX (mn) over TR ($h=10ps$);
    Spdp$^4$: Speedup of R-MATEX (mn) over TR (adpt).  
   }
    \begin{threeparttable}[b]
    \begin{tabular}{|c|r|r|r||r||r|r || r| r|| r| r ||  r|r|r|r| }
    \hline
    \multirow{2}{*}{Design} & 
    \multirow{2}{*}{DC(s)} &
    \multicolumn{2}{c||}{TR ($h=10ps$)} & 
    \multicolumn{1}{c||}{TR (adpt)} & 
    \multicolumn{2}{c||}{I-MATEX(sn)} & 
    \multicolumn{2}{c||}{I-MATEX(mn)} & 
    \multicolumn{2}{c||}{R-MATEX(sn) } & 
    \multicolumn{4}{c|}{R-MATEX(mn)} 
    \\ 
    \cline{3-15} 
			     &		&LU(s)	&Total(s)	
			    % &LU(s)	
			     &Total(s) % s I-M 
			     &LU(s)
			     &Total(s) 
			     &Total(s) 
			     &SPDP\footnote{ over TR }% (/TR.)	
			     % m I-M	
			     &LU(s)	&Total(s) % s R-M	
			     %&Spdp2\footnote{ over TR }
			     %	(/TR.)
			     &Total(s) % m R-M 
			     &SPDP\footnote{ over I-MATEX(mn) } %	(/I- (mn))
			     & SPDP\footnote{ over TR }%  (/TR.)   
			     & SPDP\footnote{ over TR (adpt)}%  (/TR.)   
			     \\ \hline  
	ibmpg1t		     & 0.13     
			      &0.13   &  6.20 	
			      & 29.48
			     & 0.12  & 27.23
			    % & 0.12  
			     & 4.68 
			     & 1.3X
			     & 0.13 & 4.93	  
		%	     & 1.3X
			     %& 0.13 
			     & 1.10
			     & 4.5X 
			     & 5.6X
			     & 26.8X
			     \\ \hline 
    	ibmpg2t 	     & 0.80  
			    & 0.83   & 28.61	
			    & 179.74 
			   & 0.76   &  124.44
			  % & 0.76 
			   & 20.19 
			   & 1.4X
			   & 0.85 
			   & 25.90  
		%	   & 1.1X
			   %& 0.85 
			   & 5.56
			    & 3.6X 
			    & 5.1X
			    & 32.3X
			    \\ \hline 
    	ibmpg3t 	     & 13.83   & 12.72  & 272.47
			    & 2792.96 
			   & 12.87 & 1246.69  %1988.98
			  % & 12.58 
			   & 213.33 %280.66 
			   & 1.3X 
			   & 12.73 
			   & 244.56  
			 %  & 1.1X
			  % & 12.81 
			   & 63.78 
			   & 3.3X 
			   & 4.3X 
			   & 43.8X
			   \\ \hline
    	ibmpg4t 	    & 16.69    & 22.50   &	368.55
			    & 1773.14 
			    & 21.09 & 484.83 %625.12
			   % & 20.99 
			    & 148.45 %231.06 
			    & 2.5X
			    & 22.09  & 140.42    
			 %   & 2.6X
			   % & 22.02  
			    & 76.55  
			    & 1.9X  
			    & 4.8X
			    & 23.1X
			    \\ \hline
    	ibmpg5t 	    & 8.16    & 10.50  & 350.33  	
			   & 2359.11
			    & 9.92 & 1821.60 
			   % & 10.03 
			    & 304.72 %413.77
			    & 1.1X
			    & 10.77 &  337.74 
			 %   & 1.0X
			   % & 10.51  
			    &  78.99 
			    &    3.9X 
			    & 4.4X
			    & 29.9X
			    \\ \hline 
    	ibmpg6t 	     & 11.17    &14.17  & 567.38
		    	& 3184.59
			    & 13.19 & 2784.46 %4059.30
			   % & 13.55 
			    & 479.79 %856.09
			   & 1.2X
			    & 14.18 & 482.42  
		%	    & 1.2X
			   % & 14.15 
			    & 108.92 
			   & 4.4X 
			    &   5.2X 
			   & 29.2X 
			    \\ 
    \hline
    \end{tabular}
   % \begin{tablenotes}\scriptsize
%	\item [1] Speedup of I-MATEX (mn) over TR ($h=10ps$) 
%	\item [2] Speedup of R-MATEX (mn)over I-MATEX (mn) 
%	\item [3] Speedup of R-MATEX (mn) over TR ($h=10ps$)
%	\item [4] Speedup of R-MATEX (mn) over TR (adpt)      
%	\vspace{-0.5in}
%    \end{tablenotes}
    \end{threeparttable}
    \label{tab:perf_mexp}
\end{table*}
\end{comment}


%Note that even though previous works\cite{Chen01_Efficient}\cite{Su03} show that using iterative approach in TR 
%could also achieve adaptive step control, long simulation time span in power grid designs make direct 
%method with fixed step size more desirable\cite{Xiong12}\cite{Yang12}\cite{Zhao02}. 
%Our benchmark includes different industrial power distribution networks (PDNs), which contain the model of 
%board and package, and also includes IBM power grid suite~\cite{Nassif08_Power}. 
%The details and simulation results of 
%both industrial designs and IBM power grid suite are presented in the following subsections. 

\begin{comment}
\begin{table*}[tp]
	\centering
	\caption{Hybrid MATEX vs. TR ($h=10ps$); 
		Max. and Avg. Err.: maximum and average differences compared
		to all output nodes' 
		solutions provided by IBM Power Grid Benchmarks;
		Spdp$^4$:$ {t_{1000}/t_{rmatex}}$ speedups of transient stepping 
		of MATEX over TR;
		Spdp$^5$: $ {t_{t\_total}/t_{r\_total}}$ total runtime speedups 
		of MATEX over TR. 
	}
	\begin{tabular}{|c||r|r||r|r| r| r|r|r|r|}
		\hline
		\multirow{2}{*}{Design}&
		\multicolumn{2}{c||}{TR}&
		\multicolumn{7}{c|}{Hybrid MATEX}
		\\
		\cline{2-10}
		& $t_{1000}(s)$
		& $t_{t\_total}(s)$
		& Group \#
		& $t_{rmatex}(s)$ 
		& $t_{r\_total}(s)$
		& Max. Err.
		& Avg. Err.
		& Spdp$^4$
		& Spdp$^5$
		\\
		\hline
		ibmpg1t & 5.94   &   6.20  &  {\color{black} 100}  &  0.50  &  0.85  &  1.4E-4 & 2.5E-5 & 11.9X  & 7.3X  
		\\
		\hline
		ibmpg2t & 26.98  &   28.61  &  100  &  2.02  &  3.72  &  1.9E-4 &  4.3E-5 & 13.4X  & 7.7X  
		\\
		\hline
		ibmpg3t & 245.92  & 272.47  &  100  &  20.15  &  45.77 &  2.0E-4 &  3.7E-5 & 12.2X  & 6.0X  
		\\
		\hline
		ibmpg4t & 329.36  & 368.55  &  {\color{black} 15} & 22.35 &  65.66 &  1.1E-4 &  3.9E-5 & 14.7X  & 5.6X  
		\\
		\hline
		%ibmpg5t & 331.67  & 350.33  &  100 & 35.67 & 54.21 &  0.7E-4 &  1.1E-5 &  9.3X  & 6.5X  
		ibmpg5t & 408.78  & 428.43  &  100 & 35.67 & 54.21 &  0.7E-4 &  1.1E-5 &  11.5X  &  7.9X  
		\\
		\hline
		ibmpg6t & 542.04  & 567.38  &  100 & 47.27 & 74.94 &  1.0E-4 &  3.4E-5 &  11.5X  & 7.6X  
		\\
		\hline
	\end{tabular}
	%	\label{tab:pg_rmatex}
\end{table*}


\end{comment}
framework.tex
\subsection{MATEX Circuit Solver} 
\label{sec:frame_motivation}
 
%Matrix exponential based time integration kernel with Krylov subspace method can solve Eq. (\ref{eqn:dae}) with larger time step and high accuracy than lower order approximation methods such as BE, TR. 
We incorporate matrix exponential based integration scheme with Krylov subspace method into our MATEX framework, which is summarized in Algorithm \ref{algo:rmatex}. 
We set $\mbf X_1$ and $\mbf X_2$ in Line \ref{algo:rmatex:init} based on the choice of Krylov subspace method as follows, 
\begin{itemize}
%\item MEXP: $\mbf X_1=\mbf C, ~\mbf X_2 = \mbf G$.
\item I-MATEX: $\mbf X_1 =  \mbf G, ~\mbf X_2=\mbf C$
\item R-MATEX: $\mbf X_1 =  \mbf C+ \gamma \mbf G, ~\mbf X_2=\mbf C$
\end{itemize} 

For linear system of PDN, the matrix factorization in line \ref{algo:rmatex:lu} is only performed once, and the matrices $\mbf L$ and $\mbf U$ are reused in the while loop from line \ref{algo:rmatex:while} to line \ref{algo:rmatex:whileend}. 
%It is the same as traditional fixed step low order schemes, BE-FTS or TR-FTS. 
Line \ref{algo:rmatex:algo1} uses Arnoldi process  with corresponding inputs to construct Krylov subspace as shown in Algorithm \ref{algo:arnoldi}.
%\ \newline 

\begin{algorithm}[t]
 	\label{algo:rmatex}
 	%\footnotesize
 	{
 		%\caption{Arnoldi process in Matrix Exponential-Based Circuit Simulation with Krylov subspace Method}
 		\caption{MATEX Circuit Solver}
 		%$\rm K_m(\mbf X_1^{-1}\mbf X_2, \mbf v)$} %~\cite{Weng12_TCAD, Weng12_ICCAD}}
 		\KwIn{ $ \mbf C, \mbf G , \mbf B, \mbf u, \epsilon$,
 			%error tolerance $E_{tol}$, 
 			and time span $T$. 
 			%with a pre-defined constant $C_{ost}$ for the selection of time stepping time-integration kernel. 
 			%A value is large than 1 and smaller than average $m$ is a reasonable choice. 
 			% $L_{imit}$ is the threshold to decide whether MexpH uses MOR or not.
 			%, it is a cost for evaluation of TR-FTS and MEXP
 		}
 		%\KwOut{ $\mbf V_m, \mbf H_m,  m$ } 
 		\KwOut{  The set of $\mbf x$ from $[0, T]$.}
 		{Set $\mbf X_1, \mbf X_2$\;
 	\label{algo:rmatex:init}
	%{\tcc*[l]{R-MATEX: $\mbf X_1= \mbf C+ \gamma \mbf G$, $\mbf X_2=\mbf C$; or I-MATEX: $\mbf X_1 =  \mbf G$, $\mbf X_2=\mbf C$. }}
		%$\mbf X_1=\mbf C, \mbf X_2 = \mbf G$ for MEXP.
		 			%	[$\mbf L, \mbf U$]=LU($\mbf X_1$)\\
 			%$s=T_{start}$\;
 			$t=0$\;
 			$\mbf x(t) =$DC\_analysis\;  
 			$[\mbf L, \mbf U] = \text{LU\_Decompose}(\mbf X_1)$\; 	\label{algo:rmatex:lu}
 			\While { $t < T$}
 			{\label{algo:rmatex:while}
 				Compute maximum allowed step size $h$\;
 				Update $ \mbf P (t,h),\mbf F (t,h)$\;
 				Obtain  $\mbf x(t+h)$ by  {\bf Algorithm  \ref{algo:arnoldi}} with inputs $[\mbf L, \mbf U, \mbf X_2, h ,t, \mbf x(t), \epsilon, \mbf P (t,h),\mbf F (t,h)]$\; \label{algo:rmatex:algo1}
 				$t=t+h$\;
 			}\label{algo:rmatex:whileend}
 		}
 	}
\end{algorithm}

\subsection{DR-MATEX (Distributed R-MATEX Framework) by Decomposition of Input Sources, Linear Superposition, and Parallel Computation Model}
%Distributed/Parallel Processing.}
\label{sec:dist_frame}
\subsubsection{Motivation} \

%With the ability of large time stepping, the input term $\mbf b$ in Eq. (\ref{eqn:discrete_sol}) actually serves a double-edged sword in MEVP computation. 
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 $t_s$ 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.\footnote{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.}
% efficient stepping with variable time steps as well as the theoretical limitation of low order accuracy 
\vspace{+0.07in}

\subsubsection{Treatment and Methodology} \ 
\label{ltsgts}

% In order to solve the frequent breakpoints problem by distributed computing, 
In matrix exponential based integration framework, %as long as the approximation of $e^{h\mbf A}\mbf v$ is accurate enough, 
we can choose any time spot  $t+h \in [t,t_s]$ with computed Krylov subspace basis. The solution of $\mbf x(t+h)$ is computed by scaling the existing Hessenberg matrix $\mbf H$ with the time step $h$ as below
\begin{eqnarray}
\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V}_m e^{h \mbf H} \mbf e_1-\mbf P(t,h).
\label{sq:snapshot}
\end{eqnarray}
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. 
%, we can still use the last Krylov subspaces. 
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:
\vspace{+0.08in}
\begin{itemize}
	\item 
	%\textbf{Definition}: 
	\emph{Local Transition Spot} ($LTS$): the set of $TS$ at an input source to the PDN.
	\item
	%  \textbf{Definition}: 
	\emph{{Global Transition Spot}} ($GTS$): the union of $LTS$ among all the input sources to the PDN. 
	\item  
	%\textbf{Definition}: 
	\emph{{Snapshot}}:  a set $GTS \setminus LTS$ at {one input source}. 
\end{itemize} 
\vspace{+0.08in}

If we simulate the PDN with respect to all the input sources, the points in the set of $GTS$ are the places where generations of Krylov subspace cannot be avoided. For example, there are three input sources in a PDN (Fig. \ref{fig:rc_ckt}). The input waveforms are shown in Fig. \ref{fig:input}. The first line is $GTS$, which is contributed by the union of $LTS$ 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 $LTS$ and keeps track of $Snapshot$ for the later usage of summation via linear superposition. 
Between two LTS points $t$ and $t+h$,  the \emph{Snapshot} points $$t+h_1 < t+h_2 < \cdots < t+h_l  \in (t, t+h]$$ can reuse the Krylov subspace generated at $t$. 
%, which is mentioned in Sec. \ref{sec:discuss}. 
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.
% when processing these subtasks in parallel. 
Besides, when subtasks are assigned, there is no communication among the computing nodes, which leads to so-called \emph{Embarrassingly Parallel} computation model.  
\vspace{+0.07in}



\begin{figure}[t] 
	\centering
	\includegraphics[width=2.3in]{./figs/rc_ckt.eps}
	\caption{Part of a PDN model with input sources from Fig. \ref{fig:input}.}
	\label{fig:rc_ckt}
\end{figure}

\begin{figure}[ht] 
	\centering
	\includegraphics[width=3.4in]%{./figs/time_spots.eps}
	{./figs/time_spots_2.eps} 
	\caption{ 
		{Illustration of input transitions. 
			$GTS$: Global Transition Spots; $LTS$: Local Transition Spots;
			{$Snapshots$: the crossing positions by dash lines and LTS \#$k$ without solid points.} 
		}
	}
	\label{fig:input}
\end{figure} 

\begin{figure}[t] 
	\centering
	\includegraphics[width=3.4in]%{./figs/ckt_ts_decomp.eps}
	{./figs/source_large.eps}
	\caption{Grouping of ``Bump'' shape transitions for sub-task simulation.
		The matrix exponential based method can  utilize  adaptive stepping in each $LTS$ and reuse the Krylov subspace basis generated at the latest point in $LTS$. 
		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:
		$t_{delay}$: delay time;
		$t_{rise}$: rise time;
		$t_{width}$: pulse width; 
		$t_{fall}$: fall time; and
		$t_{period}$: period).
	}
	\label{fig:ckt_ts_decomp}
\end{figure}



\begin{figure*}[t]
	\centering
	\includegraphics[width=0.6\textwidth]%{./figs/hierachical_ckt_sim.eps}
	{./figs/hierachical_ckt_sim_new.eps}
	\caption{DR-MATEX: The distributed MATEX framework using R-MATEX circuit solver.}
	\label{fig:drmatex}
\end{figure*} 

\subsubsection{More Aggressive Tasks Decomposition} \ 
\label{bump}
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.\footnote{IBM power grid benchmarks provide the pulse input model in SPICE format.}  We group the input sources, which have the same $$(t_{delay}, t_{rise}, t_{fall}, t_{width})$$ into one set. 
For example, the input source \#1 of Fig. \ref{fig:input} is divided to  \#1.1 and \#1.2 in Fig. \ref{fig:ckt_ts_decomp}. 
The input source \#2 in Fig. \ref{fig:input} is divided to  \#2.1 and \#2.2 in Fig. \ref{fig:ckt_ts_decomp}.  Therefore, there are four groups in Fig. \ref{fig:ckt_ts_decomp}, Group 1 contains $LTS\#1.1$. Group 2 contains $LTS\#2.1$. Group 3 contains $LTS \#2.2$. 
Group 4 contains $LTS \#1.2$ and $\#3$. Our proposed framework MATEX is shown in Fig. \ref{fig:drmatex}. After pre-computing $GTS$ and decomposing $LTS$ based on ``bump'' shape (Fig. \ref{fig:ckt_ts_decomp}), we group them and form $LTS$ $\#1 \sim \#K$.\footnote{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.} 
\vspace{+0.07in}

\subsubsection{MATEX Scheduler in DR-MATEX} \

In DR-MATEX, the role of MATEX scheduler is just to send out $GTS$ and $LTS$ 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 $GTS$, $LTS$\#$k$, $\mbf P_k$ and $\mbf F_k$, which contain the corresponding $\mbf 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. \ref{fig:drmatex}), by when all nodes complete their transient simulations. 


Within each slave node, the circuit solver (Algorithm \ref{algo:ckt_solver}) 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. 


\begin{algorithm}
	\label{algo:ckt_solver} 
	{ 
		\caption{DR-MATEX: The distributed MATEX framework using R-MATEX at Node\#$k$.} 
		\KwIn{ $LTS$\#$k$, $GTS$, $\mbf P_k$, $\mbf F_k$, error tolerance $E_{tol}$, and simulation time span $T$.}
		%\KwOut{ $\mbf V_m, \mbf H_m,  m$ }
		\KwOut{Local solution $\mbf x$ along $GTS$ in node 
			$k \in [1,\cdots,S]$, where $S$ is the number of nodes}
		{   
			%	[$\mbf L, \mbf U$]=LU($\mbf X_1$)\\
			%	Normalize the intial vector 
			$t=0$, $\mbf X_1= \mbf C + \gamma \mbf G$, and $\mbf X_2 = \mbf C$\;
			%$s=T_{start}$\;
			$\mbf x(t) = \text{Local\_Initial\_Solution}$\;
			%$\mbf x(t) = local\_DC\_analysis$; 
			%	{\tcp*{ DC for input at \#$s$}} 
			%$[\mbf L,\mbf U] = \text{LU\_Decompose}( \mbf C+\gamma \mbf G)$\;
			$[\mbf L,\mbf U] = \text{LU\_Decompose}( \mbf X_1)$\;
			\While { $t < T$}
			{
				Compute maximum allowed step size $h$ based on $GTS$\;
				% and	$ \alpha = \frac{2h}{\gamma}$\; 
				{
					%\tcp{\color{blue} R-MATEX mode: acurate large step marching}
					\If{$t \in $  LTS\#$k$}
					{
						\tcc{\small  Generate Krylov subspace for the point at $LTS\#k$ and  compute $\mbf x(t+h)$}
						%	\tcp*{Reach LTS and update $\mbf V_m, \mbf H_m, \mbf v, s$}
						$[\mbf x(t+h), \mbf V_m, \mbf H_m,\mbf v] = \text{MATEX\_Arnoldi}(\mbf L,\mbf U,\mbf X_2,$ $ h, t, \mbf x(t), \epsilon, \mbf P_k(t,h), \mbf F_k(t,h))$\;
						%\label{algo2:arnoldi}
						%	Update $\mbf x(t +h) $ based on Eq. (\ref{eqn:new_exact})\; 
						$a_{lts} = t$\; 
						%\tcp{put an LTS anchor for current $t$} 
					}
					\Else{
						\tcc{\small Obtain $\mbf x(t+h)$ at $Snapshot$ with computed Krylov subspace}
						%	    \For  {$s_i \in (t_i + h,t_{i+1}]$}
						{
							%// Take snapshots by reusing $\mbf V_m , \mbf H_m, \mbf v$ \;
							$h_{a} = t + h - a_{lts}$\;
							$\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V}_m e^{h_{a} \mbf H_m} \mbf e_1 - \mbf P_k(t,h)$\;
	     				}
					}
					$t= t+h$\;
				}
			}    }
		}
	\end{algorithm}
	
\subsection{Runtime Analysis of MATEX PDN Solver } 
\label{sec:complexity} 

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 $T_{bs}$. 
The total time of serial parts is $T_{serial}$, which includes matrix factorizations, result collection, etc. 
For $\mbf x(t+h)$, the evaluation of matrix exponential with $e^{h\mbf H_m}$ is $T_{H}$, which is in proportion to the time complexity $O(m^3)$. Besides, we need extra $T_e$ to form $  \mbf x(t+h)$, which is proportional to $O(n m^2)$ by $\beta\mbf V_m e^{h\mbf H_m} \mbf e_1$. 

Given $K$ points of $GTS$, without decomposition of input sources, the runtime is \begin{eqnarray} 
\label{eq:full}
KmT_{bs} + K (T_{H}+ T_e)  + T_{serial}.
\end{eqnarray}
After dividing the input transitions and sending to enough computing nodes, we have $k$ points of $LTS$ for each node based on feature extraction and grouping (e.g., $k=4$ for one ``bump'' shape feature). The total computation runtime is  %in the order of 
\begin{eqnarray} 
\label{eq:reduced}
kmT_{bs} + K (T_{H}+ T_e)  + T_{serial},
\end{eqnarray}
where  $K (T_{H}+ T_e)$ contains the portion of computing $Snapshot$ in DR-MATEX mode. 
% Therefore, the major runtime gain is from the $K$ of Eq. (\ref{eq:full}) to $k$ of Eq. (\ref{eq:reduced}).
The speedup of DR-MATEX over single MATEX is  
\begin{eqnarray}
\text{Speedup}=\frac{ KmT_{bs} + K(T_{H}+ T_e)+ T_{serial}} { kmT_{bs}+ K(T_{H}+ T_e)+T_{serial}}.
	%\text{Speedup}=\frac{ KmT + K(m^3+ T_e)+ T_{s} }  { kmT+ K(m^3+T_e)+T_{s}}
\end{eqnarray}

For R-MATEX, we have  small $m$. Besides, $T_{bs}$ is relatively larger than $(T_{H}+T_e)$ in our targeted problem. 
Therefore, the most dominating part is the $KmT_{bs}$ in Eq. (\ref{eq:full}). 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 $$NT_{bs}+T_{serial}.$$ 
Then, the speedup of distributed DR-MATEX over the traditional method is 
\begin{eqnarray}
	\label{eqn:tr_spd}
	\text{Speedup}'=\frac{NT_{bs}+T_{serial} }  { kmT_{bs}+ K(T_{H}+T_e)+T_{serial}}.
\end{eqnarray}
%Usually, $N$ is much larger $km$. 

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. 
%When elongating time span of simulation, $N$ will increase, while $k$ will not change too much due to its irrelevance to time span (may bring more input transition features and increase computing nodes), 
% Therefore, our MATEX is more robust and has larger speedups. 
 
intro_ljw.tex
 
\begin{comment}
Circuit simulation is an important and 
{\color{red} heavily used process}
during the cycle of integrated circuit (IC) designs.
The increasing integration between functional modules and interconnects
parasitic effects of modern VLSI design aggravates such simulation tasks,
rendering the transient simulation of the full system very inefficient.
Therefore, accurate yet fast circuit simulation methods have always been 
desired in the IC industry. 
\end{comment}
\IEEEPARstart{M}{odern} VLSI design verification relies 
heavily on the analysis of power  delivery  network (PDN)
to estimate power supply noises \cite{Najm03_DAC, Nassif00_Fast,Brooks07_DATE,Lin04_CICC,Lin01_ICCAD,Skadron12, Skadron14_ISCA,Skadron14_ASPDAC}. 
%It is crucial to cost-efficient production of VLSI chips. 
The performance of power delivery network highly impacts on the quality of 
global, detailed and mixed-size placement \cite{ljw14_DAC, Pan05_ICCAD, ljw15_TCAD},  
clock tree synthesis \cite{Xiao10_ICCAD}, global and detailed routing \cite{Zhang12_ICCAD}, as well as timing \cite{Kahng13_ICCAD} and power optimization.  %\cite{ljw12_TVLSI}.
Lowering supply voltages, increasing current densities as well as tight design margins demand more accurate large-scale PDN simulation. 
Advanced technologies \cite{Zhuo12_ICCAD,Zeng10_DAC}, three dimensional (3D) IC structures \cite{Zhuang13_ICCCAS,Lim14_ICCAD,xiang14}, 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 \cite{Zhuang15_emc, Zhuang15_DAC}, PDN is often modeled as a large-scale linear circuit with voltage supplies and time-varying current sources \cite{Nassif08_Power, zuochang08, Li12_Tau}. Those linear matrices are obtained by parasitic extraction process \cite{Zhuo14_DAC, Zhuang2012, Yu14_TCAD, Yu2013, Lin04_CICC}. 
After those processes, we need time-domain 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 \cite{Leon75}.
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.
%It also leads to the unstable convergence rate of iterative methods \cite{Saad03}, 
Due to the requirement of a robust solution,
compared to iterative methods \cite{Saad03},
direct methods \cite{Davis06} are often favored for VLSI circuit simulation, and thus adopted by state-of-the-art power grid (PG) solvers in TAU PG simulation contest \cite{Yu12, Yang12, Xiong12}. 
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 \cite{Xiong12, Yu12, Li12_Tau} in their implicit numerical integration framework.  
However, the maximum of step size choice is limited by the smallest distance  $h_{upper}$ among the breakpoints \cite{Nagel75}. Some engineering efforts are spent to break this limitation by sacrificing the accuracy. 
In our work, we always obey the upper limit $h_{upper}$ of time step to maintain the fidelity of model, which means the fixed time step $h$ cannot go beyond   $h_{upper}$ in case of missing breakpoints. 


%The size of time step is pre-defined
%The fixed step size is limited by the local tran

Beyond traditional methods, 
a  class of methods called matrix exponential
time integration has been embraced by MEXP \cite{Weng12_TCAD}.
The major complexity is caused by matrix exponential computations. 
MEXP utilizes standard Krylov subspace method \cite{Saad92} to approximate matrix exponential and vector product. MEXP can solve the DAEs with much higher order polynomial approximations than traditional ones \cite{Weng12_TCAD, Saad92}.
%Another merit of matrix exponential based circuit simulation is that, for linear circuit,  it can proceed the adaptive time stepping without factorization of matrices on-the-fly, while the traditional counterparts cannot avoid such time-consuming process during the adaptive time marching. 
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.


\begin{comment}
$\bullet$ Step size of fixed time step PG solvers \cite{Yang12, Yu12, Xiong12} is limited by the minimum size of transition distance along a input source. Meanwhile, using adaptive stepping causes very expensive matrix re-factorization overhead.

$\bullet$ At each transition of input sources, system states  are required to compute, which limits the step size of adaptive stepping in MEXP. 
%regenerate Krylov subspace basis.
$\bullet$ 
\end{comment}

Nowadays, the emerging multi-core and many-core platforms bring powerful computing resources and opportunities for parallel computing. Even more, cloud computing techniques \cite{cloud} drive distributed systems scaling to thousands of computing nodes 
\cite{mapreduce, xiang09, spark10}, etc. 
Distributed computing systems have been incorporated into products of many leading EDA companies and in-house simulators \cite{He14_emc, Jia15_TCAD, Wang13, Tan13_ICCAD, Feng08}. 
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 \cite{Zhuang14_DAC, gander2013paraexp} show  great potentials by parallelizing matrix exponential based method to achieve the runtime performance improvement and maintain high accuracy. 

%Previous works \cite{Li11, Weng12_ICCAD, Weng12_TCAD, Wang13, He14_emc} have been made in order to improve circuit simulation by novel algorithms,  parallel processing and distributed computing.

 
In this work, we develop a transient simulation framework  
using matrix exponential integration scheme, \emph{MATEX}, for  PDN simulation.
Following are the challenges we need to address. 
%We address following challenges. 
% are required to solve for the efficiency. 
First,   when the circuit is stiff, the standard Krylov subspace has convergence problem and slows down the computation of MEVP.
% is slow by standard Krylov subspace method. 
Second, the frequent time breakpoints due to the transitions of PDN current sources modeling  triggers the generations of Krylov subspace. % with the capability of large time stepping,, which, even for small time step size. 
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:
\begin{itemize} 
\item   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. 
\item 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  {%\color{red}
up to around $15\times$} runtime speedup  against the benchmarks over the traditional method  TR-FTS with good accuracy. 
%\item[(2)] HR-MATEX is proposed to  improve the performance of single computing mode MATEX by removing  Krylov subspace generations for small time stepping. HR-MATEX combines low order scheme (Trapezoidal method), so that we can have fast time stepping for small time steps instead Krylov subspace method. For large time marching, HR-MATEX reverts to original MATEX easily without any  overhead. This shows the traditional scheme can actually be a subset of our framework. This method gains $13.4\times$ runtime speedups over traditional TR method with fix step size.
\item Furthermore, DR-MATEX is designed to improve R-MATEX with distributed computing resources.% { DR-MATEX has the embarrassingly parallel computation model.} 
\begin{itemize}
\item 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.
\item
After all nodes finish the simulation computations, the results are summed up based on the linear superposition property of the PDN system.
\item 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\times$ speedup over traditional method TR-FTS. 
\end{itemize}
%\item Our DR-MATEX has an unique feature that we do not need to communicate among different nodes during the transient simulation process. It differs from the relevant work proposed by Wang and Xiong \cite{Wang13}, which has communication during the each stages during the transient simulation.   
\end{itemize}  
% In addition,  our framework gains  runtime speedups up to $50\times$ in the environment of single computing node, and $88\times$ in the distributed computing environment  
%in transient computing for the power grid simulation benchmarks. 

% \emph{Our Contributions}. 
%efficiently parallelize the simulation and bring huge speedups than the traditional methods, which are constraint by the number of time steps.
%with variant Krylov subspace-based methods. 
%to achieve large speedups over previous MEXP. 
%To accelerate matrix exponential computation in stiff circuit simulation,  
%We generalize  \cite{Zhuang13_ASICON}
%with rational Krylov subspace 
%\cite{Zhuang13_ASICON}MATEX can efficiently distributed the simulation and achieve huge speedups compared to traditional algorithmic framework using fixed time step method


\begin{comment}
\textbf{MATEX leverages the superposition property of linear system.
Through decomposing the input transition points, MATEX assigns subsets of 
input transition points to different computing nodes and  process them in parallel}.
In each computing node, the assigned input transitions decomposition \emph{locally} reduces the chances of Krylov subspace basis generations.
%and provides longer time regions (before the next transition point).
%Thus, the circuit solver can adopt large adaptive time stepping,
The intermediate time points, which is used to form the \emph{global solutions} via superposition, are obtained by computing matrix exponentials with scaled latest Krylov subspace for different time steps, instead of generating new one. 
% Therefore, re-computing the whole new one.
In contrast, traditional methods still need to derive solutions step by step, which includes costly matrix factorizations for the adaptive stepping methods or backward and forward substitutions for the fixed-step methods. 
%either performing backward and forward substitutions 
%or adaptive stepping with solving the linear system using direct or iterative methods. 
Thus, they cannot benefit from such input transitions decomposition. 

\textbf{  %Overcome the stiffness impediment of MEXP then form an efficient circuit solver within MATEX
Through incorporate inverted and rational Krylov subspace bases into MATEX,
we tackle the stiffness problem to MEXP and then form an efficient circuit solver within MATEX.
}
(1)
MATEX's circuit solver efficiently reduces the dimension of basis for matrix exponential and vector product, 
remarkably eases the memory consumptions 
and boost computational performance. 
(2) After unleashing from the constraint of stiffness, MATEX uses larger time size to expedite the whole time domain simulation.
%When dealing with singular capacitance/inductance matrix, 
(3) When dealing with the singular capacitance/inductance matrix in DAE, 
the regularization process \cite{Chen12_TCAD, Weng12_TCAD} is required 
in MEXP\cite{Weng12_TCAD} to 
%The regularization process brings extra overhead to 
make this matrix nonsingular, which brings extra computational overhead.
However, MATEX does not need this process. It is because the two new adopted Krylov subspaces do not require the direct factorization of this singular matrix.
\end{comment}
%\emph{Paper Organization}.
%The rest of paper is organized as follows.
The rest of this paper is organized as follows. Sec.  \ref{sec:method} introduces the background of
linear circuit simulation and matrix exponential formulations.
%Sec. \ref{sec:matex} presents overall framework of MATEX.  
Sec. \ref{sec:advanced} illustrates the Krylov techniques to accelerate matrix exponential and vector product computation.   
Sec. \ref{sec:matex} presents MATEX circuit solver and the parallel framework DR-MATEX.
Sec. \ref{sec:exp_results} shows numerical results 
and Sec. \ref{sec:conclusion} concludes this paper.


\begin{comment}
and its adaptive time controlling into MATEX. 
During the transient simulation, 
the input changes cause the updates of matrices. 
However, the updated region resides within only columns of the matrix.
By precomputing block LU calculation beforehand and 
updating that region step by step, MATEX does not require LU factorizations on-the-fly, 
which is very different from that of TR-like methods.
%{\color{red} Only by updating small part of matrices,
%due to the input change,
%MATEX's adaptive time stepping does not require 
%LU factorization on-the-fly.}

We also combine MATEX %together 
with \emph{inverted Krylov subspace} as an alternative technique,
which only requires decomposing the conductance matrix. 
Therefore,
it is advantageous the avoidance of inverting the 
matrix of capacitance and inductance, which may 
include many parasitic energy elements 
with large magnitude variations
after circuit extraction.
%{\color{red}{The connections of these two methods are inspired by the works 
%from model order reduction community, \emph{i.e.,}
%PRIMA \cite{PRIMA98} and SPRIM \cite{SPRIM04}.}}

%These solve the limitations (1) and (2);
Additionally, the parallelism of MATEX
is based on the superposition property of linear circuit.
Within our theoritically sound skeleton,
we split the input sources, use all of divided subset of input sources 
to stimuli the system, seperately, %and efficiently adaptive time stepping,
making the whole simulation task divided
into serveral sub-tasks. 
% and solved seperately 
These sub-tasks can be solved parallel, without sacrificing any accuracy. 
This circumvents the limitation (3).
Besides, the traditional methods (FE, BE, and TR) cannot benefit from this
splitting, due to their limitations not from the side of input sources
but from LTE and stability.
With the advent of inexpensive computer clusters, and multi/many-core technology,
this distributable algorithmic framework is a nice one for future implementations.


\begin{comment}


The major reason is that mainstream circuit simulators inherit implicit numerical
integration, e.g. backward Euler, Trapezoidal methods, 
for the sake of stability, where the explicit methods suffer. 
However, implict methods trouble the perspective of scalability
due to the requirement 
of solving large linear systems of differential equations.
On the contrary, the explict methods are more scalable and parallelizable 
for current multi-core and manycore architecture.
The merit of explict methods resides in simple time propogation scheme via 
matrix-vector multiplication who costs relative lower
computation and memory complexity about $O(n)$.
Even though the existence of stability issue, 
explicit methods still become a revived research 
area and verifies their parallelism and scalability 
within modern parallel computer architecture
~\cite{Dong08}\cite{Dong09}\cite{Weng12_ICCAD} during these years. 

One of promising and distinct method for circuit simulation
is matrix exponetial method (MATEX) proposed in 
\cite{Weng12_TCAD}. Apart from the above numerical methods, 
the DAE/ODE system (in one discritized region of certain time step) can  
be solved in analytical way via matrix expoential operator $e^{\mbf A}, 
\mbf A \in \mbf R_{n \times n}$. 
By the analytical property, this distinct class 
of numerical approach removes the local truncation  
and it is A-\emph{stable} for passive circuits. 
Therefore, it is not constraint by the stablity problem. 
{\color{red}{This property matters!}}
%(When $h >  0, h \to \infty$, the $x(h) \to 0$). Very unlike forward
%Euler method that require small time step to let fast damp compoent
%die down first, bring too much time steps. But M
However, the kernel computation of MATEX lies in computing matrix expoential is very high,  
e.g., direct computation consumes around $O(n^3)$ complexity. 
Weng \emph{et al.} \cite{Weng12_TCAD}\cite{Weng12_ICCAD} introduced a \emph{standard Krylov subspace}
method-based MATEX for circuit simulation to solve this problem. 
In addition, unlike the implicit methods (Trapezoidal and Back Euler) 
they both require fixed step size so that they reduce the LU online for linear circuit
transient simulation,
this MATEX can utlize the scaling invariant property of Krylov subspace to do the adaptive time-steping 
without compromising any accuracy.
(1) MATEX can jump further than BE and Trap, because the exponenital 
denomator in the exponential approximation, it is much better
than the Taylor expansion.
(2) Under the (PWL) assumptions, MATEX can also carry 
the input information during the stepping.  
(3) With the help of Krylov subspace scaling invariant, it can re-evaluate
the inter time point without rebuilding Krylov subspace.

%By this strategy, it helps to capture the low frequency response from circuit system and further speedups the whole simulation without fixed step size limitation.
% In recent years, the Krylov subspace method has
Nonetheless, the step size is actually limitted by the stiffness of matrix when utilizing standard Krylov subspace.
It is because the large spectrum of $\mbf A$ and the Arnoldi process focus on capturing %the extreme magnitude
% and
well-distributed eigenvalues, so it ofter misses the dominant factors for the circuit system/for
matrix epxonential operator, which are the ones with smallest magnitudes.
Therefore, the Arnoldi process requires large number of order to preserve the 
accurate of matrix exponential and vector calculation. 
Also, the standard Krylov-Based MATEX requires the process of regularization \cite{Chen12_TCAD} 
to deal with singular matrix $\mbf C$, i.e., nodes without grounded capacitance or 
inductance/voltage branches. Both drag the overall performance of MATEX for circuit simulation.  
% former MATEX~\cite{Weng12_TCAD} requires regularization. 
% In addition, when the matrix $\mbf C$ is stiff, it requires smaller 
%time steps or even larger Krylov subspace 
% order $m$, which degrades the method. 


In this paper, we integrated our previous works on MATEX for linear circuits 
using \emph{rational Krylov subspace}~\cite{Zhuang13_ASICON} and \emph{inverted Krylov subspace}.
%The rational Krylov subpsace uses 
They use 
$(\mbf I-\gamma \mbf A)^{-1}$ and $\mbf A^{-1}$, respectively, 
as the basis instead of $\mbf A$ (\emph{standard Krylov subspace}) in \cite{Weng12_TCAD}. 
The rational basis limits the spectrum of a circuit matrix into a unit circle,
and the inverted operation puts the important eigenvalues with smallest magnitudes from former spectrum
into a large magnitude ones. By using Arnoldi algorithm, the Hessenberg matrix can resemble the former matrix
$\mbf A$ projected onto the new invariant subspace better and contain the more accurate important eigenvalues
for matrix exponential operator.
This method provides a good convergent rate in Arnoldi process within MATEX framework.
The inverted Krylov subspace's convergent rate in Arnoldi process may not be as good as rational basis. However,
it only requires the LU factorization of conductance matrix $\mbf G$. It provides
a strategy dealing with $\mbf C$ is very complicated from post-layout extraction, but
$\mbf G$ is much more simple then. Similarly as rational basis, by the invertion $\mbf A$, the small magnitude eigenvalues 
become the large ones in the transformed spectrum, so that the exponential of a matrix can be accurately approximated.
A very typical and important application for circuit simulation is power grid analysis, which plays an 
important role in IC deisgn methodology. Given current stimulus and power grid structure, designers should verify
and pedict the worst-case voltage noise through the simulation before signing off their design. 
The model order reduction provides good way to reduce matrix of circuit system, e.g. 
PRIMA \cite{PRIMA98} and SPRIM \cite{SPRIM04}. However, with the large number of 
input sources, this kind of strategy cannot benefit from 
the reduction models any more. 
The majority works still use Trapezoidal method with fixed step size
for the power grid simulation.
Till now, accourding to Tau power grid simulation contest 2012 \cite{Li12_Tau}\cite{Yu12}, 
the trapezoidal method with fixed step size using one LU and 
many backward/forward substitutions lead to the best performance. 
The small and fixed time steps bring benefits for the transient power grid analysis 
when dealing the voltage variation 
with high frequency, i.e. the input is changing rapidly or constantly. 

However,
when package has an important impact on accurarte 
on-die power grid analysis~\cite{Pant06}\cite{Pant07}, 
and middle and low frequency phenomena is of interests to the designers, 
trapezoidal method with fixed time step becomes very inefficient. 
It is 
because such mehtod takes long time to represent the steady state of power grid system.
In the architecture analysis, power distribution networks (PDN) 
involves with pakcage and board 
level simulations \cite{PengLi12_ICCAD}\cite{Pant06}\cite{Pant07},
which requires longer timespan of simulation,
e.g. up to few $\mu s$, the simulation process are 
often time-consuming with the limitation of fixed step size to avoid online LU.
Our method can exploit large time step, ``steady region'' of waveform
during the simulation, meanwhile, it only requires one LU factoraization
as trapezoidal method, with backward/froward substitutions.
With equivalent error buget, the krylov subspace MATEX can jump larger time step/fewer 
number of time steps
than trapezoidal
method. Therefore, even though such fully implicit scheme can be used here,
when a sparse LU factorization is affordable,
the simulation burden is higher than our Krylov subspaces MATEX because of 
more time steps need.

The aim of this paper is to present and organize the Krylov subspace based 
MATEX solver in transient simulation of linear circuit.
We list our main contributions of this paper are as follows.
\begin{itemize}
%\item A breakthrough of the circuit simulation matrix exponential.
\item Accelerate MATEX via rational and inverted Krylov 
subspace.
\item Using adaptive time 
stepping strategy for circuit simulation and overcome the 
problem to achieve low-frequency pheonomena of circuit system.
\item Eliminating the regularization process 
for singular matrix $\mbf C$ for explicit method 
of circuit simulation.
\item By input grouping and superposition property of linear circuit 
to further accelerate circuit simulation.
\item Establish the relation of MATEX with previous works on Model Order Reduction. 
\end{itemize}

By such contributions, we are able to achieve hihgly efficient trasient 
simulation. Compared to the works~\cite{Weng12_TCAD}.
It is very nice property from linear system, which they can be 
superposition. 
By this idea, the same system with the parts of inputs send to different
machines, cores, or threads and add together without communication 
during the process, only require on precomputed time stamp/anchors.
It is a nice parallelism framework, we provide such for the further 
extension.
It is noted that Trapezodal cannot achieve this benefit, because
it is not analytical solution, and the inputs are not the 
hurdles during their simulation but the stability and local trancation
error are.

The paper is organized as follows. Section II is the fundemantals 
of circuit simulation and matrix exponential inteperation.
Section III describes the Krylov methods for fast matrix exponential evaluations. 
They are three algorithms for such computation, Standard Krylov subspace, inverted 
Krylov subspace, Rational krylov subspace.
Section IV presents the computational results of simple linear circuit cases, also contains industrial power grid designs and the famous IBM benchmark suite.

\begin{comment}

\begin{itemize}
\item Accelerate MATEX via rational and inverted Krylov subspace and its coresponding adaptive time stepping strategy for circuit simulation.
%    \item And, summarize three Krylov strategies for MATEX, the standard/inverted/rational Krylov subspace.
%  \item Use inverted Krylov subspace, enable the krylov subspace based MATEX without regularization~\cite{Chen12_TCAD}.
\item Eliminating the regularization process for singular matrix $\mbf C$ for explicit method of circuit simulation.
\item Speedups linear circuit simulation, by adaptive time stepping.%  partitioning input sources, and using superposition,
    % \item Using block matrix for Krylov subspace generation.
    % \item Establish the relation of MATEX with previous works on Model Order Reduction. 
    % As a explicit circuit solver, which is also enjoy the model order reduction framework.
    %    \item Generalize the explict method for future parallelism circuit simulation (Stablity).
\end{itemize}
For solving ordinary differential equations $x'=Ax+f(t)$,
it has an anlytical solution $x(t+h)=e^{Ah}x(t)+ \int_{0}^{h}e^{A(h-\tau)}f(t+\tau)d\tau$. 
When $A$ is large, it is very expensive to cacluate the exponential operator directly.  
The explict and implict methods are strategies to approximate the operation of $e^{Ah}$.  
Explicit methods required matrix-by-vector multiplications, however, implicit methods need to solve 
a whole linear system.  
It is known that explicit methods are constraint on the size of the time steps required to make sure stability. 
Compared to explicit methods, Implicit methods have more difficulties on utilizing current 
parallel and distributed machines.  
Recent development of Krylov subspace accelerates the exponential operations.  
The  
\begin{itemize}
\item without regularization, 
    \item Accelerate and enpower MATEX framework via rational Krylov subspace algorithm.
    \item Enabling adaptive time stepping for circuit simulation.
    \item Establish the relation of MATEX with previous works on Model Order Reduction. As a explicit circuit solver, which is also enjoy the model order reduction framework.
    \item ?? Generalize the explict method for future parallelism circuit simulation (Stablity).
\end{itemize}
The MATEX has a good characteristic that the exponential operator bares the propagation in the perspective of
time derivitate,
which the method at Forward Euler (FE), Backward Euler (BE) and Trapezoidal (TRAP) method use time step $h$ to discretize the
dynamical system with respect to time. 
They lead to explicit and implicit methods,
which bring the stability and limitted time step due to the system, e.g. 
stiff systems.
Explicit methods are relatively cheap computation; Implicit methods are usually time-consuming because they require
to solve the large linear systems. For this regard, in power grid simulation, it precomputes the matrix 
($\mbf C/h+ \mbf G$) in BE, or ($\mbf C/h + \mbf G/2$) in TRAP. This involes the time step $h$. 

Also, in linear circuit, the the mexponential operator is a explicit method, of which the matrix-vector multiplication
take the major part of computation.



Power grid simulation plays an important role in IC design methodology. Given current stimulus and the power grid structure, 
designers could verify and predict the worst-case voltage noise through the simulation before signing off their design.  
However, with the huge size of modern design, power grid simulation is a time-consuming process. Moreover, manifesting effects 
from the package and the board would require longer simulation time, e.g., up to few $\mu s$, which worsens the performance 
of the power grid simulation. Therefore, an efficient power grid simulation is always a demand from industry. 

Conventionally, the power grid simulation is based on the trapezoidal method where the major computation is to solve a linear 
system by either iterative approaches\cite{Saad03} or direct methods
%\cite{Davis06}. 
Many previous 
works~\cite{Chen01_Efficient}\cite{Feng08}\cite{Kozhaya02}\cite{Nassif00_Fast}\cite{Su03} utilize iterative approaches to tackle 
the scalability issue and enable adaptive stepping in the power grid simulation. However, the iterative methods usually suffer
from the convergence problem because of the ill-conditioned matrix from the power grid design. On the other hand, the direct 
methods, i.e., Cholesky or LU factorizations, are more general for solving a linear system. Despite the huge memory 
demanding and computational effort, with a carefully chosen step size, the power grid simulation could perform only one 
factorization at the beginning while the rest of operations are just backward/forward substitutions. Since a power grid design 
usually includes board and package models, a long simulation time is required to manifest the low-frequency response. Hence, 
the cost of expensive factorization can be amortized by many faster backward/forward substitutions. Such general factorization 
and fixed step size strategy~\cite{Xiong12}\cite{Yang12}\cite{Yu12}\cite{Zhao02} is widely adopted in industry.  

%Beyond the traditional trapezoidal method, 
Recently, Weng \emph{et al.}~\cite{Weng12_TCAD, Weng12_ICCAD} introduced 
the matrix exponential method (MATEX) for the circuit simulation recently. MEXP has better accuracy and adaptivity because of the 
analytical formulation and the scaling invariant Krylov subspace approximation. Unlike the fixed step size strategy, MATEX could 
dynamically adjust the step size to exploit the low-frequency response of the power grid without expensive computation. However, 
the step size in MATEX is limited by the stiffness of circuit. Furthermore, MEXP has to regularize~\cite{Chen12_TCAD} the 
singularity in power grid, i.e., nodes without grounded capacitance or inductance/voltage branches. Both constraints would drag 
the overall performance of MATEX for the power grid simulation. 

In this work, we tailor MATEX using \emph{rational Krylov subspace} for the power grid simulation with adaptive time stepping. 
The rational Krylov subspace uses $(\mbf{I}-\gamma \mbf{A})^{-1}$ as the basis instead of $\mbf{A}$ used in the conventional Krylov 
subspace. The rational basis limits the spectrum of a circuit so that the exponential of a matrix can be accurately approximated 
even under a large step size. Moreover, the rational avoids $\mbf{C}^{-1}$ and thus eliminates the regularization process. 
As a result, MATEX with rational Krylov subspace is free from the stiffness constraint and the regularization process. The 
simulation can purely enjoy benefits of the adaptivity and the accuracy of MATEX. Even though the rational Krylov subspace still needs 
to solve a linear system as the trapezoidal method does, MATEX can factorize the matrix only once and then constructs the rest of 
rational Krylov subspaces by backward/forward substitutions. Therefore, MATEX can utilize its capability of adaptivity to accelerate 
the simulation with the same kernel operations as the fixed step size strategy. Overall, our MATEX has the following advantages: 
\begin{itemize}
    \item enabling adaptive time stepping for the power grid simulation with only one lu factorization.
    \item eliminating the regularization process for the singular power grid design.  
    \item allowing scaling large step size without compromising the accuracy. 
\end{itemize}
The experimental results demonstrate the effectiveness of MATEX with adaptive step size. The power grid simulation for IBM benchmark 
suite~\cite{Nassif08_Power} and industrial power grid designs can be accelerated 6X and 15X on average compared to the trapezoidal method. 

The rest of paper is organized as follows. Section~\ref{sec:preliminary} presents the background of the power grid simulation and MATEX. 
Sections~\ref{sec:rational_krylov} and~\ref{sec:adaptive_step_control} show the theoretical foundation of the rational Krylov subspace 
and our adaptive step scheme for the power grid simulation, respectively.  Section~\ref{sec:exp_results} presents experimental results 
of the IBM benchmark suite~\cite{Nassif08_Power} and several industrial power grid designs. Finally, Section~\ref{sec:conclusion} concludes 
this paper. 
%Specifically, the contributions of this work are:
%\begin{itemize} 
%    \item Enabling adaptive stepping in the power grid simulation with the computational cost of only one matrix factorization 
%	  and backward/forward substitutions. 
%    \item Using rational Krylov subspace method to ease the stiffness constraint~\cite{Weng12_TCAD}\cite{Weng12_ICCAD}, 
%	  meanwhile, endowing MATEX with the capability of dealing with singular $\mbf C$ without time-consuming 
%	  regularization process.
%\end{itemize}
%
\end{comment}


krylov.tex
In PDN simulation, $\mbf A$  is usually above millions and makes the direct computation of matrix exponential $e^{\mbf A}$ infeasible. The alternative way to compute the product is through Krylov subspace method \cite{Saad92}. 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. 

\subsection{MEXP: MEVP Computation via Standard Krylov Subspace Method}
\label{sec:std}
The complexity of $e^{\mbf A}\mbf v$ can be reduced 
using Krylov subspace method and still maintained  
in a high order polynomial approximation \cite{Saad92}.
%which has been deployed by MEXP 
%In this paper,
MEXP \cite{Weng12_TCAD} uses standard Krylov subspace, 
which uses $\mbf A$ directly to generate subspace basis through Arnoldi process (Algorithm  \ref{algo:arnoldi}). 
First,
we reformulate Eq. (\ref{eqn:pwl_exact_sol}) into
\begin{comment} 
\begin{eqnarray} 
\mbf{x}(t+h)
&=&
\mbf{y}(t+h)
-(\mbf A^{-1}\mbf b(t+h) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h})
\nonumber \\
\mbf{y}(t+h)
& =& e^{h\mbf{A}}\mbf{y}(t)	%	\\
\end{eqnarray}
\end{comment}
\begin{eqnarray}
	\label{eqn:new_exact}
	\mbf{x}(t+h)
	=
	e^{h\mbf A}(\mbf x(t) + \mbf F(t,h) ) - \mbf P (t,h),
\end{eqnarray}
where 
\begin{eqnarray}
	\mbf F(t,h) =  \mbf A^{-1}\mbf b(t) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h}
\end{eqnarray} and 
\begin{eqnarray}
	\mbf P(t,h) =   \mbf A^{-1}\mbf b(t+h) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h}  .
\end{eqnarray}
%$\mbf y(t) = 
%\mbf{x}(t) + \mbf A^{-1}\mbf b(t) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h}$.
The standard Krylov subspace 
\begin{eqnarray}
	{\rm K_m}(\mbf A, \mbf v) := \text{span}
	{\{ \mbf v, \mbf A\mbf v,
		\cdots, \mbf A^{m-1} \mbf v}\}
\end{eqnarray}
obtained by Arnoldi process
has the relation
\begin{eqnarray}
	\mbf{A}\mbf{ V_m} = 
\mbf {   V_m} \mbf {  H_m}
	+ {  h_{m+1,m} }
	\mbf {  v}_{m+1}\mbf e_m^\mathsf{T}, 
\end{eqnarray}
{ where %$\bar h_{m+1,m}$ is the $(m+1,m)$ entry of 
$\mbf {H_m}$ is the upper Hessenberg matrix 
\begin{eqnarray}
\mbf {H_m } = \begin{pmatrix}
 h_{1,1} &   h_{1,2} & \cdots &   h_{1,m-1} &   h_{1,m} \\
 h_{2,1} &   h_{2,2} & \cdots &   h_{2,m-1} &   h_{2,m} \\
0 &   h_{3,2} & \cdots &  h_{3,m-1} &   h_{3,m} \\
\cdots & \cdots & \cdots & \cdots  & \cdots \\
0 & 0  & \cdots &  h_{m,m-1} &   h_{m,m} \\
\end{pmatrix},
\label{eq:hess}
\end{eqnarray}
$\mbf { V_m}$ is a $n \times m $ matrix by $(\mbf v_1, \mbf v_2, \cdots, \mbf v_m)$, and $\mbf e_m$ is the $m$-th unit vector with dimension $n\times 1$.
MEVP is computed via 
\begin{eqnarray}
	e^{h\mbf A} \mbf v \approx \beta \mbf {  V_m} e^{h\mbf {  H}_m} \mbf e_1 .
\end{eqnarray}
%$\mbf {H_m}$ is smaller than $\mbf A$. 
The posterior error term is  
\begin{eqnarray}
	\label{eqn:std_error}
   r_m(h)  =  \lVert \beta { h_{m+1,m}} \mbf v_{m+1}
  	\mbf e^{\mathsf{T}} e^{h\mbf{  H}_m} \mbf e_1 \rVert,
\end{eqnarray}
% the matrix $\mbf C$, 
where $\beta = \lVert  \mbf v \rVert $.
However, for an autonomous system $ \mbf C \dot {\mbf x}(t) = - \mbf G \mbf x(t)$  in circuit simulation, we consider the residual between  
$\mbf C \dot {\mbf x}(t) $ and $-\mbf G \mbf  x(t)$, 
which is 
$$ \mbf C \dot {\mbf x}(t) + \mbf G \mbf  x(t) ,$$
instead of  $$ \dot {\mbf x}(t) - \mbf A \mbf  x(t)$$  in
$ \dot {\mbf x}(t) = \mbf A \mbf x(t)$.
This leads to 
\begin{eqnarray}
	\label{eqn:std_err}
	   r(m,h)  =    \lVert\beta  h_{m+1,m}  \mbf C  \mbf 
v_{m+1}	\mbf e^{\mathsf{T}} e^{h\mbf{  H}_m} \mbf e_1 \rVert
\end{eqnarray}
and helps us mitigate the overestimation of the error bound. 
 %, especially for the singular $\mbf C$ system.

To generate $\mbf x(t+h)$ by Algorithm \ref{algo:arnoldi}, 
we use 
\begin{eqnarray}
\label{eq:lu}
[\mbf L, ~\mbf U] = \text{LU\_Decompose}(\mbf X_1),
\end{eqnarray}
where $$\mbf X_1 = \mbf C ~ \text{and} ~ \mbf X_2 = \mbf G$$ as inputs for standard Krylov subspace.
% In the line \ref{alg1:h} of Alg. \ref{algo:arnoldi}, $\mbf H_m = \mbf {\bar H}_m$  or $  h_{i,j} =   {\bar h_{i,j}}$.
%Besides,
The error budget $\epsilon$ and Eq. (\ref{eqn:std_error}) are used to determine the convergence when  time step is $h$ and Krylov subspace dimension is $j$  (from line \ref{alg1:conv1} to line \ref{alg1:conv2} in Algorithm \ref{algo:arnoldi}). 
%{\color{red} algorithm recap, the on-the-fly krylov basis desicion }
\begin{algorithm}
	\label{algo:arnoldi}
	%\caption{Arnoldi process in Matrix Exponential-Based Circuit Simulation with Krylov subspace Method}
	\caption{MATEX\_Arnoldi}
	%$\rm K_m(\mbf X_1^{-1}\mbf X_2, \mbf v)$} %~\cite{Weng12_TCAD, Weng12_ICCAD}}
	\KwIn{ $\mbf L, \mbf U, \mbf X_2, h ,t, \mbf x(t),  \epsilon, \mbf P (t,h),\mbf F (t,h)$  }
	%\KwOut{ $\mbf V_m, \mbf H_m,  m$ }
	\KwOut{ $\mbf x(t+h), \mbf V_m, \mbf H, \mbf v$}
	{   
		%	[$\mbf L, \mbf U$]=LU($\mbf X_1$)\\
		%	Normalize the intial vector 
		$\mbf v = \mbf x(t) + \mbf F(t,h)$\;
		$\mbf v_1 = \frac{\mbf v }{\lVert \mbf v \rVert}$\;
		\For  {$j=1:m$}
		{
			% Calculate 
			%	    $   \mbf w =  \mbf X_1^{-1} \mbf X_2 \mbf v_{j} $ by solving
			\label{arnoldi_orth} 
			$\mbf w = \mbf U \backslash  (\mbf L \backslash (\mbf X_2 \mbf v_{j}))$
			{\tcc*[l]{a pair of forward and backward substitutions.}}
			\For {$i = 1:j$}
			{
				%	Orthogonalization 
				$h_{i,j} = \mbf w^T\mbf v_{i}$\; 
				$\mbf w = \mbf w - h_{i,j} \mbf v_{i}$\;
			}
			$h_{j+1,j} = \lVert \mbf w \rVert $\; 
			$\mbf v_{j+1} = \frac{\mbf w }{h_{j+1,j}} $\;
			% \If {$ h_{i+1,j} < \epsilon $}
			\If {$r(j,h)  < \epsilon $\label{alg1:conv1}} 
			{
				$m = j$\;
				break\; 
			}\label{alg1:conv2}
		}
		%	$\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V_m} e^{h \mbf H_m} \mbf e_1 $\;
		$\mbf x(t+h) = \lVert \mbf v \rVert \mbf {V_m} e^{h \mbf H} \mbf e_1 - \mbf P(t,h)$\;
		\label{alg1:h}
		%$\lVert r_m(h) \rVert  =  \lVert \mbf v \rVert \left| h_{m+1,m} e_m^T e^{h\mbf H_m} e_1 \right|$
	}
\end{algorithm}

%\subsection{Discussions of MEXP}
\label{sec:discuss}
%There are two main challenges in MEXP for transient simulation of linear circuit. 
%, \emph{input sources} and \emph{Krylov subspace dimension}.
%First, 



The standard Krylov subspace may not be efficient when simulating stiff circuits \cite{Weng12_TCAD, Weng12_ICCAD}.
%matrix exponential approximation 
For the accuracy of approximation of $e^{\mbf A}\mbf 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 reason is that the Hessenberg matrix $\mbf H_{m}$ of standard Krylov subspace tends to approximate the large magnitude eigenvalues of $\mbf A$ \cite{Van06}. 
Due to the exponential decay of higher order terms in Taylor's expansion, such components are not the crux of circuit system's behavior 
%\cite{Botchev12, Van06, Orec:12}.
\cite{Botchev12, Van06}.
%exponential approximation 
Therefore, to simulate stiff circuit, we need to gather more vectors into subspace basis 
and increase the size of $\mbf H_{m}$ to fetch more useful components, which results in both
memory overhead and computational complexity to Krylov subspace generations for each time step.
% during time stepping.
%Direct deploying MEXP into MATEX's \emph{Circuit Solver}  is not efficient to leverage the benefits of local flexible and larger time stepping.
In the following subsections, we adopt  ideas from \emph{spectral transformation} \cite{Botchev12, Van06} to effectively capture small magnitude eigenvalues in $\mbf A$, leading to a fast and accurate MEVP computation. 
% \emph{Circuit Solver} for MATEX.

%\vspace{-0.02in}
\subsection{I-MATEX: MEVP Computation via Invert Krylov Subspace Method}
\label{sec:imatex}
%Based on the rationale that standard Krylov subspace method
%tends to approximate large magnitude eigenvalues of $\mbf A$. 
%Then, 
Instead of $\mbf A$, we use $\mbf A^{-1}$ (or $\mbf G^{-1} \mbf C$) as our target matrix to form
%then recover it 
\begin{eqnarray}
	{\rm K_m}(\mbf A^{-1} , \mbf v) := \text{span}
	{\{ \mbf v, \mbf A^{-1}\mbf v, \cdots, \mbf A^{-(m-1)} \mbf v}\}.
\end{eqnarray}
Intuitively, by inverting $\mbf A$, the small magnitude eigenvalues become the large ones
of $\mbf A^{-1}$.
The resulting $\mbf H_{m}$ is likely to capture these
eigenvalues first. Based on Arnoldi algorithm, the invert Krylov subspace has the relation 
%in pratical computation.
\begin{eqnarray}
	\mbf{A}^{-1}\mbf{V_m} =  
	\mbf V_m \mbf H_m + 
	h_{m+1,m} \mbf v_{m+1}\mbf e_m^\mathsf{T}. 
\end{eqnarray}
The matrix exponential $e^{\mbf A} \mbf v$ is calculated as 
%$\left | \beta \right| \mbf V_m \mbf e^{h\mbf H'^{-1}_m}\mbf e_1$.
\begin{eqnarray}
	e^{\mbf A} \mbf v \approx \beta \mbf V_m e^{h\mbf {H}^{-1}_m} \mbf e_1.
\end{eqnarray}
To put this method into Algorithm 1 is just by modifying the inputs 
$\mbf X_1 = \mbf G$ for the LU decomposition in Eq. (\ref{eq:lu}),
%$[\mbf L,\mbf U] =  \text{LU\_Decompose}(\mbf G)$, 
and $\mbf X_2 = \mbf C$.  
In the line \ref{alg1:h} of Algorithm \ref{algo:arnoldi}, 
$$\mbf H = \mbf {H}^{-1}_m$$ for the invert Krylov version.
The posterior error approximation \cite{Zhuang14_DAC} is
\begin{eqnarray}
	\label{eq:err_inverted_krylov}
	    r_m(h) 
	= 
	\lVert \beta  h_{m+1,m}  \mbf A  \mbf v_{m+1} 
	\mbf e^\mathsf{T}_m {\mbf H}^{-1}_{m}
	e^{h {\mbf H}}
	\mbf e_1 \rVert,
\end{eqnarray}
which is derived from residual based error approximation in
\cite{Botchev12}.
However, as mentioned in Sec. \ref{sec:std}, we consider the residual of $\left(\mbf C \dot {\mbf x}(t) + \mbf G \mbf x(t) \right)$, 
instead of $\left( \dot {\mbf x}(t) - \mbf A \mbf  x(t) \right)$, 
which leads to 
\begin{eqnarray}
	\label{eq:er_inverted_krylov}
   r(m,h) 
	= 
	\lVert   \beta h_{m+1,m} \mbf G   \mbf v_{m+1} 
	\mbf e^\mathsf{T}_m {\mbf H}^{-1}_{m}
	e^{h {\mbf H}_{m}}
	\mbf e_1  \rVert .
\end{eqnarray}
{We use Eq. (\ref{eq:er_inverted_krylov}) for the line \ref{alg1:conv1} of Alg. \ref{algo:arnoldi}.
 %, which is $r(j,h)=r'(j,h)$.}
% It mitigates the overestimation of error bound.
%When the $\mbf A$ has large condition number and densely distributed spectrum,
%When $\mbf A$ has a large condition number, the convergence rate  might be slow 
%for matrix exponential and vector approximation. 
%In the next subsection, we solve this problem 
%by confining the specturm of
%using shift-and-invert Krylov subspace \cite{Van06, Moret04},

%\vspace{-0.02in}
\subsection{R-MATEX: MEVP Computation via Rational Krylov Subspace Method}
\label{krylov_rat_err}
The shift-and-invert Krylov subspace basis 
%\cite{Van06, Moret04} 
\cite{Van06} 
is designed
to confine the spectrum of
$\mbf A$.
%\begin{eqnarray}
%    \label{eq:rat_kry}
Then, we generate Krylov subspace via
\begin{eqnarray}
	& &  \mbf {K_m} ((\mbf I- \gamma \mbf {A} )^{-1}, \mbf v) :=  \\
	& &  \text{span} \{ \mbf v, (\mbf I- \gamma\mbf { A} )^{-1}\mbf v,
	\cdots,
	%  \nonumber \\
	(  \mbf I -  \gamma \mbf{A})^{-(m-1)} \mbf v \}, \nonumber 
\end{eqnarray}
%\end{eqnarray}
where $\gamma$ is a predefined parameter.
With this shift, all the eigenvalues' magnitudes are larger than one. 
Then the inverse limits the magnitudes smaller than one.
According to \cite{Botchev12, Van06},
the shift-and-invert basis for matrix exponential based transient simulation
is not very sensitive to $\gamma$, once it is set to around the order near time steps 
used in transient simulation.
%and the $\gamma$ can be set as $h/10$,
%where $h$ is the fixed step size or resolution waveform.
The similar idea has been applied to simple power grid simulation 
with matrix exponential method
%called as rational Krylov subspace based MEXP 
\cite{Zhuang13_ASICON}.
Here, we generalize this technique and integrate into MATEX.
The Arnoldi process constructs $\mbf{V_m}$ and $\mbf{H_m}$. 
We have 
\begin{eqnarray}
	(\mbf{I}-  \gamma \mbf {A})^{-1} \mbf{ V}_m =
	\mbf{ V}_m\mbf{ {H}}_m +  {h}_{m+1,m} \mbf { v}_{m+1} \mbf e^\mathsf{T}_m . 
\end{eqnarray}
%    \label{eq:rat_kry2}
%\end{eqnarray}
%where $e_m$ is the $m$-th unit vector with dimension $m \times 1$. 
%Matrix $\mbf{H_m}$ is an upper Hessenberg matrix, and 
%$\mbf{V_m}$ consists of $[\mbf v_1, \mbf v_2, \cdots , \mbf v_m]$. 
%After re-arranging (\ref{eq:rat_kry2}),
We can project the $e^{\mbf{{A}}}$ onto the rational Krylov subspace as follows.
\begin{eqnarray}
	\label{eq:scale}
	e^{\mbf{A}h} \mbf{v}  \approx
	 \beta \mbf{  V}_m e^{h\frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1.
\end{eqnarray}
In the line \ref{alg1:h} of Algorithm \ref{algo:arnoldi},  
$${\mbf{H}} = \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}.$$  
Following the same procedure \cite{Zhuang14_DAC, Botchev12}, 
the posterior error approximation is derived as  
\begin{eqnarray}
	\label{eq:err_rational_krylov}
   r_m(h) 
	= 	\lVert  \beta
	{h}_{m+1,m}  \frac{\mbf I - \gamma \mbf A_m}{\gamma}
	\mbf {  v}_{m+1} 
	\mbf e^\mathsf{T}_m \mbf{  H}_m^{-1}  
	%{ \phi( h\mbf{\wtd H_m})} 
	e^{ h\mbf{H}} 
	\mbf e_1  \rVert .
\end{eqnarray}
Note that in practice, instead of 
computing $(\mbf I-\gamma \mbf{A})^{-1}$ directly,
%\begin{eqnarray}
$(\mbf C + \gamma \mbf G)^{-1}\mbf C $ is utilized. %\mbf V_m =  \mbf V_m \mbf H_m + h_{m+1,m}
%\mbf v_{m+1} \mbf e^T_m . 
%\end{eqnarray}
The corresponding Arnoldi process shares  
the same skeleton of Algorithm \ref{algo:arnoldi} 
%and Alg. \ref{algo:ckt_solver}
with input matrices 
%The algorithm for rational basis of Krylov subspace in our MATEX is by setting 
$$\mbf X_1 = (\mbf C+ \gamma \mbf G)$$ for the LU decomposition Eq. (\ref{eq:lu}),
%$[\mbf L , \mbf U] = \text{LU\_Decompose}( \mbf {C} + \gamma \mbf {G} )$
and $$\mbf X_2 = \mbf {C}.$$
The residual estimation is % of $\left(\mbf C \dot {\mbf x}(t) + \mbf G \mbf x(t) \right)$  leads to 
\begin{eqnarray}
	\label{eq:er_rational_krylov}
  r(m,h)  
	= 	\lVert     \beta
	{h}_{m+1,m}
	\frac{\mbf C  + \gamma \mbf G}{\gamma} 
	\mbf {  v}_{m+1} 
	\mbf e^\mathsf{T}_m \mbf{  {H}}_m^{-1}  
	%{ \phi( h\mbf{\wtd H_m})} 
	e^{ h\mbf{  H}_m} 
	\mbf e_1 \rVert .
\end{eqnarray}
%{\color{red} 
Then, we plug Eq. (\ref{eq:er_rational_krylov}) into the line \ref{alg1:conv1} of Algorithm \ref{algo:arnoldi}.
}
%, which is $r(j,h)=\wtd r(j,h)$.
% Same as the I-MATEX, it mitigates the overestimation of error bound.
%\vspace{-0.02in}

\subsection{Regularization-Free MEVP Computation}
When $\mbf C$ is a singular matrix, MEXP \cite{Weng12_TCAD} needs the regularization process
\cite{Chen12_TCAD} to remove the singularity of DAE in Eq. (\ref{eqn:dae}).
It is because MEXP needs factorize $\mbf C$ directly to form the input $\mbf X_1$ for Algorithm 1. 
This brings extra computational overhead when the case is large \cite{Chen12_TCAD}. 
%\begin{comment}
It is not necessary if we can obtain the generalized eigenvalues and corresponding eigenvectors for matrix pencil $(-\mbf G, \mbf C)$. 
Based on \cite{Wilk:79}, we derive the following lemma,
%\vspace{+0.1in}
\begin{lemma}
	Considering a homogeneous system 
	%\begin{equation}
	$$ \mbf C \mbf {  \dot{x}} = -\mbf G \mbf x.$$ 
	%\end{equation}
	$\mbf u$ and $\lambda$ are the eigenvector and eigenvalue 
	of matrix pencil $(-\mbf G, \mbf C)$, then $$\mbf x =  e^{t\lambda} \mbf u$$ is a solution of the system.
\end{lemma}
%\vspace{+0.1in}
{ 
\begin{proof}\footnote{We repeat the proof from \cite{Wilk:79} with some modifications for our formulation.} 
If $\lambda$ and $\mbf u$ are an eigenvalue and eigenvector of a generalized eigenvalue problem $$-\mbf G \mbf u = \lambda \mbf C \mbf u.$$
Then, $ \mbf x = e^{t\lambda} \mbf u$ is the solution of $\mbf C \mbf {\dot x} =  - \mbf G \mbf x$. 
\end{proof}}
\vspace{+0.07in}
Because we do not need to compute $\mbf C^{-1}$ explicitly during Krylov subspace generation, I-MATEX and R-MATEX are regularization-free. Instead, we factorize $\mbf G$ for invert Krylov subspace basis generation (I-MATEX), 
or $(\mbf C+\gamma\mbf G)$ for rational Krylov subspace basis (R-MATEX).\footnote{It is also applied to the later work of DR-MATEX in Sec. \ref{sec:dist_frame}.}}
Besides, their Hessenberg matrices Eq. (\ref{eq:hess}) are  invertible, 
which contain corresponding important generalized eigenvalues/eigenvectors from 
matrix pencil $(-\mbf G, \mbf C)$, and define the behavior of linear dynamic system in Eq. (\ref{eqn:dae}) of  interest. 
  
\subsection{Comparisons among Different Krylov Subspace Algorithms for MEVP Computation}
\label{sec:compare_krylov}
 
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 $$\frac{Re(\lambda_{min})}{Re(\lambda_{max})} = 4.7 \times 10^{6},$$ where $\lambda_{max} = -8.49\times 10^{10}$ and $\lambda_{min} = -3.98 \times 10^{17}$ are the  maximum and minimum eigenvalues of $\mbf A = -\mbf C^{-1} \mbf G$. Fig. \ref{fig:order_reduction} 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\times$ of the size to achieve around relative error $1\%$ compared to Invert and Rational Krylov subspace methods.
The relative error is 
$$\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v ||},$$ 
where $h=0.4ps$, $\gamma = 10^{-13}$. 
The matrix $\mbf A$ is a relatively small matrix and computed by MATLAB $expm$ function. The result of $e^{h\mbf A}\mbf v$ serves as the baseline for accuracy. 
The relative error is the real relative difference compared to the analytical solution $e^{h\mbf A}\mbf v$ of the ODE $$\frac{d\mbf x}{dt}= \mbf A \mbf x$$ with an initial vector $\mbf 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 \times 10^4$, respectively. The large error of FE is due to the instability issue of its low order explicit time integration scheme.
In Fig. \ref{fig:order_reduction}, 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 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 $\mbf C$ and $\mbf G$ to make the different stiffness value $4.7 \times 10^{10}$. Fig. \ref{fig:order_reduction_2system} 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.
\begin{figure}[t]
	\centering
	\includegraphics[trim = 0.7cm 0.0cm 0.0cm 0.0cm, clip, keepaspectratio, width=0.4\textwidth]{./figs/TCAD_M_ErrReduction_04p.eps}
	\caption{
		%\color{red}
	{The relative error vs. dimensional $m$ of different Krylov subspace methods. The relative error is $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v ||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. Note: The relative error is the  difference compared to analytical solution $e^{h\mbf A}\mbf v$ of the ODE $ \frac{d\mbf x}{dt}= \mbf A \mbf x$ with an initial vector $\mbf v$, 
	which is generated by MATLAB $rand$ function, and its entries are positive numbers in $(0,1]$.  			 
	}
}
\label{fig:order_reduction}
\end{figure} 
\begin{figure}[t]
\centering
\includegraphics[trim = 0.7cm 0.0cm 0.0cm 0.0cm, clip, keepaspectratio, width=0.4\textwidth]{./figs/TCAD_M_ErrReduction_04p_twoSystem.eps}
\caption{{{The relative error vs. dimension $m$ of different Krylov subspace methods. The relative error is $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. The rational Krylov subspace has very stable error reduction rate. The number in the bracket represents the stiffness value of the system.}}}
\label{fig:order_reduction_2system}
\end{figure} 	
	
Regarding the relative error distributions vs. time step $h$ and dimension $m$, 
Fig. \ref{fig:std}, Fig. \ref{fig:inv}, and  Fig. \ref{fig:rat} are computed by standard, invert, and rational Krylov subspaces ($\gamma = 5\times 10^{-13}$), respectively.
Fig. \ref{fig:std} 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. \ref{fig:std}, invert (Fig. \ref{fig:inv}) and rational (Fig. \ref{fig:rat}) Krylov subspace  methods reduce errors quickly for large $h$.
%, which is useful to simulation. 
%It is because that the first smallest magnitude eigenvalues and their corresponding invariant subspaces play more dominating role, when we use larger step 
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 \cite{Van06}, which are efficiently captured by invert and rational Krylov subspace  methods.

\begin{figure}[t]
\centering
\includegraphics[trim = 0.1cm 0.0cm 0.0cm 0cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_std.eps}
\caption{The error  of MEVP via standard Krylov Subspace: 
$\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {  V_m} e^{h\mbf {  H}_m}||}{||e^{h\mbf A}\mbf v ||}$ 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 $\mbf 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. 
}\label{fig:std}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[trim = 0.15cm 0.0cm 0.0cm 1cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_inv.eps}
\caption{
{The error  of MEVP via invert Krylov Subspace:  $\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {V}_m e^{h\mbf {H }^{-1}_m}||}{||e^{h\mbf A}\mbf v ||}$ vs. time step $h$ and dimension of invert Krylov subspace basis ($m$). Compared to Fig. \ref{fig:std}, invert Krylov subspace method reduces the errors for large $h$.}
}
\label{fig:inv}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[trim = 0.15cm 0.0cm 0.0cm  1cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_rat.eps}
\caption{
{The error  of MEVP via rational Krylov Subspace: $\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $\gamma=5\times 10^{-13}$,  vs. time step $h$ and dimension of rational Krylov subspace basis ($m$). Compared to Fig. \ref{fig:std}, rational Krylov subspace method reduces the errors for large $h$ as Fig. \ref{fig:inv}.}
}
\label{fig:rat}
\end{figure}

The error of rational Krylov subspace is relatively insensitive to  $\gamma$ when it is selected between the time-step range of interests (Fig. \ref{fig:gamma}).
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. %Fig. \ref{fig:order_reduction} and Fig. \ref{fig:order_reduction_2system}. 

\begin{figure}[t]
	\centering
	\includegraphics[trim = 0cm 0cm 0cm 0cm, clip, keepaspectratio,  width=0.45\textwidth]{./figs/TCAD_err_sweep_r_04ps.eps}
	\caption{ The error  of MEVP via rational Krylov Subspace 
	{$\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=4ps$. The flat region shows the error is actually relatively insensitive to $\gamma$, when $\gamma$ is in the range of step size $h$ of interests. 
	}
}
\label{fig:gamma}
\end{figure}
 
	 


matex_method.tex
\subsection{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),
\begin{eqnarray}
    \label{eqn:dae} 
    \mbf{C}\dot{\mbf{x}}(t) = - \mbf{G}\mbf{x}(t) + \mbf{B}\mbf{u}(t),
\end{eqnarray}
where  $\mbf{C}$ is the  matrix for  capacitive and inductive elements. $\mbf{G}$ is the matrix for conductance and resistance, and $\mbf{B}$ is the input selector matrix.  $\mbf{x}(t)$ is the vector of time-varying node voltages
and branch currents. $\mbf{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 \cite{Nassif08_Power, Li12_Tau} to represent the activities under the networks. To solve Eq. (\ref{eqn:dae}) numerically, the system is discretized with time step $h$ and transformed to a linear algebraic system. Given an initial condition $\mbf x(0)$ from DC analysis or previous time step $\mbf x(t)$ and a time step $h$, $\mbf x(t+h)$ can be obtained by traditional \emph{low order approximation} methods \cite{Leon75}. 

\subsection{Traditional Low Order Time Integration Schemes}
\subsubsection{BE} Backward Euler based time integration scheme (Eq.(\ref{eqn:be})) is a robust implicit first-order method. 
\begin{eqnarray}
    \label{eqn:be}
    \Bigg(\frac{\mbf{C}}{h} + 
     {\mbf G} \Bigg) {\mbf{x}}(t+h) 
    =  \frac{\mbf C}{h} 
    \mbf{x}(t) + \mbf{B}     \mbf{u}(t+h)  .
  %  (2{\mbf{C}} + h{\mbf G}) {\mbf{x}}(t+h )  &   = & (2{\mbf C}- h{\mbf G} ) \mbf{x}(t) 
  %  \\ \nonumber
  % &+ & h\mbf{B}(\mbf{u}(t) +\mbf{u}(t+h))
\end{eqnarray}
\subsubsection{TR} Trapezoidal based time integration scheme (Eq.(\ref{eqn:trap})) is a popular implicit second-order method.
\begin{eqnarray}
  \label{eqn:trap}
    \Bigg(\frac{\mbf{C}}{h} +  
    \frac{\mbf G}{2}\Bigg) {\mbf{x}}(t+h)    & = &  
    \Bigg(\frac{\mbf C}{h}- \frac{\mbf G}{2} \Bigg)
    \mbf{x}(t)  \\ &+& \mbf{B} \frac{  \mbf{u}(t) + \mbf{u}(t+h)} {2}. \nonumber
\end{eqnarray}
It is probably the most commonly used strategy for large-scale circuit simulation, which has higher accuracy than BE.

\subsubsection{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 \cite{Xiong12,Yu12,Yang12,Li12_Tau}. 
% In this paper, we can them as BE-FTS and TR-FTS. 
If only one $h$ is used for the entire simulation, the choice is limited by the minimum breakpoint \cite{Nagel75} distance $h_{upper}$ among all the input sources. 
Fig. \ref{fig:interleave} (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. \ref{fig:interleave} (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. % such as Eq. (\ref{eqn:be}) and Eq. (\ref{eqn:trap}). 

\begin{comment}
\begin{figure}[ht] 
    \centering
    \includegraphics[width=3in]{./figs/interleave.eps}
     \caption{Example: interleaves two input source to create smaller transition time. ({\bf a}) Before interleaving, the input sources has smallest transition time $h_{upper}=10ps$; {\bf (b)} After interleaving, the input sources has the smallest transition time $h_{upper}=5ps$.}
    \label{fig:interleave1}
\end{figure}
\end{comment}

\begin{figure}[ht]
% \vspace{-0.1in}
	\centering
	\includegraphics[trim = 3.5cm 4cm 5cm 0cm, clip,keepaspectratio, width=0.45\textwidth]{./figs/interleaving_new.eps}
	%	\vspace{-0.2in}
	\caption{% \small 
		Example: Interleave  two input sources to create smaller transition time. ({  a}) Before interleaving, the smallest transition time of the input sources is $h_{upper}=10ps$; { (b)} After interleaving, the smallest transition time of the input sources is $h_{upper}=5ps$. %\cite{Zhuang15_TCAD}
	}
	\label{fig:interleave}	
%\vspace{-0.15in}
\end{figure}

%, which is very computational demanding.
%which is a very typical and 
%important application of linear circuit simulation. 
%power grid formulation is the same as Eq. (\ref{eqn:dae}).
%TR is the integration kernel, often applied to %the important linear circuit transient
%time domain simulation 
\subsection{Matrix Exponential Time Integration Scheme}
%\subsection{Exponential Time Differencing Method}

The solution of Eq. (\ref{eqn:dae}) can be obtained analytically \cite{Leon75}. 
For a simple illustration, we convert Eq. (\ref{eqn:dae}) into 
%an system of ordinary differential equations (ODEs),
\begin{eqnarray} 
\label{eqn:ode}
    \dot{\mbf{x}}(t) = \mbf{A}\mbf{x}(t) + \mbf{b}(t),
\end{eqnarray}
when $\mbf{C}$ is not singular\footnote{The assumption is to simplify the explanation in this section. After Sec. \ref{sec:imatex}, we use I-MATEX, R-MATEX and DR-MATEX to compute the solution of DAE without inversion of $\mbf C$. Therefore, the methods are suitable for general DAE system, i.e., Eq. (\ref{eqn:dae}) without the assumption here.}, 
$$\mbf{A} = - \mbf{C}^{-1}\mbf{G}~\text{, and}~\mbf{b}(t) = \mbf{C}^{-1}\mbf{B}\mbf{u}(t).$$
%This is the fundamental framework for our work in this paper,
Given  a solution at time 
$t$ and a time step $h$, the solution at
$t + h$ is
\begin{eqnarray}  
    \label{eqn:discrete_sol}
    \mbf{x}(t+h) = e^{h\mbf{A}}\mbf{x}(t) + 
    \int^{h}_{0}e^{(h-\tau)\mbf{A}}\mbf{b}(t+\tau)d\tau.
\end{eqnarray}
%Eq. (\ref{eqn:discrete_sol}) is the fundamental integration
%kernel of our framework MATEX.

Assuming that the input $\mbf{u}(t)$ is a piecewise linear (PWL) function of $t$, 
%e.g., $\mbf u(t)$ is linear within every time step, 
we can integrate the last term of Eq. (\ref{eqn:discrete_sol})  analytically, turning the solution with matrix exponential operator:
\begin{eqnarray} 
\label{eqn:pwl_exact_sol}
    \mbf{x}(t+h)  \nonumber = 
    % &=& 
      -\left(\mbf A^{-1}\mbf b(t+h) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h}\right) +
        \\
     %&  &
       e^{h\mbf{A}}\left(\mbf{x}(t) + \mbf A^{-1}\mbf b(t) + \mbf A^{-2} \frac{\mbf b(t+h)- \mbf b(t)}{h}\right )  .
  %  e^{h\mbf{A}}\mbf y(t)
   	%	\\
	%	\nonumber
\end{eqnarray}
\begin{comment}
Input transition spot $TS$: starting from a given
time point $t$ and $h \in [t,~t+TS]$, every points the input sources 
vector $\mbf u(t+h)$ maintains a \begin{const}
content
\end{environment-name}ant
slope. 
\end{comment}
For the time step choice,   breakpoints  (also known as input transition spots (TS) \cite{Zhuang14_DAC}) are the time points where slopes of input  vector change. 
Therefore, for Eq. (\ref{eqn:pwl_exact_sol}), the maximum time step starting from $t$ is $(t_s-t)$,
where $t_s$ is the smallest one in $TS$ 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.
% e.g. $\frac{ \mbf b(TS) - \mbf b(TS-\epsilon)}{\epsilon} \neq \frac{\mbf b(TS+\epsilon) - \mbf b(TS)}{\epsilon}$ when $\epsilon >0$. 


%which costs around $O(n^3)$ for $e^{\mbf A}$,
%where $n$ is the dimension of the matrix. 
%number of unknowns.



sub-concl.tex

In this work, we propose an efficient framework MATEX for accurate PDN time-domain simulation based on the  exponential 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 time-step scheme. Compared to the commonly adopted framework TR with fixed time step (TR-FTS), our single mode framework (R-MATEX) gains runtime speedup {\color{black} up to around $15\times$}. 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\times$} speedup.
 % by the distributed computing framework.  

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 \cite{Leon75}, 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 \cite{Hochbruck10_EXP, Hochbruck98_JSC} 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.  








































TCAD_2015_MATEX_final.tex
%\documentclass[journal,twoside,10pt]{IEEEtran}
%  \documentclass[a4paper,9pt]{sig-alternate}
\documentclass[10pt]{IEEEtran}
  
\newcommand{\mbf}{\mathbf}
\newcommand{\wtd}{\widetilde}

\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage{url}
\usepackage{multirow}
%\usepackage{cite}

\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts} 
\usepackage{amsthm}
\usepackage{color}
\usepackage{epstopdf}
\usepackage{verbatim} 
\usepackage[noadjust]{cite}
\usepackage{threeparttable}

%\usepackage{caption}
%\usepackage{subcaption}
%\newtheorem{lemma}{Lemma}
\usepackage{soul}
%\usepackage{amsthm}
\usepackage{verbatim, url}
\newtheorem{mydef}{Definition}
\newcounter{theorem}
\newcounter{def}
\setcounter{theorem}{0}
\setcounter{def}{0}
\newcounter{mylea}
\newcounter{corollary}
\setcounter{corollary}{0}
\setcounter{mylea}{0}
\usepackage{fancyhdr}
\newtheorem{definition}[def]{Definition}
\newtheorem{lemma}[mylea]{Lemma}
\newtheorem{mythe}[theorem]{Theorem}
\newtheorem{myclm}[theorem]{Claim}
\newtheorem{mycol}[corollary]{Corollary}
%\pagestyle{plain}
%% \setcounter{page}{1}
%\pagenumbering{arabic}
\renewcommand{\headrulewidth}{0pt}
\fancyhead[L]{}
\fancyhead[R]{}
\fancyfoot[C]{\thepage}
\pagestyle{plain}
\pagenumbering{arabic}

%\addtolength{\textwidth}{0.05in}

\begin{document}



%\newpage 

%\title{MATEX: Transient Simulation of Power Distribution Networks Using Matrix Exponential Operators and Krylov Subspace Algorithms}
\title{{Simulation Algorithms with Exponential Integration for Time-Domain Analysis of Large-Scale Power Delivery Networks}}
\author{
Hao~Zhuang, \IEEEmembership{Student Member, IEEE,} 
Wenjian~Yu, \IEEEmembership{Senior Member, IEEE,} 
Shih-Hung~Weng, 
Ilgweon~Kang, \IEEEmembership{Student Member, IEEE,}
Jeng-Hau~Lin, \IEEEmembership{Student Member, IEEE,}   
Xiang Zhang, \IEEEmembership{Member, IEEE,} 
Ryan Coutts,  
%Jingwei Lu, \IEEEmembership{Member, IEEE,} 
and~Chung-Kuan~Cheng, \IEEEmembership{Fellow, IEEE}\\  %$Preprint~of~IEEE~Transactions~on~Computer~Aided~Design~2016$  \\
\thanks{Manuscript received in May 2015; Revised in August, 2015 and November 2015. 
Accepted in Jan, 2016. } 
\thanks{ This work was supported in part by NSF-CCF 1017864, NSFC 61422402 and NSFC 61274033. 
%This work was supported by the National Science Foundation (NSF) under Grant CCF-1017864 and also partially supported by NSFC (grant No. 61422402).
%Corresponding author: H. Zhuang (hao.zhuang@cs.ucsd.edu)
 %and Chung-Kuan Cheng (ckcheng@ucsd.edu).
}
\thanks{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)} 
\thanks{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) } 
%\thanks{Shih-Hung Weng was with the Department of Computer Science and Engineering, University of California, San Diego, CA. 
%He is now a Research Scientist with Facebook Inc, Menlo Park, CA. shweng@fb.com
%s2weng@ucsd.edu}
\thanks{X. Zhang and R. Coutts are with Department of Electrical and Computer Engineering, University of California, San Diego, CA. (e-mail:xiz110@eng.ucsd.edu,~rmcoutts@eng.ucsd.edu) }
\thanks{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: %s2weng@eng.ucsd.edu and
shweng@fb.com)
}
 

%\thanks{Cadence Design Systems, Austin, TX 78759, gxin@cadence.com}
 
%\thanks{
%with the Department of Computer Science and Engineering, 
%University of California, San Diego}
%\thanks{Manuscript received ***, 2013; Corresponding author: Hao Zhuang (email: zhuangh@ucsd.edu).}
}
% make the title area
 
%\IEEEaftertitletext{\vspace{-0.0\baselineskip}}

%\input{letter}
%\newpage
%\input{reviews}
%\input{reviews2}
%\newpage 

% \markboth{IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems,~Vol.~XX, No.~X,~XXX~2015}{\MakeLowercase{\text Zhuang {et al.}}: Simulation Algorithms with Matrix Exponential Integration for Time-Domain Analysis of Power Delivery Networks}
%\onecolumn \input{copyright}
 

\twocolumn  
\maketitle


%\onecolumn
%\vspace{-0.15in}
% As a general rule, do not put math, special symbols or citations
% in the abstract or keywords.

\begin{abstract}
%\input{new_abst.tex}
\input{abst_tcad.tex}
\end{abstract}

\begin{comment}
{
%\scriptsize
%\vspace{-0.12in}
\section*{Categories and Subject Descriptors} 
%\vspace{-0.02in}
B.7.2 [Integrated Circuits]: Design Aids, Simulation. 
J.6 [Computer-aided engineering]: Computer-aided design (CAD).
% Note that keywords are not normally used for peerreview papers.

%\vspace{-0.12in}
\section*{General Terms} 
%\vspace{-0.02in}
Algorithms,  Design, Theory, Verification, Performance
%Circuit Simulation Algorithm, Matrix Exponential, Numerical integration, 
%Linear Circuits, Transient Simulation, 
%Adaptive Time Stepping, Approximation Algorithms,
%Power Distribution Networks, Power Grid
\end{comment}
%\vspace{-0.12in}
%\index Term}
%\vspace{-0.02in}
\begin{IEEEkeywords}
Circuit simulation, power delivery/distribution networks, power grid, time-domain simulation, transient simulation, matrix exponential, Krylov subspace method, parallel processing. 
\end{IEEEkeywords} \ 
\IEEEpeerreviewmaketitle

%\vspace{-0.12in}
%Introduction
\section{Introduction}
\label{sec:intro} 
%import from intro.tex
%\input{intro/intro}
\input{intro_ljw}
%\input{new_intro}%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Preliminary
%\section{Preliminary}
%\label{sec:preliminary}

%\subsection{Power Grid Formulation}
%\label{subsec:pwr_formulation}

%import from sub-pwr_grid_formulation.tex
%\input{sub-pwr_grid_formulation}

%\section{Circuit Formulation and Matrix Exponential Method}
% \vspace{-0.15in}
\section{Background}
\label{sec:method}
%import from sub-mexp_method.tex
\input{matex_method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Rational Krylov Space 

%\vspace{-0.08in}
%\section{MATEX with Modified Krylov Subspace Algorithms}
\section{Krylov Subspace Methods for Matrix Exponential and Vector Product (MEVP)}
% Krylov Subspace-Based Matrix Computation}
\label{sec:advanced}
 \input{krylov}

%\vspace{-0.15in}
\section{MATEX Framework}
%\underline{Ma}trix \underline{Ex}ponential Based Transient Simulation Framework}
%\section{MATEX with Standard Krylov Subspace Algorithm}
%\input{matex/krylov_app}
\label{sec:matex}
\input{framework}
 
 
%\subsection{DR-MATEX (Distributed R-MATEX Framework):  Accelerated Simulation by Decomposing the Groups of Current Sources and Embarrassingly Parallel Computation Model}
%Distributed/Parallel Processing.}
%\label{sec:dist_frame}
%\input{drmatex}
 
\section{Experimental Results}
\label{sec:exp_results} 
\input{exp_results} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Conclusion

\section{Conclusions and Future Directions} 
% and Future Work}
\label{sec:conclusion}
\input{sub-concl}
%\begin{comment}
%\vspace{-0.15in}

\section*{Acknowledgment}
\input{ack}
%\ifCLASSOPTIONcaptionsoff
%  \newpage
%\fi
%\bibliographystyle{abbrv}
\bibliographystyle{ieeetr}
%\vspace{-0.05in}
%\pagestyle{abbrv}
%\bibliography{../bib/CADRefs}  
\bibliography{CADJournal_ljw} 
% sigproc.bib is the name of the Bibliography in this case
 \input{bio} 
\end{document}



captions.tex
%\documentclass[journal,twoside,10pt]{IEEEtran}
%  \documentclass[a4paper,9pt]{sig-alternate}
\documentclass[10pt]{IEEEtran}
  
\newcommand{\mbf}{\mathbf}
\newcommand{\wtd}{\widetilde}

\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage{url}
\usepackage{multirow}
%\usepackage{cite}

\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts} 
\usepackage{amsthm}
\usepackage{color}
\usepackage{epstopdf}
\usepackage{verbatim} 
\usepackage[noadjust]{cite}
\usepackage{threeparttable}

%\usepackage{caption}
%\usepackage{subcaption}
%\newtheorem{lemma}{Lemma}
\usepackage{soul}
%\usepackage{amsthm}
\usepackage{verbatim, url}
\newtheorem{mydef}{Definition}
\newcounter{theorem}
\newcounter{def}
\setcounter{theorem}{0}
\setcounter{def}{0}
\newcounter{mylea}
\newcounter{corollary}
\setcounter{corollary}{0}
\setcounter{mylea}{0}
\usepackage{fancyhdr}
\newtheorem{definition}[def]{Definition}
\newtheorem{lemma}[mylea]{Lemma}
\newtheorem{mythe}[theorem]{Theorem}
\newtheorem{myclm}[theorem]{Claim}
\newtheorem{mycol}[corollary]{Corollary}
%\pagestyle{plain}
%% \setcounter{page}{1}
%\pagenumbering{arabic}
\renewcommand{\headrulewidth}{0pt}
\fancyhead[L]{}
\fancyhead[R]{}
\fancyfoot[C]{\thepage}
\pagestyle{plain}
\pagenumbering{arabic}

%\addtolength{\textwidth}{0.05in}

\begin{document}


%\newpage 

%\title{MATEX: Transient Simulation of Power Distribution Networks Using Matrix Exponential Operators and Krylov Subspace Algorithms}
%\title{{Simulation Algorithms with Exponential Integration for Time-Domain Analysis of Large-Scale Power Delivery Networks}}
%\author{} 
%\thanks{Cadence Design Systems, Austin, TX 78759, gxin@cadence.com}
 
%\thanks{
%with the Department of Computer Science and Engineering, 
%University of California, San Diego}
%\thanks{Manuscript received ***, 2013; Corresponding author: Hao Zhuang (email: zhuangh@ucsd.edu).}

% make the title area
 
%\IEEEaftertitletext{\vspace{-0.0\baselineskip}}

%\input{letter}
%\newpage
%\input{reviews}
%\input{reviews2}
%\newpage 

% \markboth{IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems,~Vol.~XX, No.~X,~XXX~2015}{\MakeLowercase{\text Zhuang {et al.}}: Simulation Algorithms with Matrix Exponential Integration for Time-Domain Analysis of Power Delivery Networks}
 
\onecolumn   
 
\noindent 
Copyright (c) 2015 IEEE. Personal use of this material is permitted. However, permission to use
this material for any other purposes must be obtained from the IEEE by sending an email to
pubs-permissions@ieee.org.
\newpage 
%\IEEEpee 

\noindent 
Fig. 1. % \small 
Example: Interleave two input sources to create smaller transition time. ({a}) Before interleaving, the smallest transition time of the input sources is $h_{upper}=10ps$; {(b)} After interleaving, the smallest transition time of the input sources is $h_{upper}=5ps$. %\cite{Zhuang15_TCAD}
\\

\noindent   
Fig. 2. 		{The relative error vs. dimensional $m$ of different Krylov subspace methods. The relative error is $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v ||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. Note: The relative error is the  difference compared to analytical solution $e^{h\mbf A}\mbf v$ of the ODE $ \frac{d\mbf x}{dt}= \mbf A \mbf x$ with an initial vector $\mbf v$, 
	which is generated by MATLAB $rand$ function, and its entries are positive numbers in $(0,1]$.  			 
} 			 
\\

\noindent   	
Fig. 3. {{{The relative error vs. dimension $m$ of different Krylov subspace methods. The relative error is $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. The rational Krylov subspace has very stable error reduction rate. The number in the bracket represents the stiffness value of the system.}}}
\\
 
\noindent  
Fig. 4. {The error  of MEVP via standard Krylov Subspace: 
	$\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {  V_m} e^{h\mbf {  H}_m}||}{||e^{h\mbf A}\mbf v ||}$ 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 $\mbf 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. 
}
\\

\noindent  
Fig. 5.	
{The error  of MEVP via invert Krylov Subspace:  $\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {V}_m e^{h\mbf {H }^{-1}_m}||}{||e^{h\mbf A}\mbf v ||}$ 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$.}
%The error  of MEVP via invert Krylov Subspace:  $\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {V}_m e^{h\mbf {H }^{-1}_m}||}{||e^{h\mbf A}\mbf v ||}$ 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$.
 \\
 
 \noindent  
Fig. 6. 
{The error  of MEVP via rational Krylov Subspace: $\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $\gamma=5\times 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.}
%The error  of MEVP via rational Krylov Subspace: $\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $\gamma=5\times 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.
 \\
 
 \noindent  
Fig. 7. { The error  of MEVP via rational Krylov Subspace 
	{$\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=4ps$. The flat region shows the error is actually relatively insensitive to $\gamma$, when $\gamma$ is in the range of step size $h$ of interests. 
	}\\
 
 \noindent  
 Fig. 8. {Part of a PDN model with input sources from Fig. 9}
 \\
 
 \noindent  
Fig. 9. 
	{Illustration of input transitions. 
		$GTS$: Global Transition Spots; $LTS$: Local Transition Spots;
		{$Snapshots$: the crossing positions by dash lines and LTS \#$k$ without solid points.} 
	}
\\

\noindent  
Fig. 10.  {Grouping of ``Bump'' shape transitions for sub-task simulation.
	The matrix exponential based method can  utilize  adaptive stepping in each $LTS$ and reuse the Krylov subspace basis generated at the latest point in $LTS$. 
	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:
	$t_{delay}$: delay time;
	$t_{rise}$: rise time;
	$t_{width}$: pulse width; 
	$t_{fall}$: fall time; and
	$t_{period}$: period).
}
 \\
 
 \noindent  
Fig. 11. DR-MATEX: The distributed MATEX framework using R-MATEX circuit solver.
\\
 
\newpage 
 
 
\noindent \begin{table}[h]  
  	\caption{Specifications of IBM power grid benchmarks.}  
\end{table}  
 
\begin{table}[h] 
	\caption{Performance comparisons (single computing node): TR-FTS, I-MATEX, and R-MATEX. 
	{\bf DC(\lowercase{s})}: Runtime of DC analysis (seconds);
	{{\bf \lowercase{$m$}$_I$}: The maximum \lowercase{$m$} of Krylov subspace in I-MATEX.}
	{\bf Tran(\lowercase{s})}: Runtime of transient simulation after DC (seconds), excluding the matrix factorization runtime; 
	{\bf Total(\lowercase{s})}: Runtime of overall transient simulation (seconds); \
	{{\bf Df(\lowercase{u}V)}: Maximum and average voltage differences compared to provided solutions (uV); } 
	{ {\bf  \lowercase{$m$}$_R$}: The maximum \lowercase{$m$}  of Krylov subspace in R-MATEX}
	%((seconds)  
	{\bf SPDP\lowercase{$^{r}_{tr}$}}: Speedup of R-MATEX over TR-FTS with respect to {\bf Tran(\lowercase{s})}; 
	{\bf SPDP\lowercase{$^{r}_{i}$}}:  Speedup of R-MATEX over I-MATEX with respect to {\bf Tran(\lowercase{s})}. } 
	% {\bf SPDP$^{hr}_{tr}$}:  Speedup of HR-MATEX over TR-FTS.\\ 
\end{table}  
 
\begin{table}[h] 
	\caption {{The performance of DR-MATEX (Distributed R-MATEX)}. {\bf Group \#}: Group number of the testcases. 
		%			Divided by the first breakpoints. 
		This number represents the total number of simulation sub-tasks for the design;
		{\bf Tran(\lowercase{s})}: Runtime of transient simulation after DC (seconds); 
		{\bf Total(\lowercase{s})}: Runtime of overall transient simulation (seconds); 
		%vs. TR ($h=10ps$); 
		{\bf Max. Df.(V)}  and {\bf Avg. Df.(V)}: maximum and average differences compared to the solutions of all output nodes provided by IBM power grid benchmarks; 
		%by IBM Power Grid Benchmarks solutions  and our internal reference solutions for the modified cases;
		%	{\bf SPDP$_{i}$}: Speedup over R-MATEX in Table \ref{tab:perf_mexp}.
		{\bf SPDP\lowercase{$_{tr}$}}: Speedup over TR-FTS's {\bf Tran(\lowercase{s})} in Table II;
		{\bf SPDP\lowercase{$_r$}}: Speedup over R-MATEX's {\bf Tran(\lowercase{s})} in Table II;
		{\bf Peak \lowercase{$m$}}: the peak dimension used in DR-MATEX  for MEVP; 
		{\bf Mem. Ratio over TR-FTS}: The peak memory comparison between the maximum memory consumption of DR-MATEX  over TR-FTS in Table  II.
		%Spdp$^4$:$ {t_{1000}/t_{rmatex}}$ speedups of transient stepping of MATEX over TR;
		%Spdp$^5$: $ {t_{t\_total}/t_{r\_total}}$ total runtime speedups of MATEX over TR. 
	}
\end{table}  
  
\end{document}



TCAD_2015_MATEX_final_preprint.tex
%\documentclass[journal,twoside,10pt]{IEEEtran}
%  \documentclass[a4paper,9pt]{sig-alternate}
\documentclass[10pt]{IEEEtran}
  
\newcommand{\mbf}{\mathbf}
\newcommand{\wtd}{\widetilde}

\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage{url}
\usepackage{multirow}
%\usepackage{cite}

\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts} 
\usepackage{amsthm}
\usepackage{color}
\usepackage{epstopdf}
\usepackage{verbatim} 
\usepackage[noadjust]{cite}
\usepackage{threeparttable}

%\usepackage{caption}
%\usepackage{subcaption}
%\newtheorem{lemma}{Lemma}
\usepackage{soul}
%\usepackage{amsthm}
\usepackage{verbatim, url}
\newtheorem{mydef}{Definition}
\newcounter{theorem}
\newcounter{def}
\setcounter{theorem}{0}
\setcounter{def}{0}
\newcounter{mylea}
\newcounter{corollary}
\setcounter{corollary}{0}
\setcounter{mylea}{0}
\usepackage{fancyhdr}
\newtheorem{definition}[def]{Definition}
\newtheorem{lemma}[mylea]{Lemma}
\newtheorem{mythe}[theorem]{Theorem}
\newtheorem{myclm}[theorem]{Claim}
\newtheorem{mycol}[corollary]{Corollary}
%\pagestyle{plain}
%% \setcounter{page}{1}
%\pagenumbering{arabic}
\renewcommand{\headrulewidth}{0pt}
\fancyhead[L]{}
\fancyhead[R]{}
\fancyfoot[C]{\thepage}
\pagestyle{plain}
\pagenumbering{arabic}

%\addtolength{\textwidth}{0.05in}

\begin{document}



%\newpage 

%\title{MATEX: Transient Simulation of Power Distribution Networks Using Matrix Exponential Operators and Krylov Subspace Algorithms}
\title{{Simulation Algorithms with Exponential Integration for Time-Domain Analysis of Large-Scale Power Delivery Networks}}
\author{
Hao~Zhuang, \IEEEmembership{Student Member, IEEE,} 
Wenjian~Yu, \IEEEmembership{Senior Member, IEEE,} 
Shih-Hung~Weng, 
Ilgweon~Kang, \IEEEmembership{Student Member, IEEE,}
Jeng-Hau~Lin, \IEEEmembership{Student Member, IEEE,}   
Xiang Zhang, \IEEEmembership{Member, IEEE,} 
Ryan Coutts,  
%Jingwei Lu, \IEEEmembership{Member, IEEE,} 
and~Chung-Kuan~Cheng, \IEEEmembership{Fellow, IEEE}\\ \[Preprint~for~IEEE~Transactions~on~Computer~Aided~Design ~2016\] \\
\thanks{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. 
%This work was supported by the National Science Foundation (NSF) under Grant CCF-1017864 and also partially supported by NSFC (grant No. 61422402).
%Corresponding author: H. Zhuang (hao.zhuang@cs.ucsd.edu)
 %and Chung-Kuan Cheng (ckcheng@ucsd.edu).
}
\thanks{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)} 
\thanks{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) } 
%\thanks{Shih-Hung Weng was with the Department of Computer Science and Engineering, University of California, San Diego, CA. 
%He is now a Research Scientist with Facebook Inc, Menlo Park, CA. shweng@fb.com
%s2weng@ucsd.edu}
\thanks{X. Zhang and R. Coutts are with Department of Electrical and Computer Engineering, University of California, San Diego, CA. (e-mail:xiz110@eng.ucsd.edu,~rmcoutts@eng.ucsd.edu) }
\thanks{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: %s2weng@eng.ucsd.edu and
shweng@fb.com)
}
 

%\thanks{Cadence Design Systems, Austin, TX 78759, gxin@cadence.com}
 
%\thanks{
%with the Department of Computer Science and Engineering, 
%University of California, San Diego}
%\thanks{Manuscript received ***, 2013; Corresponding author: Hao Zhuang (email: zhuangh@ucsd.edu).}
}
% make the title area
 
%\IEEEaftertitletext{\vspace{-0.0\baselineskip}}

%\input{letter}
%\newpage
%\input{reviews}
%\input{reviews2}
%\newpage 

% \markboth{IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems,~Vol.~XX, No.~X,~XXX~2015}{\MakeLowercase{\text Zhuang {et al.}}: Simulation Algorithms with Matrix Exponential Integration for Time-Domain Analysis of Power Delivery Networks}
%\onecolumn \input{copyright}
 
\twocolumn  
\maketitle

%\onecolumn
%\vspace{-0.15in}
% As a general rule, do not put math, special symbols or citations
% in the abstract or keywords.
\begin{abstract}
%\input{new_abst.tex}
\input{abst_tcad.tex}
\end{abstract}

\begin{comment}
{
%\scriptsize
%\vspace{-0.12in}
\section*{Categories and Subject Descriptors} 
%\vspace{-0.02in}
B.7.2 [Integrated Circuits]: Design Aids, Simulation. 
J.6 [Computer-aided engineering]: Computer-aided design (CAD).
% Note that keywords are not normally used for peerreview papers.

%\vspace{-0.12in}
\section*{General Terms} 
%\vspace{-0.02in}
Algorithms,  Design, Theory, Verification, Performance
%Circuit Simulation Algorithm, Matrix Exponential, Numerical integration, 
%Linear Circuits, Transient Simulation, 
%Adaptive Time Stepping, Approximation Algorithms,
%Power Distribution Networks, Power Grid
\end{comment}
%\vspace{-0.12in}
%\index Term}
%\vspace{-0.02in}
\begin{IEEEkeywords}
Circuit simulation, power delivery/distribution networks, power grid, time-domain simulaion, transient simulation, matrix exponential, Krylov subspace method, parallel processing. 
\end{IEEEkeywords}  
\IEEEpeerreviewmaketitle

%\vspace{-0.12in}
%Introduction
\section{Introduction}
\label{sec:intro} 
%import from intro.tex
%\input{intro/intro}
\input{intro_ljw}
%\input{new_intro}%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Preliminary
%\section{Preliminary}
%\label{sec:preliminary}

%\subsection{Power Grid Formulation}
%\label{subsec:pwr_formulation}

%import from sub-pwr_grid_formulation.tex
%\input{sub-pwr_grid_formulation}

%\section{Circuit Formulation and Matrix Exponential Method}
% \vspace{-0.15in}
\section{Background}
\label{sec:method}
%import from sub-mexp_method.tex
\input{matex_method}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Rational Krylov Space 

%\vspace{-0.08in}
%\section{MATEX with Modified Krylov Subspace Algorithms}
\section{Krylov Subspace Methods for Matrix Exponential and Vector Product (MEVP)}
% Krylov Subspace-Based Matrix Computation}
\label{sec:advanced}
 \input{krylov}

%\vspace{-0.15in}
\section{MATEX Framework}
%\underline{Ma}trix \underline{Ex}ponential Based Transient Simulation Framework}
%\section{MATEX with Standard Krylov Subspace Algorithm}
%\input{matex/krylov_app}
\label{sec:matex}
\input{framework}
 
 
%\subsection{DR-MATEX (Distributed R-MATEX Framework):  Accelerated Simulation by Decomposing the Groups of Current Sources and Embarrassingly Parallel Computation Model}
%Distributed/Parallel Processing.}
%\label{sec:dist_frame}
%\input{drmatex}
 
\section{Experimental Results}
\label{sec:exp_results} 
\input{exp_results} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Conclusion

\section{Conclusions and Future Directions} 
% and Future Work}
\label{sec:conclusion}
\input{sub-concl}
%\begin{comment}
%\vspace{-0.15in}

\section*{Acknowledgment}
\input{ack}
%\ifCLASSOPTIONcaptionsoff
%  \newpage
%\fi
%\bibliographystyle{abbrv}
\bibliographystyle{ieeetr}
%\vspace{-0.05in}
%\pagestyle{abbrv}
%\bibliography{../bib/CADRefs}  
\bibliography{CADJournal_ljw} 
% sigproc.bib is the name of the Bibliography in this case
 \input{bio} 
\end{document}



softcopies.tex
%\documentclass[journal,twoside,10pt]{IEEEtran}
%  \documentclass[a4paper,9pt]{sig-alternate}
\documentclass[10pt]{IEEEtran}
  
\newcommand{\mbf}{\mathbf}
\newcommand{\wtd}{\widetilde}

\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage{url}
\usepackage{multirow}
%\usepackage{cite}

\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{amsfonts} 
\usepackage{amsthm}
\usepackage{color}
\usepackage{epstopdf}
\usepackage{verbatim} 
\usepackage[noadjust]{cite}
\usepackage{threeparttable}

%\usepackage{caption}
%\usepackage{subcaption}
%\newtheorem{lemma}{Lemma}
\usepackage{soul}
%\usepackage{amsthm}
\usepackage{verbatim, url}
\newtheorem{mydef}{Definition}
\newcounter{theorem}
\newcounter{def}
\setcounter{theorem}{0}
\setcounter{def}{0}
\newcounter{mylea}
\newcounter{corollary}
\setcounter{corollary}{0}
\setcounter{mylea}{0}
\usepackage{fancyhdr}
\newtheorem{definition}[def]{Definition}
\newtheorem{lemma}[mylea]{Lemma}
\newtheorem{mythe}[theorem]{Theorem}
\newtheorem{myclm}[theorem]{Claim}
\newtheorem{mycol}[corollary]{Corollary}
%\pagestyle{plain}
%% \setcounter{page}{1}
%\pagenumbering{arabic}
\renewcommand{\headrulewidth}{0pt}
\fancyhead[L]{}
\fancyhead[R]{}
\fancyfoot[C]{\thepage}
\pagestyle{plain}
\pagenumbering{arabic}

%\addtolength{\textwidth}{0.05in}

\begin{document}



%\newpage 

%\title{MATEX: Transient Simulation of Power Distribution Networks Using Matrix Exponential Operators and Krylov Subspace Algorithms}
%\title{{Simulation Algorithms with Exponential Integration for Time-Domain Analysis of Large-Scale Power Delivery Networks}}
%\author{} 
%\thanks{Cadence Design Systems, Austin, TX 78759, gxin@cadence.com}
 
%\thanks{
%with the Department of Computer Science and Engineering, 
%University of California, San Diego}
%\thanks{Manuscript received ***, 2013; Corresponding author: Hao Zhuang (email: zhuangh@ucsd.edu).}

% make the title area
 
%\IEEEaftertitletext{\vspace{-0.0\baselineskip}}

%\input{letter}
%\newpage
%\input{reviews}
%\input{reviews2}
%\newpage 

% \markboth{IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems,~Vol.~XX, No.~X,~XXX~2015}{\MakeLowercase{\text Zhuang {et al.}}: Simulation Algorithms with Matrix Exponential Integration for Time-Domain Analysis of Power Delivery Networks}
 
\onecolumn   
 

\noindent 
Copyright (c) 2015 IEEE. Personal use of this material is permitted. However, permission to use
this material for any other purposes must be obtained from the IEEE by sending an email to
pubs-permissions@ieee.org.
\newpage 
 
 \begin{figure}[ht] 
 	\centering
 	\includegraphics[trim = 3.5cm 4cm 5cm 0cm, clip,keepaspectratio, width=0.45\textwidth]{./figs/interleaving_new.pdf}
 %	\caption{Example: interleaves two input source to create smaller transition time. ({\bf a}) Before interleaving, the input sources has smallest transition time $h_{upper}=10ps$; {\bf (b)} After interleaving, the input sources has the smallest transition time $h_{upper}=5ps$.}
 	\label{fig:interleave}
 \end{figure}
 
 
  \newpage
 \begin{figure}[ht]
 	\centering
 \includegraphics[trim = 0.7cm 0.0cm 0.0cm 0.0cm, clip, keepaspectratio, width=0.4\textwidth]{./figs/TCAD_M_ErrReduction_04p_twoSystem.eps}
% 	\caption{ 	{The relative error vs. dimensional $m$ of different Krylov subspace method. The relative error  $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v ||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. Note: the relative error is the real relative difference compared to analytical solution $e^{h\mbf A}\mbf v$ of the ODE $ \frac{d\mbf x}{dt}= \mbf A \mbf x$ with an initial vector $\mbf v$,  			which is generated by MATLAB $rand$ function and its entries are positive numbers in $(0,1]$.  	} 	}
 	\label{fig:order_reduction}
 \end{figure} 
 
 
 \newpage
 \begin{figure}[ht]
 	\centering
\includegraphics[trim = 0.7cm 0.0cm 0.0cm 0.0cm, clip, keepaspectratio, width=0.4\textwidth]{./figs/TCAD_M_ErrReduction_04p_twoSystem.eps}
 	%\caption{{{The relative error vs. dimension $m$ of different Krylov subspace method. The relative error $\frac{||e^{h\mbf A}\mbf v - \beta \mbf V_m e^{h\mbf H_m}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. The rational Krylov subspace has very stable error reduction rate. The number in the bracket represents the stiffness of the system.}}}
 	\label{fig:order_reduction_2system}
 \end{figure} 	
 
 
 
 \newpage
 \begin{figure}[h]
 	\centering
\includegraphics[trim = 0.1cm 0.0cm 0.0cm 0cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_std.eps}% 	\caption{		{The error  of MEVP via Rational Krylov Subspace: $\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $\gamma=5\times 10^{-13}$,  vs. time step $h$ and dimension of Rational Krylov subspace basis ($m$). Compared to Fig. \ref{fig:std}, Rational Krylov subspace method reduces the errors for large $h$ as Fig. \ref{fig:inv}.}}
 	\label{fig:rat}
 \end{figure}
 
  \newpage
 \begin{figure}[ht]
 	\centering
\includegraphics[trim = 0.15cm 0.0cm 0.0cm 1cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_inv.eps}
 %	\caption{The error  of MEVP via Standard Krylov Subspace:  		$\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {  V_m} e^{h\mbf {  H}_m}||}{||e^{h\mbf A}\mbf v ||}$ vs. time step $h$ and dimension of standard Krylov subspace basis ($m$). The standard Krylov subspace approximate well in extremely small $h$, since it captures the important eigenvalues of $\mbf A$ at that region. However, the small $h$ is not useful to simulation. For large $h$, it costs large amount of $m$ to reduce the error. }
 	\label{fig:std}
 \end{figure}
 
 \newpage 
 \begin{figure}[ht]
 	\centering
\includegraphics[trim = 0.15cm 0.0cm 0.0cm  1cm, clip, keepaspectratio, width=0.45\textwidth]{./figs/TCAD_err_distribution_rat.eps}
% 	\caption{ 		{The error  of MEVP via Invert Krylov Subspace:  $\frac{||e^{h\mbf A}\mbf v -  \beta \mbf {V}_m e^{h\mbf {H }^{-1}_m}||}{||e^{h\mbf A}\mbf v ||}$ vs. time step $h$ and dimension of Invert Krylov subspace basis ($m$). Compared to Fig. \ref{fig:std}, Invert Krylov subspace method reduces the errors for large $h$.} 	}
 	\label{fig:inv}
 \end{figure}
 
 
 
  \newpage
 
 \begin{figure}[h]
 	\centering
 	\includegraphics[trim = 0cm 0cm 0cm 0cm, clip, keepaspectratio,  width=0.45\textwidth]{./figs/TCAD_err_sweep_r_04ps.eps}
% 	\caption{ The error  of MEVP via Rational Krylov Subspace  {$\frac{||e^{h\mbf A}\mbf v - \beta \mbf {  V}_m e^{h \frac{\mbf I - \mbf { {H}}_m^{-1}}{\gamma}}\mbf e_1||}{||e^{h\mbf A}\mbf v||}$, where $h=4ps$. The flat region shows the error actually insensitive what $\gamma$ chosen in rational Krylov subspace method when $\gamma$ is in the range of step size $h$ of interests. 	}
 	\label{fig:gamma}
 \end{figure}
 
 
 
 
 
 
 
  \newpage
  \begin{figure}[ht] 
  	\centering
  	\includegraphics[width=2.3in]{./figs/rc_ckt.eps}
  	% 	\caption{Part of a PDN model with input sources from Fig. \ref{fig:input}.}
  	\label{fig:rc_ckt}
  \end{figure}
 
 \newpage
 
 \begin{figure*}[ht]
 	\centering
	\includegraphics[width=3.4in]%{./figs/time_spots.eps}
	{./figs/time_spots_2.eps} 
 %	\caption{DR-MATEX: The distributed MATEX framework using R-MATEX circuit solver}
 	\label{fig:drmatex}
 \end{figure*} 
 \newpage
 \begin{figure}[ht] 
 	\centering
 	\includegraphics[width=3.4in]%{./figs/ckt_ts_decomp.eps}
 	{./figs/source_large.eps}
 	\label{fig:ckt_ts_decomp}
 \end{figure}

 \newpage
 \begin{figure}[ht] 
 	\centering
	\includegraphics[width=0.6\textwidth]%{./figs/hierachical_ckt_sim.eps}
	{./figs/hierachical_ckt_sim_new.eps}
  	\label{fig:input}
 \end{figure} 

 
 
 
 
 
 
 
 
 
 \newpage

\begin{table}[h]
	\centering
	\caption{}
	\begin{tabular}{|c|r|r|r|r|r|r|r|} 
		\hline
		Design    & \#R     & \#C      & \#L  & \#I    & \#V  & \#Nodes  \\ \hline
		ibmpg1t   & 41K   & 11K    & 277  & 11K  & 14K  & 54K  \\ \hline
		ibmpg2t   & 245K  & 37K    & 330  & 37K  & 330    & 165K  \\  \hline
		ibmpg3t   & 1.6M & 201K   & 955  & 201K  & 955    & 1.0M \\ \hline
		ibmpg4t   & 1.8M & 266K   & 962  & 266K  & 962    & 1.2M \\ \hline
		ibmpg5t   & 1.6M & 473K   & 277  & 473K  & 539K   & 2.1M \\ \hline
		%092148 
		ibmpg6t   & 2.4M & 761K   & 381  & 761K &836K  & 3.2M \\ \hline
	\end{tabular}
	\label{tab:spec_ibmpg}
\end{table} 



\newpage
\begin{table}[h]
	\centering
	\caption{%\color{red}
	}
	\begin{threeparttable}[b]
		\begin{tabular}{|c|r|r|r||r|r|r|r||r|r|r| r||r|r| 
			}
			\hline
			\multirow{2}{*}{Design} & 
			\multirow{2}{*}{DC(s)} &
			\multicolumn{2}{c||}{TR-FTS} & 
			\multicolumn{4}{c||}{I-MATEX} & 
			\multicolumn{4}{c||}{R-MATEX} & 
			\multicolumn{2}{c|}{Speedups} 
			%	}
			\\ 
			\cline{3-14} 
			&	
			& Tran(s)	&Total(s)	
			% &LU(s)
			& $m_I$	
			& Tran(s)
			& Total(s) 
			& Df(uV)
			& $m_R$	
			& Tran(s)
			& Total(s) 
			& Df(uV)
			% s I-M 
			%	&Tran(s)
			%	& Total(s)		
			& SPDP$^{r}_{tr}$ %\footnote{ over TR }%  (/TR.)   
			& SPDP$^{r}_{i}$ %\footnote{ over TR }% (/TR.)	
			%			& SPDP$^{r}_{tr}$ %\footnote{ over TR (adpt)}%  (/TR.) 
			% m I-M	
			%	&LU(s)	&Total(s) % s R-M	
			%&Spdp2\footnote{ over TR }
			%	(/TR.)
			%			&Total(s) % m R-M 
			%&SPDP\footnote{ over I-MATEX(mn) } %	(/I- (mn))
			\\ \hline  \hline 
			ibmpg1t		     
			& 0.2  & 5.7 	& 6.00 
			& 30  	& 28.8 	& 28.9	& $58 \backslash 9.8$  
			& 5 	& 10.1 & 10.3  & $45 \backslash 6.8 $ & $0.6\times$ 
			& $2.9\times$
			\\ \hline 
			ibmpg2t    
			& 0.8  & 40.0 & 41.9 	
			& 28 & 130.0 &  130.9	& $92  \backslash 10.5$
			& 5  & 35.6 & 37.4 & $45 \backslash 6.8$ 
			&  $1.1\times$	&  $ 3.7\times$ 
			%	&  $1.9\times$	 
			\\ \hline 
			ibmpg3t 	     
			&  16.4  &  263.2  	& 295.0  
			& 29 & 1102.5  & 1115.1 	& $95 \backslash 20.4$ 
			&  5 & 	275.5  & 301.0 	&  $95 \backslash 18.5$
			& 	$1.0\times$ 	&	$4.0 \times$		% & $ 1.6\times$  
			\\ \hline
			ibmpg4t    
			& 13.5 & 460.8 	& 501.9
			&  29 
			& 433.8	& 458.2 &  $101  \backslash 39.3$
			&  5& 200.5  & 239.1   & $99 \backslash 34.2$ 
			& $ 2.3\times$
			& 	 $ 2.2\times$ %	& $ 6.4\times$
			\\ \hline
			ibmpg5t   
			&  9.0 	&  476.6 	&  	498.0 
			&  30 &  1934.4	& 1944.5  	&  $29  \backslash 5.6$		
			&  5 & 383.1   & 401.9 &  $29  \backslash 4.4$    
			&  $ 1.2\times$
			& $ 5.0\times$	% &$ 2.1\times$
			\\ \hline 
			ibmpg6t    
			&  15.3 	&  716.0 	& 749.1  	
			& 25  &  2698.9	& 2713.7  	&   $39\backslash 8.6$
			& 5 & 773.5 & 800.5 	&   $33\backslash 5.6$	
			%			& 323.20 & 357.10
			&  $0.9\times$ 
			& $3.5\times$ %	&  $2.2\times$
			\\ 
			\hline
			% \hline 
			ibmpg1t\_new	   
			&  0.2	&  51.3 & 51.7
			& 30 &  27.2  & 27.4 &  $58  \backslash 9.8$	  
			& 5 & 11.7   & 12.1  & $53 \backslash 6.9 $	 
			& $4.4\times$ 	& $2.3\times$% & $13.4\times$	 
			\\ \hline 
			ibmpg2t\_new 	      
			&  0.9   &  431.4 	&  	433.5
			&  28 	& 114.9  & 115.7 	&  $49 \backslash 10.5$	   
			& 5 	& 43.3 	& 44.9  &  $33\backslash 5.6$	  
			& $10.0\times$ 		& $2.7\times$	 %	&  $20.5\times$
			\\ \hline 
			ibmpg3t\_new	        
			&  16.3 &  3716.5  	&  	3749.0
			&  29	&  1219.3 & 1232.6  	&  $95 \backslash 20.4$	
			&  5	& 481.7   & 508.2 & $95 \backslash 18.9$	    
			&  $7.7\times$ 
			&  $2.5\times$	 
			\\ \hline
			ibmpg4t\_new 	     
			& 18.3 	& 5044.6 & 5085.3 
			&  29   
			& 753.5 	& 776.4 	&  $101 \backslash 39.3$		 
			& 6	& 350.9 & 387.2 & $99 \backslash 34.2$		
			& $ { 14.4}\times$ & 	 $2.1 \times$ 
			\\ \hline
			ibmpg5t\_new 	     
			& 10.5  & 5065.9 & 5110.1
			& 30  	& 2494.0	&  2504.7
			& $ 30 \backslash 5.6 $ 
			& 5  & 746.2	&  766.4 & $ 30 \backslash 4.4 $   &  $6.8\times$			
			&  $3.3\times$	  %	&  $20.8\times$
			\\ \hline 
			ibmpg6t\_new 	       
			&  13.1	&  7015.3 	&  	7059.7
			& 25 & 3647.9 & 3663.1 &  $39 \backslash 8.6 $ 
			& 6	& 895.1	&  923.1 &  $33\backslash 7.3 $ 
			&  $7.8\times$  &	$4.1\times$	 %	&  $19.4\times$
			\\ 
			\hline	
			\hline
			Average   & --- & ---
			& --- & --- & --- & --- 
			& $65\backslash 15.7$ &  --- & ---
			&  --- & $57\backslash 12.8$ 
			& $5\times$
			&$3\times$ % & $13.4\times$
			\\
			\hline 		
		\end{tabular} 
	\end{threeparttable}
	\label{tab:perf_mexp}
\end{table}


\newpage

\begin{table}[h] 
	\caption{
		} 
	\centering
	\begin{tabular}{|c||r|r|r|r|r||r|r||c|c|}
		\hline
		\multirow{2}{*}{Design}&
		\multicolumn{5}{c||}{DR-MATEX}&
		\multicolumn{2}{c||}{Speedups}
		& Peak   &
		{Mem. Ratio} 
		\\
		\cline{2-8}
		& Group \#
		& Tran(s)
		& Total(s)
		& Max Df.(V)
		& Avg Df.(V)
		& SPDP$_{tr}$
		%	& SPDP$_{hr}$
		& SPDP$_{r}$
		& $m$
		& over TR-FTS
		%& Total(s)
		\\
		\hline \hline 
		ibmpg1t & 100  &1.4  &  1.9 
		&  5.3e-5   & 8.6e-6   
		& $4.0\times$ 
		&  $7.1\times$      & 6 & 1.9
		\\
		\hline
		ibmpg2t & 100 & 8.9
		& 11.4 & 
		4.6e-5   & 8.6e-6    
		&  $4.5\times$  
		%	&  $11.7\times$   
		&   $4.0\times$  
		& 7 & 1.9 
		\\
		\hline
		ibmpg3t & 100 & 91.7  & 129.9
		%  69.71 & 124.41    
		&  9.6e-5 & 19.7e-6   &  $2.9\times$ 
		%	&  $9.5\times$  
		&   $4.4\times$  & 6  &  1.5 
		\\
		\hline
		ibmpg4t & 15  & 52.3 & 112.2	&  9.9e-5  & 27.9e-6    &  $8.8\times$   &    $3.8\times$  &  8 & 1.4  
		\\
		\hline   
		ibmpg5t & 100  	& 148.4 & 178.9  		&  9.0e-5  &  1.1e-6     
		&    $3.2\times$ 
		%&   $6.4\times$  
		&   $2.6\times$   &  7 % 1.4% 6.3/4.6
		& 1.5
		\\
		\hline
		ibmpg6t  & 100    
		& 189.9  & 234.2	&  3.4e-5  &  7.2e-6   
		&  $3.8\times$    
		& $4.1\times$   & 7 
		& 1.5 % 8.5/6.1~3.613.3
		\\
		\hline  
		ibmpg1t\_new  &  {\color{black} 100} 
		&  2.4   & 2.8 
		&  5.3e-5   & 8.6e-6     
		&  $21.8\times$  
		% & 3.8$\times$   
		& $5.0\times$  & 6  & 1.9
		\\
		\hline
		ibmpg2t\_new &{\color{black} 100}  
		% & 2.27 & 5.59
		& 5.6 & 7.0
		%	& 2.27 & 4.39   
		& 4.6e-5 & 8.6e-6  
		& $61.6\times$   
		% & 9.3$\times$  
		& $6.2\times$ &  7 & 1.9
		\\
		\hline
		ibmpg3t\_new  &  {\color{black} 100} & 	103.0 & 140.9 
		%  69.71 & 124.41
		% 47.61 & 113.87 % 136.96 & 186.28
		&  9.8e-5 & 19.9e-6 
		&  $25.6\times$  
		% &  $ 6.4\times$  
		&  $3.3\times$  
		%  &  $4.2\times$  &  $3.7\times$   &  $183.6\times$  
		& 7  & 1.5 
		\\
		\hline
		ibmpg4t\_new & {\color{black} 15}   
		& 51.5 & 108.4
		%		& 40.11 & 82.01  
		&  9.9e-5  & 27.6e-6 
		& {$ 98.0\times$}  
		% & $2.5\times$ 
		&  $6.8\times$  & 8   & 1.4 %  7.5/5.7 
		\\
		\hline
		%ibmpg5t & 331.67  & 350.33  &  100 & 35.67 & 54.21 &  0.7E-4 &  1.1E-5 &  9.3X  & 6.5X  
		ibmpg5t\_new  &  {\color{black} 100}  
		& 185.6 & 227.8
		%			& 45.56  & 69.20 
		& 9.9e-5  &  2.2e-6     
		& $27.3\times$ 
		% & $5.3\times$ 
		& $4.0\times$ 
		& 7 &  1.5
		\\
		\hline
		ibmpg6t\_new  & {\color{black} 100}   
		&  274.8  
		&  317.7
		%		& 62.92 	& 98.57
		&  3.4e-5  &  7.1e-6
		&  $25.5\times$  %&  $5.7\times$
		&  $3.3\times$   &  7  & 1.5   
		\\
		\hline 	\hline
		{Average}
		& --- & --- & --- &  {7.1e-5} &  {12.3e-6}   
		& { $26\times$}
		% & $6.4\times$ 
		&  { $5 \times$} &  {6.7} 
		& 
		{1.6}
		\\
		\hline 
	\end{tabular}
	\label{tab:pg_rmatex}
\end{table}



%computing runtimes for the transient part exculding LU.
%For R-MATEX, we use the maximum values among this 
%to define the $t_{adpt}$ in Table. \ref{tab:pg_rmatex}.
%MATEX with grouping is parallelizable by simulating each group in 
%multicore or cluster system.  
%The average number of pairs of backward and forward substitutions for Krylov subspace generation is around 60 ($km$ in Eq. (\ref{eqn:tr_spd})),
%while TR($h=10ps$) has 1000 pairs ($N$ in Eq. \ref{eqn:tr_spd}) on each cases.
%The reduction of these substitutions brings large speedups in the transient computing (the part of $K(T_H+T_e)$ in Eq. (\ref{eqn:tr_spd}) also costs portions of computational efforts). 
%and drags the number of speedups. 
%With the huge reduction of simulation time on this part, the serial LU and DC parts in MATEX  expose dominating roles, 





\begin{comment}
\begin{figure}[ht]
\centering
\includegraphics[width=3.4in]{./figs/result_ibmpg6t.eps}
\caption{Result of IMBPG6T}
\label{fig:result_ibmpg6t}
\end{figure}
\end{comment}


%which can be further improved using advance methods.
%~\cite{nicslu},
%which is another important study area.
%There are other way to decompose the input transitions, we use one simple
%strategy based on ``hump'' shapes, and achieve such large speedups.  

%In order to demonstrate the idea, we simply device a naive grouping scheme, which will take identical input sources together as a group. 
\begin{comment}
We assign each group to one computing nodes. 
In the experiment, we notice the message passing between 
our workstations are extreme small.
Time anchors are generated as a wall-lock, 
that visible to each node, so there are not other
communication overhead.
For example, for the message receiving only takes 0.0007 seconds, which is negligible. 
Table~\ref{tab:perf_mexp_group} shows the performance of MATEX with grouping. 
In this table, ``Groups'' is the number of groups for each 
designs, ``Min. Total`` and ``Max. Total'' is the maximum runtime among different nodes, 
%and ``Ideal Spdp'' records the ideal speedup with parallelization. 
Design ibmpg4t has only 5 groups because the alignment is better. 
The performance of MEXP has a significant improvement after grouping. 
The average ideal speedup over TR is 6 times. 
MATEX shows another parallelism in the power grid simulation from the perspective of 
input sources. 
\end{comment}

\begin{comment}
\begin{table*}[ht]
\centering
\caption{Performance of R-MATEX w/ Grouping}
\begin{tabular}{|c||r |r|r|r|r|r|}
\hline
\multirow{2}{*}{Design} & \multirow{2}{*}{TR Trans(s)} & \multicolumn{5}{c|} {R-MATEX w/ grouping}		      \\ \cline{3-7} 
&			     & Groups	    & Max. Trans (s)	& Min. Trans (s)  &Spdp1  &Spdp2\\ 
\hline ibmpg1t 	    & 5.94		     & 25	   & 0.84		& 0.48		  & 7.1X  & 12.4X \\
\hline ibmpg2t 	    & 26.98                   & 25	   & 3.93		& 2.38		  & 6.7X  & 11.3X \\
\hline ibmpg3t 	    & 245.92                  & 25	   & 38.09     	&21.96		  & 6.5X & 11.2X \\
\hline ibmpg4t 	    & 329.36                  & 5	   & 33.55		&30.12		  & 9.8X  & 10.9X \\
\hline ibmpg5t 	    & 331.67                  & 25	   & 58.54 		& 33.59		  & 5.7X & 9.9X \\
\hline ibmpg6t 	    & 542.04                  & 25	   & 81.40 		&51.68		  & 6.7X & 28.0X \\
\hline
\end{tabular}
\label{tab:perf_mexp_group}
\end{table*}


\subsection{Complexity}

The performance of MATEX and TR are dominated by the number of backward/forward substitutions. This is because the time span 
is relatively long in a power grid, and therefore, the expensive LU fatorization is amortized during the simulation. In MATEX, 
the number of backward/forward substitutions at each time step is actually $m$, which is the number of vectors forming rational 
Krylov subspace. Since $m$ is small, e.g., $5\sim 6$ in IBM benchmark suite, and MATEX uses large step size than TR, MATEX has 
less numbers of substitution operations in general. Table~\ref{tab:subs} presents number of backward/forward substitions for 
each method. ``Sub.'' is the total number of substituion operations along the simulation. MATEX has less operations but close to
TR because the ill alignment of input sources leads to many tiny step sizes whereas MATEX with grouping has much less operations 
due to well aligned inputs. These results also validate Tables~\ref{tab:perf_mexp} and \ref{tab:perf_mexp_group} where MATEX would 
has less speedup due to ill alignment issue, and the grouping scheme could improve the performance. 

\begin{table}[h]
\centering
\caption{Number of backward/forward substitutions}
\begin{tabular} {|c|c|c|c|}
\hline
\multirow{2}{*}{Design} &TR	 &MATEX		  &MATEX w/ grouping \\ \cline{2-4} 
& \# Sub.	 & \# Sub.        &  \# Max Sub.   \\ \hline
ibmpg1t   		    & 1000	 &  830		  &  120	    \\ \hline 
ibmpg2t   		    & 1000  	 &  830   	  &  120 	    \\ \hline
ibmpg3t   		    & 1000  	 &  874   	  &  120 	    \\ \hline 
ibmpg4t   		    & 1000  	 &  349   	  &  143 	    \\ \hline  
ibmpg5t   		    & 1000  	 &  982   	  &  139 	    \\ \hline
ibmpg6t   		    & 1000  	 &  932   	  &  134 	    \\ 
\hline 
\end{tabular}
\label{tab:subs}
\end{table}
\end{comment}

\begin{comment}

\subsection{Scalable IBM Power Grid Suite with Input Grouping}

Following the previous subsection, we not the input is the limitation for 
adaptive steping in our MATEX framework. 
We utlize input grouping~\ref{sec:group}.
Ideal one input sources per group.

\subsubsection{2 groups}

\subsubsection{4 groups}

\subsubsection{8 groups}

\subsubsection{16 groups}

\subsubsection{32 groups}

\subsubsection{64 groups}




Secondly, we consider a single input to run ibmpg1t-ibmpg6t 
to further illustrate our method. We see there are huge speedups than using trapezoidal method.
\end{comment}
%\subsubsection{Full Solver}
%This is the ideal speedups when computing resources are available.

%\subsubsection{Model Order Reduction Preprocessing}
%In this part, we compare the MOR method and put them into trapezoidal method
%and our MATEX framework, in order to compare our method accurate.



\begin{comment}
\begin{table}[h]
\caption{Comparisons among MEXP, I-MATEX and R-MATEX with RC cases}
\centering
\begin{tabular}{|c||c|c|c|c|}
\hline
Method  &   $m_{p} $  & Err(\%)
%	&  Time(s)
& Speedup 
& Stiffness\\
%& {$\frac{Re{\lambda_{smallest}(\mbf A)}}{Re({\lambda_{largest}(A)}}$}  \\ 
\hline
MEXP	  	&   229  	
&  0.510
%	&  601.87  
& -- &
\multirow{3}{*}{$2.1 \times 10^{16}$ } 
\\
\cline{1-4}
I-MATEX	 	&14       &  0.004
%	& 0.23
& 2616X   &   \\
\cline{1-4}
R-MATEX	   & 12   &  0.004  
%	& 0.22   
&  2735X   & \\ 
\hline
\hline
MEXP	  	&   224  	&  0.004  
%&  134.32
& --  
& 	    \multirow{3}{*}{$2.1 \times 10^{12}$ } \\
\cline{1-4}
I-MATEX		& 14       &0.004   
%&   0.23
&  583X  &   \\
\cline{1-4}
R-MATEX	&  12   &  0.004  
%	& 0.20 
&  611X  & \\ 
\hline
\hline
MEXP	   	&   223  & 0.004  
%	& 52.87  
&  -- & 	    \multirow{3}{*}{$2.1 \times 10^{8}$ } 
\\
\cline{1-4}
I-MATEX  	& 14      & 0.004   
%	&  0.23 
& 229X  &   \\
\cline{1-4}
R-MATEX	  & 12   &  0.004 
%	& 		0.21  
&  252X  & \\ 
%\cline{1-6}
\hline
\end{tabular}
\label{tab:rc_cases}
\end{table}
\end{comment}


\begin{comment}
\begin{table*}[htp]
%\scriptsize
\centering
\caption{Performance comparison TR, I-MATEX and R-MATEX of single-node computing (sn)/multi-node simulation (mn) versions,
Spdp$^1$: Speedup of I-MATEX (mn) over TR ($h=10ps$); 
Spdp$^2$: Speedup of R-MATEX (mn)over I-MATEX (mn); 
Spdp$^3$: Speedup of R-MATEX (mn) over TR ($h=10ps$);
Spdp$^4$: Speedup of R-MATEX (mn) over TR (adpt).  
}
\begin{threeparttable}[b]
\begin{tabular}{|c|r|r|r||r||r|r || r| r|| r| r ||  r|r|r|r| }
\hline
\multirow{2}{*}{Design} & 
\multirow{2}{*}{DC(s)} &
\multicolumn{2}{c||}{TR ($h=10ps$)} & 
\multicolumn{1}{c||}{TR (adpt)} & 
\multicolumn{2}{c||}{I-MATEX(sn)} & 
\multicolumn{2}{c||}{I-MATEX(mn)} & 
\multicolumn{2}{c||}{R-MATEX(sn) } & 
\multicolumn{4}{c|}{R-MATEX(mn)} 
\\ 
\cline{3-15} 
&		&LU(s)	&Total(s)	
% &LU(s)	
&Total(s) % s I-M 
&LU(s)
&Total(s) 
&Total(s) 
&SPDP\footnote{ over TR }% (/TR.)	
% m I-M	
&LU(s)	&Total(s) % s R-M	
%&Spdp2\footnote{ over TR }
%	(/TR.)
&Total(s) % m R-M 
&SPDP\footnote{ over I-MATEX(mn) } %	(/I- (mn))
& SPDP\footnote{ over TR }%  (/TR.)   
& SPDP\footnote{ over TR (adpt)}%  (/TR.)   
\\ \hline  
ibmpg1t		     & 0.13     
&0.13   &  6.20 	
& 29.48
& 0.12  & 27.23
% & 0.12  
& 4.68 
& 1.3X
& 0.13 & 4.93	  
%	     & 1.3X
%& 0.13 
& 1.10
& 4.5X 
& 5.6X
& 26.8X
\\ \hline 
ibmpg2t 	     & 0.80  
& 0.83   & 28.61	
& 179.74 
& 0.76   &  124.44
% & 0.76 
& 20.19 
& 1.4X
& 0.85 
& 25.90  
%	   & 1.1X
%& 0.85 
& 5.56
& 3.6X 
& 5.1X
& 32.3X
\\ \hline 
ibmpg3t 	     & 13.83   & 12.72  & 272.47
& 2792.96 
& 12.87 & 1246.69  %1988.98
% & 12.58 
& 213.33 %280.66 
& 1.3X 
& 12.73 
& 244.56  
%  & 1.1X
% & 12.81 
& 63.78 
& 3.3X 
& 4.3X 
& 43.8X
\\ \hline
ibmpg4t 	    & 16.69    & 22.50   &	368.55
& 1773.14 
& 21.09 & 484.83 %625.12
% & 20.99 
& 148.45 %231.06 
& 2.5X
& 22.09  & 140.42    
%   & 2.6X
% & 22.02  
& 76.55  
& 1.9X  
& 4.8X
& 23.1X
\\ \hline
ibmpg5t 	    & 8.16    & 10.50  & 350.33  	
& 2359.11
& 9.92 & 1821.60 
% & 10.03 
& 304.72 %413.77
& 1.1X
& 10.77 &  337.74 
%   & 1.0X
% & 10.51  
&  78.99 
&    3.9X 
& 4.4X
& 29.9X
\\ \hline 
ibmpg6t 	     & 11.17    &14.17  & 567.38
& 3184.59
& 13.19 & 2784.46 %4059.30
% & 13.55 
& 479.79 %856.09
& 1.2X
& 14.18 & 482.42  
%	    & 1.2X
% & 14.15 
& 108.92 
& 4.4X 
&   5.2X 
& 29.2X 
\\ 
\hline
\end{tabular}
% \begin{tablenotes}\scriptsize
%	\item [1] Speedup of I-MATEX (mn) over TR ($h=10ps$) 
%	\item [2] Speedup of R-MATEX (mn)over I-MATEX (mn) 
%	\item [3] Speedup of R-MATEX (mn) over TR ($h=10ps$)
%	\item [4] Speedup of R-MATEX (mn) over TR (adpt)      
%	\vspace{-0.5in}
%    \end{tablenotes}
\end{threeparttable}
\label{tab:perf_mexp}
\end{table*}
\end{comment}

 
 
 
\end{document}



TCAD_2015_MATEX_final.bbl
\begin{thebibliography}{10}

\bibitem{Najm03_DAC}
D.~Kouroussis and F.~N. Najm, ``A static pattern-independent technique for
  power grid voltage integrity verification,'' in {\em Proc. IEEE/ACM Design
  Autom. Conf.}, pp.~99--104, 2003.

\bibitem{Nassif00_Fast}
S.~R. Nassif and J.~N. Kozhaya, ``Fast power grid simulation,'' in {\em Proc.
  IEEE/ACM Design Autom. Conf.}, pp.~156--161, 2000.

\bibitem{Brooks07_DATE}
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 {\em Proc. DATE}, pp.~1--6, 2007.

\bibitem{Lin04_CICC}
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 {\em Proc. IEEE
  CICC}, pp.~509--512, 2004.

\bibitem{Lin01_ICCAD}
S.~Lin and N.~Chang, ``Challenges in power-ground integrity,'' in {\em Proc.
  IEEE/ACM Int. Conf. Comput.-Aided Design}, pp.~651--654, 2001.

\bibitem{Skadron12}
R.~Zhang, B.~H. Meyer, W.~Huang, K.~Skadron, and M.~R. Stan, ``Some limits of
  power delivery in the multicore era,'' {\em WEED}, 2012.

\bibitem{Skadron14_ISCA}
R.~Zhang, K.~Wang, B.~H. Meyer, M.~R. Stan, and K.~Skadron, ``Architecture
  implications of pads as a scarce resource,'' in {\em Proc. IEEE/ACM Int’l
  Symp. Computer Architecture}, pp.~373--384, 2014.

\bibitem{Skadron14_ASPDAC}
K.~Wang, B.~H. Meyer, R.~Zhang, K.~Skadron, and M.~R. Stan, ``Walking pads:
  Fast power-supply pad-placement optimization.,'' in {\em Proc. IEEE/ACM Asia
  South Pac. Design Autom. Conf.}, vol.~20, p.~4, 2014.

\bibitem{ljw14_DAC}
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
  {\em Proc. IEEE/ACM Design Autom. Conf.}, 2014.

\bibitem{Pan05_ICCAD}
M.~Pan, N.~Viswanathan, and C.~Chu, ``An efficient and effective detailed
  placement algorithm,'' in {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided
  Design}, pp.~48--55, 2005.

\bibitem{ljw15_TCAD}
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, ``{ePlace-MS}: Electrostatics based
  placement for mixed-size circuits,'' {\em IEEE Trans. Comput.-Aided Design
  Integr. Circuits Syst.}, vol.~34, no.~5, pp.~685--698, 2015.

\bibitem{Xiao10_ICCAD}
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 {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided Design},
  pp.~458--462, 2010.

\bibitem{Zhang12_ICCAD}
Y.~Zhang and C.~Chu, ``{GDRouter}: Interleaved global routing and detailed
  routing for ultimate routability,'' in {\em Proc. IEEE/ACM Design Autom.
  Conf.}, pp.~597--602, 2012.

\bibitem{Kahng13_ICCAD}
A.~B. Kahng, S.~Kang, H.~Lee, I.~L. Markov, and P.~Thapar, ``High-performance
  gate sizing with a signoff timer,'' in {\em Proc. IEEE/ACM Int. Conf.
  Comput.-Aided Design}, pp.~450--457, 2013.

\bibitem{Zhuo12_ICCAD}
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 {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided Design},
  pp.~255--262, 2012.

\bibitem{Zeng10_DAC}
Z.~Zeng, X.~Ye, Z.~Feng, and P.~Li, ``Tradeoff analysis and optimization of
  power delivery networks with on-chip voltage regulation,'' in {\em Proc.
  IEEE/ACM Design Autom. Conf.}, pp.~831--836, 2010.

\bibitem{Zhuang13_ICCCAS}
H.~Zhuang, J.~Lu, K.~Samadi, Y.~Du, and C.~K. Cheng, ``Performance-driven
  placement for design of rotation and right arithmetic shifters in monolithic
  {3D ICs},'' in {\em Proc. IEEE Intl. Conf. Communications, Circuits and
  Systems}, vol.~2, pp.~509--513, 2013.

\bibitem{Lim14_ICCAD}
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 {\em
  Proc. IEEE/ACM Int. Conf. Comput.-Aided Design}, pp.~565--572, 2014.

\bibitem{xiang14}
D.~Xiang and K.~Shen, ``A thermal-driven test application scheme for pre-bond
  and post-bond testing of three-dimensional {IC}s,'' {\em ACM Journal on
  Emerging Technologies of Computing Systems}, vol.~10, no.~2, p.~Art. no. 18,
  2014.

\bibitem{Zhuang15_emc}
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 {\em IEEE Intl. Symp. EMC\&SI}, 2015.

\bibitem{Zhuang15_DAC}
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 {\em Proc. IEEE/ACM Design Autom. Conf.}, 2015.

\bibitem{Nassif08_Power}
S.~R. Nassif, ``Power grid analysis benchmarks,'' in {\em Proc. IEEE/ACM Asia
  South Pac. Design Autom. Conf.}, pp.~376--381, 2008.

\bibitem{zuochang08}
Z.~Ye, Z.~Zhu, and J.~R. Phillips, ``Sparse implicit projection ({SIP}) for
  reduction of general many-terminal networks,'' in {\em Proc. Intl. Conf.
  Computer-Aided Design}, pp.~736--741, 2008.

\bibitem{Li12_Tau}
Z.~Li, R.~Balasubramanian, F.~Liu, and S.~Nassif, ``2012 tau power grid
  simulation contest: benchmark suite and results,'' in {\em Proc. IEEE/ACM
  Int. Conf. Comput.-Aided Design}, pp.~643--646, 2012.

\bibitem{Zhuo14_DAC}
C.~Zhuo, H.~Gan, and W.-K. Shih, ``Early-stage power grid design: Extraction,
  modeling and optimization,'' in {\em Proc. IEEE/ACM Design Autom. Conf.},
  pp.~1--6, 2014.

\bibitem{Zhuang2012}
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 {G}reen's functions,'' in {\em Proc. IEEE/ACM Asia South
  Pac. Design Autom. Conf.}, pp.~377--382, 2012.

\bibitem{Yu14_TCAD}
C.~Zhang and W.~Yu, ``Efficient space management techniques for large-scale
  interconnect capacitance extraction with floating random walks,'' {\em IEEE
  Trans. Comput.-Aided Design Integr. Circuits Syst}, vol.~32, no.~10,
  pp.~1633--1637, 2013.

\bibitem{Yu2013}
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,'' {\em IEEE Trans. Comput.-Aided Design Integr. Circuits
  Syst}, vol.~32, no.~3, pp.~353--366, 2013.

\bibitem{Leon75}
L.~O. Chua and P.-M. Lin, {\em Computer Aided Analysis of Electric Circuits:
  Algorithms and Computational Techniques}.
\newblock Prentice-Hall, 1975.

\bibitem{Saad03}
Y.~Saad, {\em Iteravite Methods for Sparse Linear Systems}.
\newblock SIAM, 2003.

\bibitem{Davis06}
T.~A. Davis, {\em Direct Method for Sparse Linear Systems}.
\newblock SIAM, 2006.

\bibitem{Yu12}
T.~Yu and M.~D.-F. Wong, ``{PGT\_SOLVER}: {A}n efficient solver for power grid
  transient analysis,'' in {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided
  Design}, pp.~647--652, 2012.

\bibitem{Yang12}
J.~Yang, Z.~Li, Y.~Cai, and Q.~Zhou, ``Powerrush: Efficient transient
  simulation for power grid analysis,'' in {\em Proc. IEEE/ACM Int. Conf.
  Comput.-Aided Design}, pp.~653--659, 2012.

\bibitem{Xiong12}
X.~Xiong and J.~Wang, ``Parallel forward and back substitution for efficient
  power grid simulation,'' in {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided
  Design}, pp.~660--663, 2012.

\bibitem{Nagel75}
L.~Nagel, {\em {SPICE2: A computer program to simulate semiconductor
  circuits}}.
\newblock Ph.D. dissertation, 1975.

\bibitem{Weng12_TCAD}
S.-H. Weng, Q.~Chen, and C.~K. Cheng, ``Time-domain analysis of large-scale
  circuits by matrix exponential method with adaptive control,'' {\em IEEE
  Trans. Comput.-Aided Design Integr. Circuits Syst}, vol.~31, no.~8,
  pp.~1180--1193, 2012.

\bibitem{Saad92}
Y.~Saad, ``Analysis of some krylov subspace approximations to the matrix
  exponential operator,'' {\em SIAM J. Numer. Anal.}, vol.~29, no.~1,
  pp.~209--228, 1992.

\bibitem{cloud}
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,'' {\em Communications of the ACM}, vol.~53, no.~4, pp.~50--58,
  2010.

\bibitem{mapreduce}
J.~Dean and S.~Ghemawat, ``Mapreduce: simplified data processing on large
  clusters,'' {\em CACM}, vol.~51, no.~1, pp.~107--113, 2008.

\bibitem{xiang09}
D.~Xiang, Y.~Zhang, and Y.~Pan, ``Practical deadlock-free fault-tolerant
  routing in meshes based on the planar network fault model,'' {\em IEEE Trans.
  Computers}, vol.~58, no.~5, pp.~620--633, 2009.

\bibitem{spark10}
M.~Zaharia, M.~Chowdhury, M.~J. Franklin, S.~Shenker, and I.~Stoica, ``Spark:
  cluster computing with working sets,'' in {\em HotCloud}, pp.~10--10, 2010.

\bibitem{He14_emc}
Q.~He, W.~Au, A.~Korobkov, and S.~Venkateswaran, ``Parallel power grid analysis
  using distributed direct linear solver,'' in {\em IEEE Intl. Symp. EMC},
  pp.~866--871, 2014.

\bibitem{Jia15_TCAD}
N.~Gupte and J.~Wang, ``Secure power grid simulation on cloud,'' {\em IEEE
  Trans. Comput.-Aided Design Integr. Circuits Syst}, vol.~34, no.~3,
  pp.~422--432, 2015.

\bibitem{Wang13}
J.~Wang and X.~Xiong, ``Scalable power grid transient analysis via
  {MOR}-assisted time-domain simulations,'' in {\em Proc. IEEE/ACM Int. Conf.
  Comput.-Aided Design}, 2013.

\bibitem{Tan13_ICCAD}
X.-X. Liu, H.~Wang, and S.~X.-D. Tan, ``Parallel power grid analysis using
  preconditioned gmres solver on {CPU-GPU} platforms,'' in {\em Proc. IEEE/ACM
  Int. Conf. Comput.-Aided Design}, pp.~561--568, IEEE, 2013.

\bibitem{Feng08}
Z.~Feng and P.~Li, ``Multigrid on {GPU}: tackling power grid analysis on
  parallel simt platforms,'' in {\em Proc. IEEE/ACM Int. Conf. Comput.-Aided
  Design}, pp.~647--654, 2008.

\bibitem{Zhuang14_DAC}
H.~Zhuang, S.-H. Weng, J.-H. Lin, and C.~K. Cheng, ``{MATEX}: {A} distributed
  framework of transient simulation of power distribution networks,'' in {\em
  Proc. IEEE/ACM Design Autom. Conf.}, 2014.

\bibitem{gander2013paraexp}
M.~J. Gander and S.~Güttel, ``{PARAEXP}: A parallel integrator for linear
  initial-value problems,'' {\em SIAM Journal on Scientific Computing},
  vol.~35, no.~2, pp.~C123--C142, 2013.

\bibitem{Weng12_ICCAD}
S.-H. Weng, Q.~Chen, N.~Wong, and C.~K. Cheng, ``Circuit simulation via matrix
  exponential method for stiffness handling and parallel processing,'' in {\em
  Proc. IEEE/ACM Int. Conf. Comput.-Aided Design}, pp.~407--414, 2012.

\bibitem{Van06}
J.~van~den Eshof and M.~Hochbruck, ``Preconditioning {Lanczos} approximations
  to the matrix exponential,'' {\em SIAM J. Sci. Comput.}, vol.~27, no.~4,
  pp.~1438--1457, 2006.

\bibitem{Botchev12}
M.~A. Botchev, ``A short guide to exponential krylov subspace time integration
  for maxwell's equations,'' Dept. of Applied Mathematics, Univ. of Twente,
  2012.

\bibitem{Zhuang13_ASICON}
H.~Zhuang, S.-H. Weng, and C.~K. Cheng, ``Power grid simulation using matrix
  exponential method with rational krylov subspaces,'' in {\em Proc. IEEE Intl.
  Conf. ASIC}, 2013.

\bibitem{Chen12_TCAD}
Q.~Chen, S.-H. Weng, and C.~K. Cheng, ``A practical regularization technique
  for modified nodal analysis in large-scale time-domain circuit simulation,''
  {\em IEEE Trans. Comput.-Aided Design Integr. Circuits Syst}, vol.~31, no.~7,
  pp.~1031--1040, 2012.

\bibitem{Wilk:79}
J.~Wilkinson, ``Kronecker's canonical form and the {QZ} algorithm,'' {\em
  Linear Algebra and its Applications}, vol.~28, pp.~285--303, 1979.

\bibitem{Li11}
P.~Li, ``Parallel circuit simulation: A historical perspective and recent
  developments,'' {\em Foundations and Trends in Electronic Design Automation},
  vol.~5, no.~4, pp.~211--318, 2012.

\bibitem{Ye08_ISQED}
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 {\em ISQED}, pp.~627--632, 2008.

\bibitem{mesos}
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 {\em NSDI}, pp.~22--22, 2011.

\bibitem{nicslu}
X.~Chen, Y.~Wang, and H.~Yang, ``{NICSLU}: An adaptive sparse matrix solver for
  parallel circuit simulation,'' {\em IEEE Trans. Comput.-Aided Design Integr.
  Circuits Syst}, vol.~32, no.~2, pp.~261--274, 2013.

\bibitem{Hochbruck10_EXP}
M.~Hochbruck and A.~Ostermann, ``Exponential integrators,'' {\em Acta
  Numerica}, vol.~19, pp.~209--286, 2010.

\bibitem{Hochbruck98_JSC}
M.~Hochbruck, C.~Lubich, and H.~Selhofer, ``Exponential integrators for large
  systems of differential equations,'' {\em SIAM J. Sci. Comput.}, vol.~19,
  no.~5, pp.~1552--1574, 1998.

\end{thebibliography}
