last

[book_gpu.git] / BookGPU / Chapters / chapter16 / gpu.tex

\section{New parallel envelope-following method}
\label{sec:gmres}

In this section, we explain how to efficiently
use matrix-free GMRES to solve
the Newton update problems with implicit sensitivity calculation,
i.e., the steps enclosed by the double dashed block
in Figure~\ref{fig:ef_flow}.
Then implementation issues of GPU acceleration
will be discussed in detail. 
Finally,  the Gear-2 integration is briefly introduced.

\subsection{GMRES solver for Newton update equation}

\underline{G}eneralized \underline{M}inimum \underline{Res}idual,
or GMRES method is an iterative method for solving
systems of linear equations ($A x=b$) with dense matrix $A$.
The standard GMRES\index{iterative method!GMRES} is given in Algorithm~\ref{alg:GMRES}.
It constructs a Krylov subspace\index{iterative method!Krylov subspace} with order $m$,
\[ \mathcal{K}_m = \mathrm{span}( b, A^{} b, A^2 b,\ldots, A^{m-1} b ),\]
where the approximate solution $x_m$ resides.
In practice, an orthonormal basis $V_m$ that spans the
subspace $\mathcal{K}_{m}$ can be generated by
the Arnoldi iterations\index{iterative method!Arnoldi iterations}.
The goal of GMRES is to search for an optimal coefficient $y$
such that the linear combination $x_m = V_m y$ will minimize
its residual $\| b-Ax_m \|_2$.
The Arnoldi iteration also creates a by-product,
an upper Hessenberg matrix\index{Hessenberg matrix} $\tilde{H}_m \in R^{(m+1)\times m}$.
Thus, the projection of $A$ on the orthonormal basis $V_m$ is described by
the Arnoldi decomposition
$
 AV_m = V_{m+1} \tilde{H}_{m},
$
which is useful to check the residual at each iteration
without forming $x_m$,
and to solve for coefficient $y$ when residual is smaller than
a preset tolerance~\cite{Golub:Book'96}.

% Givens rotations are generated to triangularize
% the Hessenberg matrix $\tilde{H}_m$ in order to solve the overdetermined system.
% These same Givens rotations are also applied to the right-hand-side $h_{1,0}e_1$,
% whose last element provides residual information on whether the tolerance
% has been met in the current iteration.

%$\left\| h_{1,0}e_1-\tilde{H}_m y_m \right\|_2$.


%% \begin{algorithm}
%% \caption{Standard GMRES\index{GMRES} algorithm.} \label{alg:GMRES}
%% \begin{algorithmic}[1]
%%   \REQUIRE $ A \in \mathbb{R}^{N \times N}$, $b \in \mathbb{R}^N$,
%%       and initial guess $x_0 \in \mathbb{R}^N$
%%   \ENSURE $x \in \mathbb{R}^N$: $\| b - A x\|_2 < tol$

%%   \STATE $r = b - A x_0$
%%   \STATE $h_{1,0}=\left \| r \right \|_2$
%%   \STATE $m=0$

%%   \WHILE{$m < max\_iter$}
%%     \STATE $m = m+1$
%%     \STATE $v_{m} = r / h_{m,m-1}$
%%     \STATE \label{line:mvp} $r = A v_m$
%%     \FOR{$i = 1\ldots m$}
%%       \STATE $h_{i,m} = \langle v_i, r \rangle$
%%       \STATE $r = r - h_{i,m} v_i$
%%     \ENDFOR
%%     \STATE $h_{m+1,m} = \left\| r \right\|_2$\label{line:newnorm}
%%     %\STATE Generate Givens rotations to triangularize $\tilde{H}_m$
%%     %\STATE Apply Givens rotations on $h_{1,0}e_1$ to get residual $\epsilon$
%%     \STATE Compute the residual $\epsilon$
%%     \IF{$\epsilon < tol$}
%%     \STATE Solve the problem: minimize $\|b-Ax_m\|_2$
%%         \STATE Return $x_m = x_0 + V_m y_m$
%%     \ENDIF
%%   \ENDWHILE
%% \end{algorithmic}
%% \end{algorithm}

\begin{algorithm}
\caption{Standard GMRES\index{iterative method!GMRES} algorithm.} \label{alg:GMRES}
  \KwIn{ $ A \in \mathbb{R}^{N \times N}$, $b \in \mathbb{R}^N$,
      and initial guess $x_0 \in \mathbb{R}^N$}
  \KwOut{ $x \in \mathbb{R}^N$: $\| b - A x\|_2 < tol$}

  $r = b - A x_0$\;
  $h_{1,0}=\left \| r \right \|_2$\;
  $m=0$\;

  \While{$m < max\_iter$} {
    $m = m+1$;
    $v_{m} = r / h_{m,m-1}$\;
    \label{line:mvp} $r = A v_m$\;
    \For{$i = 1\ldots m$} {
      $h_{i,m} = \langle v_i, r \rangle$\;
      $r = r - h_{i,m} v_i$\;
    }
    $h_{m+1,m} = \left\| r \right\|_2$\label{line:newnorm} \;
    %\STATE Generate Givens rotations to triangularize $\tilde{H}_m$
    %\STATE Apply Givens rotations on $h_{1,0}e_1$ to get residual $\epsilon$
    Compute the residual $\epsilon$\;
    \If{$\epsilon < tol$} {
      Solve the problem: minimize $\|b-Ax_m\|_2$\;
      Return $x_m = x_0 + V_m y_m$\;
    }
  }
\end{algorithm}


At a first glance, the cost of using standard GMRES directly to
solve the Newton update in Eq.~\eqref{eq:Newton}
seems to come mainly from two parts: the
formulation of the sensitivity matrix $J = \ud x_M/\ud x_0$
in Eq.~\eqref{eq:sensM} in Section~\ref{sec:gear},
and the iteration inside the standard GMRES,
especially the matrix-vector multiplication and the orthonormal
basis construction (Line~\ref{line:mvp} through
Line~\ref{line:newnorm} in Algorithm~\ref{alg:GMRES}).
Based on the observation that only the matrix-vector product is
required in GMRES, the work in~\cite{Telichevesky:DAC'95} introduces
an efficient matrix-free algorithm in the shooting-Newton method,
where the equation solving part also involves with a sensitivity
matrix.
The matrix-free method does not take an explicit matrix as
input, but directly passes the saved capacitance matrices $C_i$ and
the LU factorizations of $J_i$, $i=0,\ldots,M$, into the Arnoldi
iteration for Krylov subspace generation.
Therefore, it avoids the
cost of forming the dense sensitivity matrix and focuses on
subspace construction.
Briefly speaking, Line~\ref{line:mvp} will be replaced
by a procedure without explict $A$, and we will talk
about the flow of matrix-free generation of new basis vectors
in later sections.


\subsection{Parallelization on GPU platforms}
\label{sec:gpu}
% In this subsection, we describe how to parallelize
% the proposed matrix-free GMRES solver with Gear-2 sensitivity
% on GPU platforms.
There exist many GPU-friendly computing operations
in GMRES, such as the vector addition (\verb|axpy|),
2-norm of vector (\verb|nrm2|), and
sparse matrix-vector multiplication (\verb|csrmv|),
which have been parallelized in the CUDA Basic Linear
Algebra Subroutine (CUBLAS) Library and the CUDA Sparse Linear Algebra
Library (CUSPARSE)~\cite{CUDASDK}.

GPU programming is typically limited by the data transfer bandwidth as
GPU favors computationally intensive algorithms~\cite{Kirk:Book'10}.
So how to efficiently transfer the data and wisely partition
the data between CPU memory and GPU memory
is crucial for GPU programming.
In the following, we discuss these two issues in our implementation.

As noted in Section~\ref{sec:ef}, the envelope-following method
requires the matrices gathered from all the time steps in a
period in order to solve a Newton update.
At each time step, SPICE\index{SPICE} has
to linearize device models, stamp matrix elements
into MNA (short for modified nodal analysis\index{modified nodal analysis, or MNA}) matrices,
and solve circuit equations in its inner Newton iteration\index{iterative method!Newton iteration}.
When convergence is attained,
circuit states are saved and then next time step begins.
This is also the time when we store the needed matrices
for the envelope-following computation.
Since these data are involved in the calculation of Gear-2 sensitivity
matrix in the generation of Krylov subspace vectors in
Algorithm~\ref{alg:mf_Gear}, it is desirable that all of these
matrices are transferred to GPU for its data parallel capability.

To efficiently transfer those large data, we explore asynchronous
memory copy between host and device in the recent GPUs (Fermi architecture),
so that the copying overlaps with the host's computing of
the next time step's circuit state. The implementation of asynchronous matrices
copy includes two parts: allocating page-locked memory, also known as
pinned memory, where we save matrices for one time step, and using
asynchronous memory transfer to copy these data to GPU memory.  While
it is known that page-locked host memory is a scarce resource and
should not be overused, the demand of memory size of the data
generated at one time step can be generously accommodated by today's
mainstream server configurations.

\begin{figure}[htbp]
  \centering
  \resizebox{.9\textwidth}{!}{\input{./Chapters/chapter16/figures/gmres_flow.pdf_t}}
  \caption{GPU parallel solver for envelope-following update.}
  \label{fig:gmres}
\end{figure}


%When SPICE simulation reaches the end of the period, all data needed by
%Newton update equation are ready in GPU memory, and the matrix-free
%GMRES can be applied directly on them to solve for the update vector
%in \eqref{eq:Newton}. 

% Knowing the dimensions or the sizes of the
%matrices in a computation helps to boost GPU computing efficiency.
The second issue is to decide the location of data between CPU and GPU memories.
Therefore let us first make a rough sketch of the quantities in the
GMRES Algorithm~\ref{alg:GMRES}.  Although GMRES tends to converge
quickly for most circuit examples, i.e., the iteration number $m \ll
N$, the space for storing all the subspace basis $V_m$ of $N$-by-$m$,
i.e., $m$ column vectors with $N$-length, is still big.  In addition,
every newly generated matrix-vector product needs to be orthogonalized
with respect to all its previous basis vectors in
the Arnoldi iterations\index{Arnoldi iterations}.
Hence, keeping all the vectors of $V_m$ in GPU global
memory allows GPU to handle those operations, such as inner-product of
basis vectors (\verb|dot|) and vector subtraction (\verb|axpy|), in
parallel.  

On the other hand, it is better to keep the Hessenberg matrix
$\tilde{H}$, where intermediate results of
the orthogonalization are stored, at the host side.
This comes with the following reasons.
First, its size is $(m+1)$-by-$m$ at most, rather small
if compared to circuit matrices and Krylov basis.
Besides, it is also necessary to triangularize $\tilde{H}$
and check the residual regularly in
each iteration so the GMRES can return the approximate solution as
soon as the residual is below a preset tolerance.
Hence, in consideration of the serial nature of the trianularization,
the small size of Hessenberg matrix,
and the frequent inspection of values by host, it is
preferable to allocate $\tilde{H}$ in CPU (host) memory.
As shown in Figure~\ref{fig:gmres}, the memory copy from device to host
is called each time when Arnoldi iteration generates a new vector
and the orthogonalization produces the vector $h$.

% \begin{figure}[!tb]
% \centering
% \resizebox{.45\textwidth}{!}{\input{./matlab/drawing.ps_tex}}
% \caption{The flow of envelope-following method.}
% \label{fig:gmres}
% %\end{center}
% \end{figure}

% talk about gear 2 method here.
\input{Chapters/chapter16/bdf.tex}