text:
\section{Appendix}

\subsubsection{Motivation}
With the ability of large step time-marching, the input term $\mathbf 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
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 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
approximation of $e^{\mathbf A}\mathbf v$
is accurate enough. The Krylov subspace can be reused
when $t+h \in [t,t_s]$, only by scaling $\mathbf H_m$ with $h$
in
$$\mathbf x(t+h) = \lVert \mathbf v \rVert \mathbf {V_m} e^{h \mathbf H_m} \mathbf e_1-\mathbf 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
\emph{Local Transition Spot} ($LTS$): the set of $TS$ at an input source to the PDN.
\item
\emph{{Global Transition Spot}} ($GTS$): the union of $LTS$ among all the input sources to the PDN.
\item
\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_new.eps}
\centering{
\caption{DR-MATEX: the distributed MATEX framework using R-MATEX circuit solver}
}
\label{fig:drmatex}
}
\end{figure*}

\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]
    \centering
    \includegraphics[width=3.7in
    {./figs/time_spots_2.eps}

    \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.}

	}
	}
    }
    \label{fig:input}
\end{figure}

\begin{figure}[ht]
    \centering
    \includegraphics[width=3.65in
    {./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 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$, $\mathbf X_1$, $\mathbf X_2$, and $\mathbf P_k$, $\mathbf F_k$, which contain the corresponding $\mathbf b$ for node $k$; error
tolerance $E_{tol}$, and simulation time $T$.}
\KwOut{Local solution $\mathbf x$ along $GTS$ in node
$k \in [1,\cdots,S]$, where $S$ is the number of nodes}
    {
	$t=T_{start}$\;

	$\mathbf x(t) = \text{Local\_Initial\_Solution}$\;

	$[\mathbf L,\mathbf U] = \text{LU\_Decompose}( \mathbf X_1)$\;
	\While { $t \leq T_{end}$}
	{
	Compute maximum allowed step size $h$ based on $GTS$\;

	 {

	\If{$t \in $  LTS \#$k$}
	{
	\tcc{Generate the Krylov subspace for the time point of $LTS$ and compute $\mathbf x$}
	$[\mathbf x(t+h), \mathbf V_m, \mathbf H_m,\mathbf v] = \text{MATEX\_Arnoldi}(\mathbf L,\mathbf U,\mathbf X_2,$ $ h, t, \mathbf x(t), \epsilon, \mathbf P_k(t,h), \mathbf F_k(t,h))$\;

	$a_{lts} = t$\;

	}
	\Else{
	\tcc{Compute the $\mathbf x$ at $Snapshot$ by reusing the latest Krylov subspace}
	    {

	    $h_{a} = t + h - a_{lts}$\;
	    $\mathbf x(t+h) = \lVert \mathbf v \rVert \mathbf {V_m} e^{h_{a} \mathbf H_m} \mathbf e_1 - \mathbf P_k(t,h)$\;

	    }
	}
	$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 $\mathbf H_m$ is $T_{H}$ which costs time complexity $O(m^3)$, plus extra $T_e$ to form $\mathbf x$, which costs $O(n m^2)$.
Given $K$ points of $GTS$, without decomposition of input transitions,
the time complexity is
$KmT_{bs} + K (T_{H}+ T_e)  + T_{serial}$.
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 complexity is
$kmT_{bs} + K (T_{H}+ T_e)  + T_{serial}$, where  $K (T_{H}+ T_e)$ contains the portion of computing $Snapshot$.
The speedup of distributed computing over single MATEX,
\begin{eqnarray}
\text{Speedup}=\frac{ KmT_{bs} + K(T_{H}+ T_e)+ T_{serial}}
{ kmT_{bs}+ K(T_{H}+ T_e)+T_{serial}}.
\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{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}
\subsubsection{Simple $\mathbf C$ matrix cases}

\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
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
				& 0.13    & 33.67s
				& 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
				& 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)
			     &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
			   & 2.2X
			   & 244.56
			   & 11.4X
			   & 5.1X
			   \\ \hline
    	ibmpg4t 	    & 16.69
			    & 1773.14
			    & 484.83
			   & 3.7X
			    & 140.42
			    & 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
			   & 1.1X
			    & 482.42
			    &   6.6X
			   &5.8X
			    \\
    \hline
    \end{tabular}
    \label{tab:perf_mexp}

\end{table*}
\end{comment}

\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

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.
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 $\mathbf 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.

\begin{table*}[t]
	\centering
	 	\caption
		 {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}

		 {\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)}.

	}
	}
	\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)

			& $m_I$
			& Tran(s)
			& Total(s)
			& Df(uV)
			& $m_R$
			& Tran(s)
			& Total(s)
			& Df(uV)

			& SPDP$^{r}_{tr}$
			& SPDP$^{r}_{i}$

			\\ \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$

			\\ \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$
			\\ \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$
			\\ \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$
			\\ \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$
	&  $0.9\times$
			& $3.5\times$
			\\
			\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
			\\ \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$
			\\ \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$
			\\ \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$
			\\
			\hline
			\hline
			 Average   & --- & ---
			  & --- & --- & --- & ---
			  & $65\backslash 15.7$ &  --- & ---
			 &  --- & $57\backslash 12.8$
		 & $5\times$
			 &$3\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.
We decompose the input sources by the approach discussed in Sec. \ref{bump} and obtain much fewer transitions of $LTS$ for computing nodes.
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.)
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.
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.
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.

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.
			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);

			{\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.

			{\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}.

		}
	}
	\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$_{r}$
		& $m$
		& over TR-FTS

		\\
		\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$

		&   $4.0\times$
		& 7 & 1.9
		\\
		\hline
		ibmpg3t & 100 & 91.7  & 129.9

		&  9.6e-5 & 19.7e-6   &  $2.9\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$

		&   $2.6\times$   &  7
		& 1.5
		\\
		\hline
		ibmpg6t  & 100
	& 189.9  & 234.2	&  3.4e-5  &  7.2e-6
		&  $3.8\times$
		& $4.1\times$   & 7
		& 1.5
		\\
		\hline
		ibmpg1t\_new  &  {\color{black} 100}
		&  2.4   & 2.8
		&  5.3e-5   & 8.6e-6
		&  $21.8\times$

		& $5.0\times$  & 6  & 1.9
		\\
		\hline
		ibmpg2t\_new &{\color{black} 100}

		& 5.6 & 7.0

		& 4.6e-5 & 8.6e-6
		& $61.6\times$

		 & $6.2\times$ &  7 & 1.9
		\\
		\hline
		ibmpg3t\_new  &  {\color{black} 100} & 	103.0 & 140.9

		&  9.8e-5 & 19.9e-6
		&  $25.6\times$

		 &  $3.3\times$

		& 7  & 1.5
		\\
		\hline
		ibmpg4t\_new & {\color{black} 15}
	& 51.5 & 108.4

		&  9.9e-5  & 27.6e-6
		& {$98.0\times$}

		  &  $6.8\times$  & 8   & 1.4
		\\
		\hline

		ibmpg5t\_new  &  {\color{black} 100}
		& 185.6 & 227.8
	 & 9.9e-5  &  2.2e-6
		& $27.3\times$

		 & $4.0\times$
		& 7 &  1.5
		\\
		\hline
		ibmpg6t\_new  & {\color{black} 100}
	&  274.8
&  317.7
		&  3.4e-5  &  7.1e-6
         &  $25.5\times$
		 &  $3.3\times$   &  7  & 1.5
		\\
		\hline 	\hline
		 {Average}
		& --- & --- & --- &  {7.1e-5} &  {12.3e-6}
		& { $26\times$}

		 &  { $5 \times$} &  {6.7}
		&
		 {1.6}
		\\
		\hline
	\end{tabular}
	\label{tab:pg_rmatex}
\end{table*}

\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}

\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,
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}

\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(\%)
	& Speedup
	& Stiffness\\

	\hline
	MEXP	  	&   229
	&  0.510
	& -- &
	\multirow{3}{*}{$2.1 \times 10^{16}$ }
	\\
	\cline{1-4}
	I-MATEX	 	&14       &  0.004
	& 2616X   &   \\
	\cline{1-4}
	R-MATEX	   & 12   &  0.004
	&  2735X   & \\
	\hline
	\hline
	MEXP	  	&   224  	&  0.004

	& --
	& 	    \multirow{3}{*}{$2.1 \times 10^{12}$ } \\
	\cline{1-4}
	I-MATEX		& 14       &0.004

	&  583X  &   \\
	\cline{1-4}
	R-MATEX	&  12   &  0.004
	&  611X  & \\
	\hline
	\hline
	MEXP	   	&   223  & 0.004
	&  -- & 	    \multirow{3}{*}{$2.1 \times 10^{8}$ }
	\\
	\cline{1-4}
	I-MATEX  	& 14      & 0.004
	& 229X  &   \\
	\cline{1-4}
	R-MATEX	  & 12   &  0.004
	&  252X  & \\

	\hline
    \end{tabular}
    \label{tab:rc_cases}
\end{table}
\end{comment}

\begin{comment}
\begin{table*}[htp]
    \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)

			     &Total(s)
			     &LU(s)
			     &Total(s)
			     &Total(s)
			     &SPDP\footnote{ over TR

			     &LU(s)	&Total(s)

			     &Total(s)
			     &SPDP\footnote{ over I-MATEX(mn) }
			     & SPDP\footnote{ over TR
			     & SPDP\footnote{ over TR (adpt)
			     \\ \hline
	ibmpg1t		     & 0.13
			      &0.13   &  6.20
			      & 29.48
			     & 0.12  & 27.23

			     & 4.68
			     & 1.3X
			     & 0.13 & 4.93

			     & 1.10
			     & 4.5X
			     & 5.6X
			     & 26.8X
			     \\ \hline
    	ibmpg2t 	     & 0.80
			    & 0.83   & 28.61
			    & 179.74
			   & 0.76   &  124.44

			   & 20.19
			   & 1.4X
			   & 0.85
			   & 25.90

			   & 5.56
			    & 3.6X
			    & 5.1X
			    & 32.3X
			    \\ \hline
    	ibmpg3t 	     & 13.83   & 12.72  & 272.47
			    & 2792.96
			   & 12.87 & 1246.69

			   & 213.33
			   & 1.3X
			   & 12.73
			   & 244.56

			   & 63.78
			   & 3.3X
			   & 4.3X
			   & 43.8X
			   \\ \hline
    	ibmpg4t 	    & 16.69    & 22.50   &	368.55
			    & 1773.14
			    & 21.09 & 484.83

			    & 148.45
			    & 2.5X
			    & 22.09  & 140.42

			    & 76.55
			    & 1.9X
			    & 4.8X
			    & 23.1X
			    \\ \hline
    	ibmpg5t 	    & 8.16    & 10.50  & 350.33
			   & 2359.11
			    & 9.92 & 1821.60

			    & 304.72
			    & 1.1X
			    & 10.77 &  337.74

			    &  78.99
			    &    3.9X
			    & 4.4X
			    & 29.9X
			    \\ \hline
    	ibmpg6t 	     & 11.17    &14.17  & 567.38
		    	& 3184.59
			    & 13.19 & 2784.46

			    & 479.79
			   & 1.2X
			    & 14.18 & 482.42

			    & 108.92
			   & 4.4X
			    &   5.2X
			   & 29.2X
			    \\
    \hline
    \end{tabular}

    \end{threeparttable}
    \label{tab:perf_mexp}
\end{table*}
\end{comment}

\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 & 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}

\end{table*}

\end{comment}
\subsection{MATEX Circuit Solver}
\label{sec:frame_motivation}

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 $\mathbf X_1$ and $\mathbf X_2$ in Line \ref{algo:rmatex:init} based on the choice of Krylov subspace method as follows,
\begin{itemize}
\item I-MATEX: $\mathbf X_1 =  \mathbf G, ~\mathbf X_2=\mathbf C$
\item R-MATEX: $\mathbf X_1 =  \mathbf C+ \gamma \mathbf G, ~\mathbf X_2=\mathbf C$
\end{itemize}

For linear system of PDN, the matrix factorization in line \ref{algo:rmatex:lu} is only performed once, and the matrices $\mathbf L$ and $\mathbf U$ are reused in the while loop from line \ref{algo:rmatex:while} to line \ref{algo:rmatex:whileend}.
Line \ref{algo:rmatex:algo1} uses Arnoldi process  with corresponding inputs to construct Krylov subspace as shown in Algorithm \ref{algo:arnoldi}.

\begin{algorithm}[t]
 	\label{algo:rmatex}

 	{

 		\caption{MATEX Circuit Solver}

 		\KwIn{ $ \mathbf C, \mathbf G , \mathbf B, \mathbf u, \epsilon$,

 			and time span $T$.

 		}

 		\KwOut{  The set of $\mathbf x$ from $[0, T]$.}
 		{Set $\mathbf X_1, \mathbf X_2$\;
 	\label{algo:rmatex:init}

 			$t=0$\;
 			$\mathbf x(t) =$DC\_analysis\;
 			$[\mathbf L, \mathbf U] = \text{LU\_Decompose}(\mathbf X_1)$\; 	\label{algo:rmatex:lu}
 			\While { $t < T$}
 			{\label{algo:rmatex:while}
 				Compute maximum allowed step size $h$\;
 				Update $ \mathbf P (t,h),\mathbf F (t,h)$\;
 				Obtain  $\mathbf x(t+h)$ by  {\bf Algorithm  \ref{algo:arnoldi}} with inputs $[\mathbf L, \mathbf U, \mathbf X_2, h ,t, \mathbf x(t), \epsilon, \mathbf P (t,h),\mathbf 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}
\label{sec:dist_frame}
\subsubsection{Motivation} \

There are usually many input sources in PDNs as well as their transition activities, which might narrow the regions for the stepping of matrix exponential based method due to the unaligned breakpoints.
In other words, the region before the next transition $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.}
\vspace{+0.07in}

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

In matrix exponential based integration framework,
we can choose any time spot  $t+h \in [t,t_s]$ with computed Krylov subspace basis. The solution of $\mathbf x(t+h)$ is computed by scaling the existing Hessenberg matrix $\mathbf H$ with the time step $h$ as below
\begin{eqnarray}
\mathbf x(t+h) = \lVert \mathbf v \rVert \mathbf {V}_m e^{h \mathbf H} \mathbf e_1-\mathbf 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.
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

	\emph{Local Transition Spot} ($LTS$): the set of $TS$ at an input source to the PDN.
	\item

	\emph{{Global Transition Spot}} ($GTS$): the union of $LTS$ among all the input sources to the PDN.
	\item

	\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$.
For each node, the chances of generation of Krylov subspaces are reduced.
The time periods of reusing latest Krylov subspaces are enlarged locally and bring the runtime improvement.
Besides, when subtasks are assigned, there is no communication among the computing nodes, which leads to so-called \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_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/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_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$, $\mathbf P_k$ and $\mathbf F_k$, which contain the corresponding $\mathbf 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$, $\mathbf P_k$, $\mathbf F_k$, error tolerance $E_{tol}$, and simulation time span $T$.}

		\KwOut{Local solution $\mathbf x$ along $GTS$ in node
			$k \in [1,\cdots,S]$, where $S$ is the number of nodes}
		{

			$t=0$, $\mathbf X_1= \mathbf C + \gamma \mathbf G$, and $\mathbf X_2 = \mathbf C$\;

			$\mathbf x(t) = \text{Local\_Initial\_Solution}$\;

			$[\mathbf L,\mathbf U] = \text{LU\_Decompose}( \mathbf X_1)$\;
			\While { $t < T$}
			{
				Compute maximum allowed step size $h$ based on $GTS$\;

				{

					\If{$t \in $  LTS\#$k$}
					{
						\tcc{\small  Generate Krylov subspace for the point at $LTS\#k$ and  compute $\mathbf x(t+h)$}

						$[\mathbf x(t+h), \mathbf V_m, \mathbf H_m,\mathbf v] = \text{MATEX\_Arnoldi}(\mathbf L,\mathbf U,\mathbf X_2,$ $ h, t, \mathbf x(t), \epsilon, \mathbf P_k(t,h), \mathbf F_k(t,h))$\;

						$a_{lts} = t$\;

					}
					\Else{
						\tcc{\small Obtain $\mathbf x(t+h)$ at $Snapshot$ with computed Krylov subspace}

						{

							$h_{a} = t + h - a_{lts}$\;
							$\mathbf x(t+h) = \lVert \mathbf v \rVert \mathbf {V}_m e^{h_{a} \mathbf H_m} \mathbf e_1 - \mathbf 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 $\mathbf x(t+h)$, the evaluation of matrix exponential with $e^{h\mathbf H_m}$ is $T_{H}$, which is in proportion to the time complexity $O(m^3)$. Besides, we need extra $T_e$ to form $  \mathbf x(t+h)$, which is proportional to $O(n m^2)$ by $\beta\mathbf V_m e^{h\mathbf H_m} \mathbf 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
\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.
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}}.

\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}

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.

\subsection{MEXP: MEVP Computation via Standard Krylov Subspace Method}
\label{sec:std}
The complexity of $e^{\mathbf A}\mathbf v$ can be reduced
using Krylov subspace method and still maintained
in a high order polynomial approximation \cite{Saad92}.
MEXP \cite{Weng12_TCAD} uses standard Krylov subspace,
which uses $\mathbf 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}
\mathbf{x}(t+h)
&=&
\mathbf{y}(t+h)
-(\mathbf A^{-1}\mathbf b(t+h) + \mathbf A^{-2} \frac{\mathbf b(t+h)- \mathbf b(t)}{h})
\nonumber \\
\mathbf{y}(t+h)
& =& e^{h\mathbf{A}}\mathbf{y}(t)
\end{eqnarray}
\end{comment}
\begin{eqnarray}
	\label{eqn:new_exact}
	\mathbf{x}(t+h)
	=
	e^{h\mathbf A}(\mathbf x(t) + \mathbf F(t,h) ) - \mathbf P (t,h),
\end{eqnarray}
where
\begin{eqnarray}
	\mathbf F(t,h) =  \mathbf A^{-1}\mathbf b(t) + \mathbf A^{-2} \frac{\mathbf b(t+h)- \mathbf b(t)}{h}
\end{eqnarray} and
\begin{eqnarray}
	\mathbf P(t,h) =   \mathbf A^{-1}\mathbf b(t+h) + \mathbf A^{-2} \frac{\mathbf b(t+h)- \mathbf b(t)}{h}  .
\end{eqnarray}
The standard Krylov subspace
\begin{eqnarray}
	{\rm K_m}(\mathbf A, \mathbf v) := \text{span}
	{\{ \mathbf v, \mathbf A\mathbf v,
		\cdots, \mathbf A^{m-1} \mathbf v}\}
\end{eqnarray}
obtained by Arnoldi process
has the relation
\begin{eqnarray}
	\mathbf{A}\mathbf{ V_m} =
\mathbf {   V_m} \mathbf {  H_m}
	+ {  h_{m+1,m} }
	\mathbf {  v}_{m+1}\mathbf e_m^\mathsf{T},
\end{eqnarray}
{ where
$\mathbf {H_m}$ is the upper Hessenberg matrix
\begin{eqnarray}
\mathbf {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}
$\mathbf { V_m}$ is a $n \times m $ matrix by $(\mathbf v_1, \mathbf v_2, \cdots, \mathbf v_m)$, and $\mathbf e_m$ is the $m$-th unit vector with dimension $n\times 1$.
MEVP is computed via
\begin{eqnarray}
	e^{h\mathbf A} \mathbf v \approx \beta \mathbf {  V_m} e^{h\mathbf {  H}_m} \mathbf e_1 .
\end{eqnarray}
The posterior error term is
\begin{eqnarray}
	\label{eqn:std_error}
   r_m(h)  =  \lVert \beta { h_{m+1,m}} \mathbf v_{m+1}
  	\mathbf e^{\mathsf{T}} e^{h\mathbf{  H}_m} \mathbf e_1 \rVert,
\end{eqnarray}
where $\beta = \lVert  \mathbf v \rVert $.
However, for an autonomous system $ \mathbf C \dot {\mathbf x}(t) = - \mathbf G \mathbf x(t)$  in circuit simulation, we consider the residual between
$\mathbf C \dot {\mathbf x}(t) $ and $-\mathbf G \mathbf  x(t)$,
which is
$$ \mathbf C \dot {\mathbf x}(t) + \mathbf G \mathbf  x(t) ,$$
instead of  $$ \dot {\mathbf x}(t) - \mathbf A \mathbf  x(t)$$  in
$ \dot {\mathbf x}(t) = \mathbf A \mathbf x(t)$.
This leads to
\begin{eqnarray}
	\label{eqn:std_err}
	   r(m,h)  =    \lVert\beta  h_{m+1,m}  \mathbf C  \mathbf
v_{m+1}	\mathbf e^{\mathsf{T}} e^{h\mathbf{  H}_m} \mathbf e_1 \rVert
\end{eqnarray}
and helps us mitigate the overestimation of the error bound.

To generate $\mathbf x(t+h)$ by Algorithm \ref{algo:arnoldi},
we use
\begin{eqnarray}
\label{eq:lu}
[\mathbf L, ~\mathbf U] = \text{LU\_Decompose}(\mathbf X_1),
\end{eqnarray}
where $$\mathbf X_1 = \mathbf C ~ \text{and} ~ \mathbf X_2 = \mathbf G$$ as inputs for standard Krylov subspace.
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}).
\begin{algorithm}
	\label{algo:arnoldi}

	\caption{MATEX\_Arnoldi}

	\KwIn{ $\mathbf L, \mathbf U, \mathbf X_2, h ,t, \mathbf x(t),  \epsilon, \mathbf P (t,h),\mathbf F (t,h)$  }

	\KwOut{ $\mathbf x(t+h), \mathbf V_m, \mathbf H, \mathbf v$}
	{

		$\mathbf v = \mathbf x(t) + \mathbf F(t,h)$\;
		$\mathbf v_1 = \frac{\mathbf v }{\lVert \mathbf v \rVert}$\;
		\For  {$j=1:m$}
		{

			\label{arnoldi_orth}
			$\mathbf w = \mathbf U \backslash  (\mathbf L \backslash (\mathbf X_2 \mathbf v_{j}))$
			{\tcc*[l]{a pair of forward and backward substitutions.}}
			\For {$i = 1:j$}
			{

				$h_{i,j} = \mathbf w^T\mathbf v_{i}$\;
				$\mathbf w = \mathbf w - h_{i,j} \mathbf v_{i}$\;
			}
			$h_{j+1,j} = \lVert \mathbf w \rVert $\;
			$\mathbf v_{j+1} = \frac{\mathbf w }{h_{j+1,j}} $\;

			\If {$r(j,h)  < \epsilon $\label{alg1:conv1}}
			{
				$m = j$\;
				break\;
			}\label{alg1:conv2}
		}

		$\mathbf x(t+h) = \lVert \mathbf v \rVert \mathbf {V_m} e^{h \mathbf H} \mathbf e_1 - \mathbf P(t,h)$\;
		\label{alg1:h}

	}
\end{algorithm}

\label{sec:discuss}

The standard Krylov subspace may not be efficient when simulating stiff circuits \cite{Weng12_TCAD, Weng12_ICCAD}.
For the accuracy of approximation of $e^{\mathbf A}\mathbf 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 $\mathbf H_{m}$ of standard Krylov subspace tends to approximate the large magnitude eigenvalues of $\mathbf 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}.
Therefore, to simulate stiff circuit, we need to gather more vectors into subspace basis
and increase the size of $\mathbf H_{m}$ to fetch more useful components, which results in both
memory overhead and computational complexity to Krylov subspace generations for each time step.
In the following subsections, we adopt  ideas from \emph{spectral transformation} \cite{Botchev12, Van06} to effectively capture small magnitude eigenvalues in $\mathbf A$, leading to a fast and accurate MEVP computation.

\subsection{I-MATEX: MEVP Computation via Invert Krylov Subspace Method}
\label{sec:imatex}
Instead of $\mathbf A$, we use $\mathbf A^{-1}$ (or $\mathbf G^{-1} \mathbf C$) as our target matrix to form
\begin{eqnarray}
	{\rm K_m}(\mathbf A^{-1} , \mathbf v) := \text{span}
	{\{ \mathbf v, \mathbf A^{-1}\mathbf v, \cdots, \mathbf A^{-(m-1)} \mathbf v}\}.
\end{eqnarray}
Intuitively, by inverting $\mathbf A$, the small magnitude eigenvalues become the large ones
of $\mathbf A^{-1}$.
The resulting $\mathbf H_{m}$ is likely to capture these
eigenvalues first. Based on Arnoldi algorithm, the invert Krylov subspace has the relation
\begin{eqnarray}
	\mathbf{A}^{-1}\mathbf{V_m} =
	\mathbf V_m \mathbf H_m +
	h_{m+1,m} \mathbf v_{m+1}\mathbf e_m^\mathsf{T}.
\end{eqnarray}
The matrix exponential $e^{\mathbf A} \mathbf v$ is calculated as
\begin{eqnarray}
	e^{\mathbf A} \mathbf v \approx \beta \mathbf V_m e^{h\mathbf {H}^{-1}_m} \mathbf e_1.
\end{eqnarray}
To put this method into Algorithm 1 is just by modifying the inputs
$\mathbf X_1 = \mathbf G$ for the LU decomposition in Eq. (\ref{eq:lu}),
and $\mathbf X_2 = \mathbf C$.
In the line \ref{alg1:h} of Algorithm \ref{algo:arnoldi},
$$\mathbf H = \mathbf {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}  \mathbf A  \mathbf v_{m+1}
	\mathbf e^\mathsf{T}_m {\mathbf H}^{-1}_{m}
	e^{h {\mathbf H}}
	\mathbf 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(\mathbf C \dot {\mathbf x}(t) + \mathbf G \mathbf x(t) \right)$,
instead of $\left( \dot {\mathbf x}(t) - \mathbf A \mathbf  x(t) \right)$,
which leads to
\begin{eqnarray}
	\label{eq:er_inverted_krylov}
   r(m,h)
	=
	\lVert   \beta h_{m+1,m} \mathbf G   \mathbf v_{m+1}
	\mathbf e^\mathsf{T}_m {\mathbf H}^{-1}_{m}
	e^{h {\mathbf H}_{m}}
	\mathbf e_1  \rVert .
\end{eqnarray}
{We use Eq. (\ref{eq:er_inverted_krylov}) for the line \ref{alg1:conv1} of Alg. \ref{algo:arnoldi}.

\subsection{R-MATEX: MEVP Computation via Rational Krylov Subspace Method}
\label{krylov_rat_err}
The shift-and-invert Krylov subspace basis
\cite{Van06}
is designed
to confine the spectrum of
$\mathbf A$.
Then, we generate Krylov subspace via
\begin{eqnarray}
	& &  \mathbf {K_m} ((\mathbf I- \gamma \mathbf {A} )^{-1}, \mathbf v) :=  \\
	& &  \text{span} \{ \mathbf v, (\mathbf I- \gamma\mathbf { A} )^{-1}\mathbf v,
	\cdots,

	(  \mathbf I -  \gamma \mathbf{A})^{-(m-1)} \mathbf v \}, \nonumber
\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.
The similar idea has been applied to simple power grid simulation
with matrix exponential method
\cite{Zhuang13_ASICON}.
Here, we generalize this technique and integrate into MATEX.
The Arnoldi process constructs $\mathbf{V_m}$ and $\mathbf{H_m}$.
We have
\begin{eqnarray}
	(\mathbf{I}-  \gamma \mathbf {A})^{-1} \mathbf{ V}_m =
	\mathbf{ V}_m\mathbf{ {H}}_m +  {h}_{m+1,m} \mathbf { v}_{m+1} \mathbf e^\mathsf{T}_m .
\end{eqnarray}
We can project the $e^{\mathbf{{A}}}$ onto the rational Krylov subspace as follows.
\begin{eqnarray}
	\label{eq:scale}
	e^{\mathbf{A}h} \mathbf{v}  \approx
	 \beta \mathbf{  V}_m e^{h\frac{\mathbf I - \mathbf { {H}}_m^{-1}}{\gamma}}\mathbf e_1.
\end{eqnarray}
In the line \ref{alg1:h} of Algorithm \ref{algo:arnoldi},
$${\mathbf{H}} = \frac{\mathbf I - \mathbf { {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{\mathbf I - \gamma \mathbf A_m}{\gamma}
	\mathbf {  v}_{m+1}
	\mathbf e^\mathsf{T}_m \mathbf{  H}_m^{-1}

	e^{ h\mathbf{H}}
	\mathbf e_1  \rVert .
\end{eqnarray}
Note that in practice, instead of
computing $(\mathbf I-\gamma \mathbf{A})^{-1}$ directly,
$(\mathbf C + \gamma \mathbf G)^{-1}\mathbf C $ is utilized.
The corresponding Arnoldi process shares
the same skeleton of Algorithm \ref{algo:arnoldi}
with input matrices
$$\mathbf X_1 = (\mathbf C+ \gamma \mathbf G)$$ for the LU decomposition Eq. (\ref{eq:lu}),
and $$\mathbf X_2 = \mathbf {C}.$$
The residual estimation is
\begin{eqnarray}
	\label{eq:er_rational_krylov}
  r(m,h)
	= 	\lVert     \beta
	{h}_{m+1,m}
	\frac{\mathbf C  + \gamma \mathbf G}{\gamma}
	\mathbf {  v}_{m+1}
	\mathbf e^\mathsf{T}_m \mathbf{  {H}}_m^{-1}

	e^{ h\mathbf{  H}_m}
	\mathbf e_1 \rVert .
\end{eqnarray}
Then, we plug Eq. (\ref{eq:er_rational_krylov}) into the line \ref{alg1:conv1} of Algorithm \ref{algo:arnoldi}.
}

\subsection{Regularization-Free MEVP Computation}
When $\mathbf 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 $\mathbf C$ directly to form the input $\mathbf X_1$ for Algorithm 1.
This brings extra computational overhead when the case is large \cite{Chen12_TCAD}.
It is not necessary if we can obtain the generalized eigenvalues and corresponding eigenvectors for matrix pencil $(-\mathbf G, \mathbf C)$.
Based on \cite{Wilk:79}, we derive the following lemma,
\begin{lemma}
	Considering a homogeneous system

	$$ \mathbf C \mathbf {  \dot{x}} = -\mathbf G \mathbf x.$$

	$\mathbf u$ and $\lambda$ are the eigenvector and eigenvalue
	of matrix pencil $(-\mathbf G, \mathbf C)$, then $$\mathbf x =  e^{t\lambda} \mathbf u$$ is a solution of the system.
\end{lemma}
{
\begin{proof}\footnote{We repeat the proof from \cite{Wilk:79} with some modifications for our formulation.}
If $\lambda$ and $\mathbf u$ are an eigenvalue and eigenvector of a generalized eigenvalue problem $$-\mathbf G \mathbf u = \lambda \mathbf C \mathbf u.$$
Then, $ \mathbf x = e^{t\lambda} \mathbf u$ is the solution of $\mathbf C \mathbf {\dot x} =  - \mathbf G \mathbf x$.
\end{proof}}
\vspace{+0.07in}
Because we do not need to compute $\mathbf C^{-1}$ explicitly during Krylov subspace generation, I-MATEX and R-MATEX are regularization-free. Instead, we factorize $\mathbf G$ for invert Krylov subspace basis generation (I-MATEX),
or $(\mathbf C+\gamma\mathbf 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 $(-\mathbf G, \mathbf 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 $\mathbf A = -\mathbf C^{-1} \mathbf 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\mathbf A}\mathbf v - \beta \mathbf V_m e^{h\mathbf H_m}\mathbf e_1||}{||e^{h\mathbf A}\mathbf v ||},$$
where $h=0.4ps$, $\gamma = 10^{-13}$.
The matrix $\mathbf A$ is a relatively small matrix and computed by MATLAB $expm$ function. The result of $e^{h\mathbf A}\mathbf v$ serves as the baseline for accuracy.
The relative error is the real relative difference compared to the analytical solution $e^{h\mathbf A}\mathbf v$ of the ODE $$\frac{d\mathbf x}{dt}= \mathbf A \mathbf x$$ with an initial vector $\mathbf 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 $\mathbf C$ and $\mathbf 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{

	{The relative error vs. dimensional $m$ of different Krylov subspace methods. The relative error is $\frac{||e^{h\mathbf A}\mathbf v - \beta \mathbf V_m e^{h\mathbf H_m}\mathbf e_1||}{||e^{h\mathbf A}\mathbf v ||}$, where $h=0.4ps$, $\gamma = 10^{-13}$. Note: The relative error is the  difference compared to analytical solution $e^{h\mathbf A}\mathbf v$ of the ODE $ \frac{d\mathbf x}{dt}= \mathbf A \mathbf x$ with an initial vector $\mathbf 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\mathbf A}\mathbf v - \beta \mathbf V_m e^{h\mathbf H_m}\mathbf e_1||}{||e^{h\mathbf A}\mathbf 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$.
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\mathbf A}\mathbf v -  \beta \mathbf {  V_m} e^{h\mathbf {  H}_m}||}{||e^{h\mathbf A}\mathbf 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 $\mathbf 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\mathbf A}\mathbf v -  \beta \mathbf {V}_m e^{h\mathbf {H }^{-1}_m}||}{||e^{h\mathbf A}\mathbf 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\mathbf A}\mathbf v - \beta \mathbf {  V}_m e^{h \frac{\mathbf I - \mathbf { {H}}_m^{-1}}{\gamma}}\mathbf e_1||}{||e^{h\mathbf A}\mathbf 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.

\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\mathbf A}\mathbf v - \beta \mathbf {  V}_m e^{h \frac{\mathbf I - \mathbf { {H}}_m^{-1}}{\gamma}}\mathbf e_1||}{||e^{h\mathbf A}\mathbf 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}

\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}
    \mathbf{C}\dot{\mathbf{x}}(t) = - \mathbf{G}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t),
\end{eqnarray}
where  $\mathbf{C}$ is the  matrix for  capacitive and inductive elements. $\mathbf{G}$ is the matrix for conductance and resistance, and $\mathbf{B}$ is the input selector matrix.  $\mathbf{x}(t)$ is the vector of time-varying node voltages
and branch currents. $\mathbf{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 $\mathbf x(0)$ from DC analysis or previous time step $\mathbf x(t)$ and a time step $h$, $\mathbf 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{\mathbf{C}}{h} +
     {\mathbf G} \Bigg) {\mathbf{x}}(t+h)
    =  \frac{\mathbf C}{h}
    \mathbf{x}(t) + \mathbf{B}     \mathbf{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{\mathbf{C}}{h} +
    \frac{\mathbf G}{2}\Bigg) {\mathbf{x}}(t+h)    & = &
    \Bigg(\frac{\mathbf C}{h}- \frac{\mathbf G}{2} \Bigg)
    \mathbf{x}(t)  \\ &+& \mathbf{B} \frac{  \mathbf{u}(t) + \mathbf{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}.
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.

\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]
	\centering
	\includegraphics[trim = 3.5cm 4cm 5cm 0cm, clip,keepaspectratio, width=0.45\textwidth]{./figs/interleaving_new.eps}

	\caption
		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$.
	}
	\label{fig:interleave}
\end{figure}

\subsection{Matrix Exponential Time Integration Scheme}

The solution of Eq. (\ref{eqn:dae}) can be obtained analytically \cite{Leon75}.
For a simple illustration, we convert Eq. (\ref{eqn:dae}) into
\begin{eqnarray}
\label{eqn:ode}
    \dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{b}(t),
\end{eqnarray}
when $\mathbf{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 $\mathbf C$. Therefore, the methods are suitable for general DAE system, i.e., Eq. (\ref{eqn:dae}) without the assumption here.},
$$\mathbf{A} = - \mathbf{C}^{-1}\mathbf{G}~\text{, and}~\mathbf{b}(t) = \mathbf{C}^{-1}\mathbf{B}\mathbf{u}(t).$$
Given  a solution at time
$t$ and a time step $h$, the solution at
$t + h$ is
\begin{eqnarray}
    \label{eqn:discrete_sol}
    \mathbf{x}(t+h) = e^{h\mathbf{A}}\mathbf{x}(t) +
    \int^{h}_{0}e^{(h-\tau)\mathbf{A}}\mathbf{b}(t+\tau)d\tau.
\end{eqnarray}

Assuming that the input $\mathbf{u}(t)$ is a piecewise linear (PWL) function of $t$,
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}
    \mathbf{x}(t+h)  \nonumber =

      -\left(\mathbf A^{-1}\mathbf b(t+h) + \mathbf A^{-2} \frac{\mathbf b(t+h)- \mathbf b(t)}{h}\right) +
        \\

       e^{h\mathbf{A}}\left(\mathbf{x}(t) + \mathbf A^{-1}\mathbf b(t) + \mathbf A^{-2} \frac{\mathbf b(t+h)- \mathbf b(t)}{h}\right )  .

\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 $\mathbf 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.

\section*{Categories and Subject Descriptors}
B.7.2 [Integrated Circuits]: Design Aids, Simulation.
J.6 [Computer-aided engineering]: Computer-aided design (CAD).

\section*{General Terms}
Algorithms,  Design, Theory, Verification, Performance
\end{comment}
\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

\section{Introduction}
\label{sec:intro}
\input{intro_ljw}

\section{Background}
\label{sec:method}
\input{matex_method}

\section{Krylov Subspace Methods for Matrix Exponential and Vector Product (MEVP)}
\label{sec:advanced}
 \input{krylov}

\section{MATEX Framework}
\label{sec:matex}
\input{framework}

\section{Experimental Results}
\label{sec:exp_results}
\input{exp_results}

\section{Conclusions and Future Directions}
\label{sec:conclusion}
\input{sub-concl}

\section*{Acknowledgment}
\input{ack}
\bibliographystyle{ieeetr}

\section*{Categories and Subject Descriptors}
B.7.2 [Integrated Circuits]: Design Aids, Simulation.
J.6 [Computer-aided engineering]: Computer-aided design (CAD).

\section*{General Terms}
Algorithms,  Design, Theory, Verification, Performance
\end{comment}
\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

\section{Introduction}
\label{sec:intro}
\input{intro_ljw}

\section{Background}
\label{sec:method}
\input{matex_method}

\section{Krylov Subspace Methods for Matrix Exponential and Vector Product (MEVP)}
\label{sec:advanced}
 \input{krylov}

\section{MATEX Framework}
\label{sec:matex}
\input{framework}

\section{Experimental Results}
\label{sec:exp_results}
\input{exp_results}

\section{Conclusions and Future Directions}
\label{sec:conclusion}
\input{sub-concl}

\section*{Acknowledgment}
\input{ack}
\bibliographystyle{ieeetr}

\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}

\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(\%)
& Speedup
& Stiffness\\
\hline
MEXP	  	&   229
&  0.510
& -- &
\multirow{3}{*}{$2.1 \times 10^{16}$ }
\\
\cline{1-4}
I-MATEX	 	&14       &  0.004
& 2616X   &   \\
\cline{1-4}
R-MATEX	   & 12   &  0.004
&  2735X   & \\
\hline
\hline
MEXP	  	&   224  	&  0.004
& --
& 	    \multirow{3}{*}{$2.1 \times 10^{12}$ } \\
\cline{1-4}
I-MATEX		& 14       &0.004
&  583X  &   \\
\cline{1-4}
R-MATEX	&  12   &  0.004
&  611X  & \\
\hline
\hline
MEXP	   	&   223  & 0.004
&  -- & 	    \multirow{3}{*}{$2.1 \times 10^{8}$ }
\\
\cline{1-4}
I-MATEX  	& 14      & 0.004
& 229X  &   \\
\cline{1-4}
R-MATEX	  & 12   &  0.004
&  252X  & \\
\hline
\end{tabular}
\label{tab:rc_cases}
\end{table}
\end{comment}

\begin{comment}
\begin{table*}[htp]
\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)
&Total(s)
&LU(s)
&Total(s)
&Total(s)
&SPDP\footnote{ over TR
&LU(s)	&Total(s)
&Total(s)
&SPDP\footnote{ over I-MATEX(mn) }
& SPDP\footnote{ over TR
& SPDP\footnote{ over TR (adpt)
\\ \hline
ibmpg1t		     & 0.13
&0.13   &  6.20
& 29.48
& 0.12  & 27.23
& 4.68
& 1.3X
& 0.13 & 4.93
& 1.10
& 4.5X
& 5.6X
& 26.8X
\\ \hline
ibmpg2t 	     & 0.80
& 0.83   & 28.61
& 179.74
& 0.76   &  124.44
& 20.19
& 1.4X
& 0.85
& 25.90
& 5.56
& 3.6X
& 5.1X
& 32.3X
\\ \hline
ibmpg3t 	     & 13.83   & 12.72  & 272.47
& 2792.96
& 12.87 & 1246.69
& 213.33
& 1.3X
& 12.73
& 244.56
& 63.78
& 3.3X
& 4.3X
& 43.8X
\\ \hline
ibmpg4t 	    & 16.69    & 22.50   &	368.55
& 1773.14
& 21.09 & 484.83
& 148.45
& 2.5X
& 22.09  & 140.42
& 76.55
& 1.9X
& 4.8X
& 23.1X
\\ \hline
ibmpg5t 	    & 8.16    & 10.50  & 350.33
& 2359.11
& 9.92 & 1821.60
& 304.72
& 1.1X
& 10.77 &  337.74
&  78.99
&    3.9X
& 4.4X
& 29.9X
\\ \hline
ibmpg6t 	     & 11.17    &14.17  & 567.38
& 3184.59
& 13.19 & 2784.46
& 479.79
& 1.2X
& 14.18 & 482.42
& 108.92
& 4.4X
&   5.2X
& 29.2X
\\
\hline
\end{tabular}
\end{threeparttable}
\label{tab:perf_mexp}
\end{table*}
\end{comment}

\end{document}


Brittle Bones and Weak Knees, Oh! My . . .
Like most teens of the 60's and even the teens of today, I had a devine feeling of invulnerability during those early years of my life.
It never seemed possible that I could somehow injure myself to the extent that I couldn't participate in high school sports. I always had trouble imagining that I could not do just about anything physically I wanted to do or at least give it a try.
Now it is fully realized how fortunate it was that I was able to remain injury free long enough to play football, baseball, basketball and run track for the Dragons during my four years of high school. I must have played 1000 basketball games in high school, college, in the recreation leagues and in various pick-up games in gyms all over Jackson and Barrow County. As an old-time pole vaulter I continuously fell 10 - 12 feet into sawdust pits and hard ground on many occasions. The cushioney foam landing pit was a luxury of the future. When I think back on the things I did in high school I am amazed that I never received an injury that kept me from a scheduled contest or event. I served as a high school basketball referee for 26 years, carried the down marker during all Jefferson home football games for 24 years, and never missed a game due to an injury --- or illness for that matter. It was true that several of my teammates did receive injuries from their high school sports activities, but I never saw that as a ditriment to my continued participation. I always considered it a blessing that I was able to avoid this very possible situation. Still, I don't think I ever placed much thought into what could have happened until "the chickens came home to roost" in September of 2012. It is sad that it came about when I became a 67-year old man that is probably about as athletic now as a milk cow.
It was on September 23, during a golf match that has been a twice weekly activity with me for a number of years now. Much like when I was a high school athlete, it never dawned on me that I could be seriously injured while doing something that was as normal to me as getting up and going to work had been for over 38 years before my retirement.
On this fateful day I took a swing with a three-wood from a slightly down-hill lie. My weight shift was too much for the ligaments and meniscus in my left knee to withstand. I went down in tremendous pain that left me unable to walk. I had suffered a torn meniscus and some serious ligament damage.
My Athens, Georgia bone and joint physician told me after an MRI that it would have to be repaired if I ever wanted to walk without a major limp. Then he really caught my attention when he said I could possibly never play the game of golf again due to the constant pain and stiffness in my left leg without major surgery.
I don't mind admitting that having my knee cut on and repaired gave me much time to think. I questioned why something like this would come up when I reached the age of 67 years old. How was I able to escape this type of thing earlier in my life.
As I mentioned, I never had any type of serious injury. No broken bones or sprains that kept me limited as to what I could do. Why, Now? Dr. Ryan, whom I might add didn't sugar coat my injury as he carefully explained that the chances for a complete recovery, at my age, would be extremely difficult. Some of the factors that would make a complete recovery difficult were age, worn Cartilage between my knee joints, arthritus and other factors that would work against me along the road to complete recovery. My rehab would have to be taken seriously with a total dedication to improvement if I were to regain the former flexibility and strength in my left knee.
On November 27, I underwent orthroscopic surgery on my knee at the Athens Surgical Center. Three days later, I was able to place some weight on the knee, but there was tremendous swelling of my leg from toes to hip. On my post-op visit to my doctor on December 4, I had two six-inch vials of fluid ( mostly blood ) drawn off the knee. The last several days I have experienced less pain than at anytime since the accident.
I start rehab next Friday for two hours. I will undergo this process twice weekly for two months. Hopefully by March, I'll be able to tee a golf ball up and give it a try on the links again. I miss my weekly golf matches. I'll rehab as necessary will full intensity to regain full use of my knee. Will I ever trust the knee to allow me to replicate the same kind of monster swing at the ball as I did in the past? I really won't know how successful my return to golf will be until I give it a try. I hope I'll be back to normal or close to normal by this time.
One thing I do know. I plan on being more careful than in the past to protect myself from self injury. It's been the most difficult challenge I have ever faced --- not being able to get out and do the things I enjoy doing is really frustrating. It is also a humbling and unforgettable experience.
Would I have preferred to have an injury such as this years ago while in my high school/college days if I had the choice? I don't believe so. Now I have all the time in the world to recover without losing the sports elgibility I would have lost if injured in high school or college. I guess there was just too much going on in the life of this 1960's Jefferson teenager to want to lose precious time due to an injury.
Thank you, Lord. In my mind, this was a case of better late than early.
Golf ~ The safe game . . . You can play it all your life . . . Sure you can!
Happier golf days. Shortly after making a hole-in-one at the Winder course in mid-July. It was with a seven iron on the 153-yard par 3, number 3 hole.
This is how the meniscus should look in a photo taken of the left side of my left knee. It was not injured in the accident.
Another view of the injury to the inside of my left knee.
The trimmed meniscus after surgery.
This is a view of my torn meniscus from my September golf injury. I underwent arthroscopic knee surgery on November 27 at Athens Surgical Center by Dr. Ryan of Athens Bone and Joint Service. The surgery went well. It seems recovery is well under way at this point.
The injured meniscus after trimming away the injured portion and preparing it for forming and compression into the knee joint.
A view of my knee ( Leg ) after removal of bandaging 72 hours after surgery. Tremendous swelling and toothache-like pain was customary without strong pain killers.
A before and after surgery view of my torn meniscus.


text:
West Indies Cricket unveil jersey for World Cup 2019
For Many, West Indies Aka Windies Could Be The Dark Horse Of The ICC World Cup 2019 Which Is Set To Begin On May 30 As The Host-England Take On South Africa At The Oval.
News Nation Bureau | Edited By : Vedant Sharma | Updated on: 23 May 2019, 11:08:41 AM
West Indies Cricket unveil jersey for World Cup 2019 (Image Credit: Twitter)
For many, West Indies aka Windies could be the dark horse of the ICC World Cup 2019 which is set to begin on May 30 as the host-England take on South Africa at the Oval. Even though the two-time champions had to qualify for the tournament this year and got through as the last team after Afghanistan, the Windies could be the team to look for this time. Their recent performance as a team against England at home and the individual brilliance in the Indian Premier League (IPL) 2019 have changed the way the cricketing world is looking at them now.
On Wednesday, the official World Cup kit of the Windies was launched and it was a complete maroon one with the logo of the team and the World Cup printed on it. The name 'West Indies' is written in yellow in front of the jersey while there is a shadow of the palm tree in front, signifying the Caribbean flavour. The gear is also available at the official Cricket West Indies website.
The West Indies, led by all-rounder Jason Holder for the second consecutive time, will feel boosted by the fact that the World Cup this year is taking place in England. It was in the same country that they had reached the final of the World Cup three times on the trot and won it twice (1975 and 1979). Interestingly, no other team has played so many finals in the World Cup in England.
The team also features some special talents in the form of Andre Russell and Chris Gayle who on their days can decimate any opposition and that has made the fans of the team hopeful about a grand turnaround. The Windies are currently ranked eighth in the ICC's ODI rankings and have not made the World Cup semi-final since 1996.
The West Indies will begin their campaign at the World Cup against Pakistan who they will play at Trent Bridge, Nottingham, on May 31. Although Pakistan is ranked higher than the Caribbeans at the moment, the latter have a 7-3 lead over them in the World Cup. They had beaten them in the semi-finals of the 1979 and 1983 editions played in England.
Here's the Twitter post:
Coming soon! 👀 Swipe 👈🏽 to see OFFICIAL WEST INDIES 😍🌴#CWC19 gear! #MenInMaroon #WIAllin pic.twitter.com/XkJy4Ufs9y
— Windies Cricket (@windiescricket) May 21, 2019
ICC World Cup 2019 | Indian team touchdown London; here are pictures
Sachin Tendulkar makes BIG statement on Indian team
Chris Gayle makes BIG statement on bowlers around the world
West Indies Windies Cricket ICC Cricket World Cup 2019 Sports Jason Holder Jersey

redpajama_set_name: RedPajamaCommonCrawl


text:
Family finds over 100 snakes in their home
by: FOX 5 Digital Team
Posted: Oct 22, 2014 / 02:37 PM PDT
Updated: Oct 22, 2014 / 02:37 PM PDT
REGINA, Saskatchewan — It's a home Indiana Jones wouldn't want to live in and it might give most of us an anxiety attack.
More than a 100 snakes were found inside a home, and they're all being rescued.
"We were finding them basically anywhere – in the cracks of the stone foundation, under the floor, under some boxes and other things – and we just were picking them up in pillow cases and putting them into a bucket just to count them," Megan Lawrence of Salthaven West Wildlife Rehab said on Tuesday.
Lawrence doesn't usually have to search for the animals she rescues, but the call she got this weekend was a bit unusual.
"The family contacted us when they found a few garter snakes in their basement, and then they started finding more and more. And then they were finding them in their kitchen and the bedrooms. And they decided it wasn't a good idea to have them there anymore."
Lawrence helps rehabilitate animals at the rehab center. She took stock of all the serpents they gathered and brought them back to shelter.
"And the final tally — we had 102 by the end of the weekend," she said.
"I get a feeling that's about a normal count, so 100 garter snakes in the basement would not surprise me," Ray Poulin of Royal Saskatchewan Museum said.
Poulin, who is knowledgeable about snakes, said it's far from unusual for people to find plains garter snakes that came inside. It's actually the most common call he receives when snakes are hunting for a warm place to spend the winter.
"Usually the snakes at this time are going down, right? So they're coming up to your house and going straight down one of the cracks in the soil around your house and finding a way in that way usually," Poulin said.
In this case, that's exactly what these garter snakes did at this house.
Now Lawrence has 102 new mouths to feed and take care of, possibly for the winter unless she can clear them for release in the wild before they start hibernation.
"They eat things in the wild such as fish, earthworms, frogs and small rodents and insects, so we could use things like donations toward earth worms and minnows to feed them for the winter if we're going to keep them," she said.
Border Report / 1 hour ago