\subsection{Gear-2 based sensitivity calculation}
\label{sec:gear}

The Gear-2 integration method is a backward differentiation
formula (BDF), which approximates the derivative of a function
using information from past few steps.  Gear-2 method has been
shown to be more suitable for many practical problems such as
stiff problems where circuit behavior is affected by different time constants
(fast ones with large poles and slow ones with small poles)~\cite{Vlach:Book'94}.
Switching power converters and RF systems are typically stiff systems
as waveforms changing in two different time scales are involved.

% Although on one time step, Gear-2
% involves more computation than backward-Euler and trapezoidal method,
% it has better convergence and stability, and usually allows larger
% time step.  As a result, fewer time points are needed for a given
% transient sweep range, and the simulation is speeded up with the same
% accuracy

Given the DAE~(\ref{eq:dae}), at the $i$-th time step,
Gear-2 integration\index{Gear-2 integration}
approximates the derivative $\ud q(x(t))/\ud t$
by a two-step backward finite difference,
\[
   \frac{\ud q(x(t_i))}{\ud t}
   %\ud q(x(t_i)) / \ud t
   = \frac{1}{h_i}
     \left[ q(x(t_{i})) - \alpha_1^i q(x(t_{i-1})) - \alpha_2^i q(x(t_{i-2}))\right],
\]
where the coefficients $\alpha$'s are chosen to satisfy Gear's
backward differentiation formula~\cite{Gear:book'71}.
For uniformly discretised time steps,
the index $i$ in $h_i$, $\alpha_1^i$ and $\alpha_2^i$ can be omitted.

For the first step $t_1$,
only the initial condition at $t_0$ is available.
Therefore backward Euler is used, i.e.,
\begin{equation}\label{eq:BE}
   \tfrac{1}{h_1}[q(x_1)-q(x_0)] + f(x_1) = b_1.
\end{equation}
Since $x_0$ directly affects the solution of $x_1$,
the sensitivity matrix\index{sensitivity matrix} up to now is
\begin{equation}\label{eq:sens1}
   \frac{\ud x_1}{\ud x_0}
   = \left[\frac{1}{h_1}C_1+G_1\right]^{-1} \frac{C_0}{h_1}
   = J_1^{-1} \frac{C_0}{h_1}.
\end{equation}
%For the sake of simplicity,
Let $J_i$ denote $(1/h_i)C_i + G_i$
in the remaining part of this chapter.
In addition, the LU\index{LU} factorizations of $J_i$ are already calculated
since they are used to solve the circuit state at
each time step before we calculate the sensitivity.

Next, for the solution at $t_2$,
with the previous two steps information available,
Gear-2 integration can be applied,
\begin{equation}\label{eq:Gear_t2}
   \tfrac{1}{h_2}[q(x_2) - \alpha_1 q(x_1) - \alpha_2 q(x_0)] + f(x_2) = b_2.
\end{equation}
In view of the sensitivity of $x_2$ with respect to the changes of $x_0$,
\eqref{eq:Gear_t2} indicates that $x_2$'s perturbation can be
traced back to $x_0$ along two paths:
$x_2$ is directly affected by $x_0$,
and $x_2$ is also affected indirectly by $x_0$ via $x_1$.
That is,
\[ %begin{equation}
  \frac{\ud x_2}{\ud x_0} = \frac{\partial x_2}{\partial x_1}\frac{\ud x_1}{\ud x_0}
                          + \frac{\partial x_2}{\partial x_0},
\] %end{equation}
where the two partial derivatives are
\[
   \frac{\partial x_2}{\partial x_1} = J_2^{-1} \frac{\alpha_1}{h_2} C_1, \qquad
   \frac{\partial x_2}{\partial x_0} = J_2^{-1} \frac{\alpha_2}{h_2} C_0,
\]
and $\ud x_1/\ud x_0$ %, the sensitivity of $x_1$ with respect to $x_0$,
is calculated previously in \eqref{eq:sens1}.
Thus, the sensitivity matrix deduced from \eqref{eq:Gear_t2} is
\begin{equation}\label{eq:sens2}
   \frac{\ud x_2}{\ud x_0}
   = J_2^{-1} \frac{\alpha_1}{h_2} C_1 \cdot J_1^{-1} \frac{C_0}{h_1}
   + J_2^{-1} \frac{\alpha_2}{h_2} C_0.
\end{equation}

Likewise, for the third time step $t_3$,
apply the Gear-2 integration formula to DAE~(\ref{eq:dae}),
\begin{equation}\label{eq:Gear_t3}
   \tfrac{1}{h_3} [q(x_3) - \alpha_1 q(x_2) - \alpha_2 q(x_1)] + f(x_3) = b_3,
\end{equation}
and the chain rule for sensitivity is
\[ %begin{equation}\label{eq:sens3}
% \begin{split}
%    \frac{\ud x_3}{\ud x_0} &= \frac{\partial x_3}{\partial x_2}\frac{\ud x_2}{\ud x_0}
%                             + \frac{\partial x_3}{\partial x_1}\frac{\ud x_1}{\ud x_0} \\
%    &=  J_3^{-1} \frac{\alpha_1}{h_3} C_2 \frac{\ud x_2}{\ud x_0}
%      + J_3^{-1} \frac{\alpha_2}{h_3} C_1 \frac{\ud x_1}{\ud x_0}.
% \end{split}
   \frac{\ud x_3}{\ud x_0} \!=\! \frac{\partial x_3}{\partial x_2}\frac{\ud x_2}{\ud x_0}
                            + \frac{\partial x_3}{\partial x_1}\frac{\ud x_1}{\ud x_0}
   \!=\!  J_3^{-1} \!\left[ \frac{\alpha_1}{h_3} C_2 \frac{\ud x_2}{\ud x_0}
                   \!+\!\frac{\alpha_2}{h_3} C_1 \frac{\ud x_1}{\ud x_0}\right],
\] %end{equation}
where both $\ud x_2/\ud x_0$ and $\ud x_1/\ud x_0$
are ready from the previous two time steps, i.e., Eqs.~(\ref{eq:sens2}) and~(\ref{eq:sens1}).
Follow this chain rule in the aforementioned recursive style,
the sensitivity matrix for Gear-2 integration is computed
along the remaining time steps up to the $M$-th step,
\begin{equation}\label{eq:sensM}
% \begin{split}
%    J = \frac{\ud x_M}{\ud x_0}
%      & =  J_M^{-1} \frac{\alpha_1}{h_M} C_{M-1} \frac{\ud x_{M-1}}{\ud x_0}  \\
%      &\quad + J_M^{-1} \frac{\alpha_2}{h_M} C_{M-2} \frac{\ud x_{M-2}}{\ud x_0}.
% \end{split}
J\! =\! \frac{\ud x_M}{\ud x_0}
 \! =\!  J_M^{-1} \! \left[ \frac{\alpha_1}{h_M} C_{M-1} \frac{\ud x_{M-1}}{\ud x_0}
     \! + \! \frac{\alpha_2}{h_M} C_{M-2} \frac{\ud x_{M-2}}{\ud x_0}\right].
\end{equation}

\begin{algorithm}
\caption{the matrix-free\index{matrix-free} method for
 Krylov subspace\index{iterative method!Krylov subspace} construction}
\label{alg:mf_Gear}
  \KwIn{ current Krylov subspace basis vector $v$,
           time step lengths $h_i$,
           saved $C_i$ matrices and LU\index{LU} factors of $J_i$, $i=0,\ldots,M$}
  \KwOut{ matrix-vector product $r$, such that $r = A v$}
  solve $J_1 p_2 = h_1^{-1} C_0 v$ for $p_2$\;
  solve $J_2 p_1 = h_2^{-1} \left( \alpha_1 C_1 p_2 + \alpha_2 C_0 v \right)$ for $p_1$\;
  \For{$i=3 \ldots M$}{
     solve
       $J_i p_0 = \alpha_1 h_i^{-1} C_{i-1} p_1 + \alpha_2 h_i^{-1} C_{i-2} p_2$
       for $p_0$ \label{line:mf_Gear_loop}\;
     $p_2 = p_1$\;
     $p_1 = p_0$\;
  }
   $r = (k-1) p_0 - k v$ \label{line:shift}\;
\end{algorithm}


%% \begin{algorithm}
%% \caption{The matrix-free\index{matrix-free} method for
%%  Krylov subspace\index{Krylov subspace} construction.}
%% \label{alg:mf_Gear}
%% \begin{algorithmic}[1]
%%   \REQUIRE current Krylov subspace basis vector $v$,
%%            time step lengths $h_i$,
%%            saved $C_i$ matrices and LU\index{LU} factors of $J_i$, $i=0,\ldots,M$
%%   \ENSURE matrix-vector product $r$, such that $r = A v$
%%   \STATE solve $J_1 p_2 = h_1^{-1} C_0 v$ for $p_2$
%%   \STATE solve $J_2 p_1 = h_2^{-1} \left( \alpha_1 C_1 p_2 + \alpha_2 C_0 v \right)$ for $p_1$
%%   \FOR{$i=3 \ldots M$}
%%     \STATE solve
%%        $J_i p_0 = \alpha_1 h_i^{-1} C_{i-1} p_1 + \alpha_2 h_i^{-1} C_{i-2} p_2$
%%        for $p_0$ \label{line:mf_Gear_loop}
%%     \STATE $p_2 = p_1$
%%     \STATE $p_1 = p_0$
%%   \ENDFOR
%%   \STATE $r = (k-1) p_0 - k v$ \label{line:shift}
%% \end{algorithmic}
%% \end{algorithm}

%It is easy to judge that, as $J$ is dense, explicitly
%computing $J$ is expensive and impractical for large problems.
%In addition, with a dense matrix,
%both direct or iterative solvers will be inefficient.
%Fortunately, one iterative solver, GMRES,
%only requires $J x$ (actually $Ax$),
%instead of an explicit $J$, to fulfill the solving process.
%This virtue greatly saves computation cost.
Note that as the matrix-free\index{matrix-free} GMRES method is applied,
which only requires matrix-vector multiplication and no explicit $J$ is required,
as explained in Algorithm~\ref{alg:mf_Gear}.
% Next section will explain how to use this technique for
% our proposed Gear-2 envelope-following.

With matrix-free method, the matrix-vector multiplication
in Line~\ref{line:mvp} of Algorithm~\ref{alg:GMRES}
is replaced by the iteration shown in Algorithm~\ref{alg:mf_Gear},
which calculates the multiplication product of the Gear-2 sensitivity we
encounter in envelope-following and a basis vector in the Krylov
subspace.
Note that the scaling and shift of $J$ in $A=(k-1)J-kI$,
as described in Eq.~\eqref{eq:A},
is taken into consideration by Line~\ref{line:shift}
of Algorithm~\ref{alg:mf_Gear}.

For the matrix-free generation of new basis vectors, it is
straightforward to apply CUBLAS and CUSPARSE routines, and some
customized GPU kernel functions to implement
Algorithm~\ref{alg:mf_Gear} with the stored LU matrices of $J_i$ and
$C_i$ mentioned earlier. For example in Line~\ref{line:mf_Gear_loop},
as the iteration index $i$ traverses all the $M$ time points'
matrix information,
sparse matrix-vector multiplication \verb|csrmv|
is first called to calculate $C_{i-1} p_1$ and $C_{i-2} p_2$.
And after the two resulted vectors are combined by \verb|axpy| of CUBLAS,
$p_0$ is solved for by \verb|trsv|, since as we noted before,
$J_i$ is already in LU factorization form when transient simulation
at step $i$ finished.
