new

[GMRES2stage.git] / paper.tex

diff --git a/paper.tex b/paper.tex
index 62927056fdd97597284db09396876b30ce329db5..fe7fa3952ac90a99424b026e688b0cb5b2d13ac7 100644 (file)
--- a/paper.tex
+++ b/paper.tex
@@ -369,11 +369,8 @@
 %
 % paper title
 % can use linebreaks \\ within to get better formatting as desired
 %
 % paper title
 % can use linebreaks \\ within to get better formatting as desired

-\title{TSIRM: A Two-Stage Iteration with least-square Residual Minimization algorithm to solve large sparse linear systems}
-%où
-%\title{A two-stage algorithm with error minimization to solve large sparse linear systems}
-%où
-%\title{???}
+\title{TSIRM: A Two-Stage Iteration with least-squares Residual Minimization algorithm to solve large sparse linear systems}
+

 
 
 
 
 
 


@@ -428,16 +425,17 @@ Email: lilia.ziane@inria.fr}
 
 
 \begin{abstract}
 
 
 \begin{abstract}

-In  this article,  a  two-stage  iterative algorithm is proposed to improve  the
-convergence of Krylov based iterative methods,  typically those of GMRES variants. The
-principle of  the proposed approach  is to  build an external  iteration over  the Krylov
-method, and to  frequently store its current  residual   (at  each
-GMRES restart for instance). After a given number of outer iterations, a minimization
-step  is applied  on the  matrix composed by the  saved residuals,  in  order to
-compute a better solution and to make  new iterations if required.  It is proven that
-the proposal has  the same convergence properties than the  inner embedded method itself. 
-Experiments using up  to 16,394 cores also show that the proposed algorithm
-runs around 5 or 7 times faster than GMRES.
+In  this article, a  two-stage iterative  algorithm is  proposed to  improve the
+convergence  of  Krylov  based  iterative  methods,  typically  those  of  GMRES
+variants.  The  principle of  the  proposed approach  is  to  build an  external
+iteration over the  Krylov method, and to frequently  store its current residual
+(at each GMRES restart for instance).  After a given number of outer iterations,
+a least-squares minimization step is applied on the matrix composed by the saved
+residuals, in order  to compute a better solution and to  make new iterations if
+required.  It  is proven that the  proposal has the  same convergence properties
+than the  inner embedded  method itself.  Experiments  using up to  16,394 cores
+also  show that the  proposed algorithm  runs around  5 or  7 times  faster than
+GMRES.

 \end{abstract}
 
 \begin{IEEEkeywords}
 \end{abstract}
 
 \begin{IEEEkeywords}


@@ -584,7 +582,7 @@ performances.
 
 The present  article is  organized as follows.   Related works are  presented in
 Section~\ref{sec:02}. Section~\ref{sec:03} details the two-stage algorithm using
 
 The present  article is  organized as follows.   Related works are  presented in
 Section~\ref{sec:02}. Section~\ref{sec:03} details the two-stage algorithm using

-a  least-square  residual   minimization,  while  Section~\ref{sec:04}  provides
+a  least-squares  residual   minimization,  while  Section~\ref{sec:04}  provides

 convergence  results  regarding this  method.   Section~\ref{sec:05} shows  some
 experimental  results  obtained  on  large  clusters  using  routines  of  PETSc
 toolkit. This research work ends by  a conclusion section, in which the proposal
 convergence  results  regarding this  method.   Section~\ref{sec:05} shows  some
 experimental  results  obtained  on  large  clusters  using  routines  of  PETSc
 toolkit. This research work ends by  a conclusion section, in which the proposal


@@ -607,7 +605,7 @@ is summarized while intended perspectives are provided.
 
 %%%*********************************************************
 %%%*********************************************************
 
 %%%*********************************************************
 %%%*********************************************************

-\section{Two-stage algorithm with least-square residuals minimization}
+\section{Two-stage iteration with least-squares residuals minimization algorithm}

 \label{sec:03}
 A two-stage algorithm is proposed  to solve large  sparse linear systems  of the
 form  $Ax=b$,  where  $A\in\mathbb{R}^{n\times   n}$  is  a  sparse  and  square
 \label{sec:03}
 A two-stage algorithm is proposed  to solve large  sparse linear systems  of the
 form  $Ax=b$,  where  $A\in\mathbb{R}^{n\times   n}$  is  a  sparse  and  square


@@ -618,7 +616,7 @@ inner-outer iteration  solver based  on iterative Krylov  methods. The  main
 key-points of the proposed solver are given in Algorithm~\ref{algo:01}.
 It can be summarized as follows: the
 inner solver is a Krylov based one. In order to accelerate its convergence, the 
 key-points of the proposed solver are given in Algorithm~\ref{algo:01}.
 It can be summarized as follows: the
 inner solver is a Krylov based one. In order to accelerate its convergence, the 

-outer solver periodically applies a least-square minimization  on the residuals computed by  the inner one. %Tsolver which does not required to be changed.
+outer solver periodically applies a least-squares minimization  on the residuals computed by  the inner one. %Tsolver which does not required to be changed.

 
 At each outer iteration, the sparse linear system $Ax=b$ is partially 
 solved using only $m$
 
 At each outer iteration, the sparse linear system $Ax=b$ is partially 
 solved using only $m$


@@ -629,8 +627,8 @@ inner solver. The current approximation of the Krylov method is then stored insi
 $S$ composed by the successive solutions that are computed during inner iterations.
 
 At each $s$ iterations, the minimization step is applied in order to
 $S$ composed by the successive solutions that are computed during inner iterations.
 
 At each $s$ iterations, the minimization step is applied in order to

-compute a new  solution $x$. For that, the previous  residuals are computed with
-$(b-AS)$. The minimization of the residuals is obtained by 
+compute a new  solution $x$. For that, the previous  residuals of $Ax=b$ are computed by
+the inner iterations with $(b-AS)$. The minimization of the residuals is obtained by  

 \begin{equation}
    \underset{\alpha\in\mathbb{R}^{s}}{min}\|b-R\alpha\|_2
 \label{eq:01}
 \begin{equation}
    \underset{\alpha\in\mathbb{R}^{s}}{min}\|b-R\alpha\|_2
 \label{eq:01}


@@ -639,7 +637,7 @@ with $R=AS$. Then the new solution $x$ is computed with $x=S\alpha$.
 
 
 In  practice, $R$  is a  dense rectangular  matrix belonging in  $\mathbb{R}^{n\times s}$,
 
 
 In  practice, $R$  is a  dense rectangular  matrix belonging in  $\mathbb{R}^{n\times s}$,

-with $s\ll n$.   In order  to minimize~(\eqref{eq:01}), a  least-square method  such as
+with $s\ll n$.   In order  to minimize~\eqref{eq:01}, a  least-squares method  such as

 CGLS ~\cite{Hestenes52}  or LSQR~\cite{Paige82} is used. Remark that these  methods are more
 appropriate than a single direct method in a parallel context.
 
 CGLS ~\cite{Hestenes52}  or LSQR~\cite{Paige82} is used. Remark that these  methods are more
 appropriate than a single direct method in a parallel context.
 


@@ -657,7 +655,7 @@ appropriate than a single direct method in a parallel context.
     \State $S_{k \mod s}=x^k$ \label{algo:store}
     \If {$k \mod s=0$ {\bf and} error$>\epsilon_{kryl}$}
       \State $R=AS$ \Comment{compute dense matrix} \label{algo:matrix_mul}
     \State $S_{k \mod s}=x^k$ \label{algo:store}
     \If {$k \mod s=0$ {\bf and} error$>\epsilon_{kryl}$}
       \State $R=AS$ \Comment{compute dense matrix} \label{algo:matrix_mul}

-            \State Solve least-square problem $\underset{\alpha\in\mathbb{R}^{s}}{min}\|b-R\alpha\|_2$ \label{algo:}
+            \State $\alpha=Solve\_Least\_Squares(R,b,max\_iter_{ls})$ \label{algo:}

       \State $x^k=S\alpha$  \Comment{compute new solution}
     \EndIf
   \EndFor
       \State $x^k=S\alpha$  \Comment{compute new solution}
     \EndIf
   \EndFor


@@ -675,7 +673,7 @@ $\epsilon_{tsirm}$).  Line~\ref{algo:store}, $S_{k~ mod~ s}=x^k$ consists in cop
 solution  $x_k$  into the  column  $k~ mod~ s$ of  the  matrix  $S$. After  the
 minimization, the matrix $S$ is reused with the new values of the residuals.  To
 solve the minimization problem, an  iterative method is used. Two parameters are
 solution  $x_k$  into the  column  $k~ mod~ s$ of  the  matrix  $S$. After  the
 minimization, the matrix $S$ is reused with the new values of the residuals.  To
 solve the minimization problem, an  iterative method is used. Two parameters are

-required for that: the maximum number of iteration and the threshold to stop the
+required for that: the maximum number of iterations and the threshold to stop the

 method.
 
 Let us summarize the most important parameters of TSIRM:
 method.
 
 Let us summarize the most important parameters of TSIRM:


@@ -683,22 +681,22 @@ Let us summarize the most important parameters of TSIRM:
 \item $\epsilon_{tsirm}$: the threshold to stop the TSIRM method;
 \item $max\_iter_{kryl}$: the maximum number of iterations for the Krylov method;
 \item $s$: the number of outer iterations before applying the minimization step;
 \item $\epsilon_{tsirm}$: the threshold to stop the TSIRM method;
 \item $max\_iter_{kryl}$: the maximum number of iterations for the Krylov method;
 \item $s$: the number of outer iterations before applying the minimization step;

-\item $max\_iter_{ls}$: the maximum number of iterations for the iterative least-square method;
-\item $\epsilon_{ls}$: the threshold used to stop the least-square method.
+\item $max\_iter_{ls}$: the maximum number of iterations for the iterative least-squares method;
+\item $\epsilon_{ls}$: the threshold used to stop the least-squares method.

 \end{itemize}
 
 
 The  parallelisation  of  TSIRM  relies   on  the  parallelization  of  all  its
 \end{itemize}
 
 
 The  parallelisation  of  TSIRM  relies   on  the  parallelization  of  all  its

-parts. More  precisely, except  the least-square step,  all the other  parts are
+parts. More  precisely, except  the least-squares step,  all the other  parts are

 obvious to  achieve out in parallel. In  order to develop a  parallel version of
 our   code,   we   have   chosen  to   use   PETSc~\cite{petsc-web-page}.    For
 line~\ref{algo:matrix_mul} the  matrix-matrix multiplication is  implemented and
 efficient since the  matrix $A$ is sparse and since the  matrix $S$ contains few
 colums in  practice. As explained  previously, at least  two methods seem  to be
 obvious to  achieve out in parallel. In  order to develop a  parallel version of
 our   code,   we   have   chosen  to   use   PETSc~\cite{petsc-web-page}.    For
 line~\ref{algo:matrix_mul} the  matrix-matrix multiplication is  implemented and
 efficient since the  matrix $A$ is sparse and since the  matrix $S$ contains few
 colums in  practice. As explained  previously, at least  two methods seem  to be

-interesting to solve the least-square minimization, CGLS and LSQR.
+interesting to solve the least-squares minimization, CGLS and LSQR.

 
 In the following  we remind the CGLS algorithm. The LSQR  method follows more or
 
 In the following  we remind the CGLS algorithm. The LSQR  method follows more or

-less the same principle but it take more place, so we briefly explain the parallelization of CGLS which is similar to LSQR.
+less the same principle but it takes more place, so we briefly explain the parallelization of CGLS which is similar to LSQR.

 
 \begin{algorithm}[t]
 \caption{CGLS}
 
 \begin{algorithm}[t]
 \caption{CGLS}


@@ -725,7 +723,7 @@ less the same principle but it take more place, so we briefly explain the parall
 
 
 In each iteration  of CGLS, there is two  matrix-vector multiplications and some
 
 
 In each iteration  of CGLS, there is two  matrix-vector multiplications and some

-classical operations:  dots, norm, multiplication  and addition on  vectors. All
+classical operations:  dot product, norm, multiplication  and addition on  vectors. All

 these operations are easy to implement in PETSc or similar environment.
 
 
 these operations are easy to implement in PETSc or similar environment.
 
 


@@ -757,7 +755,7 @@ In order to see the influence of our algorithm with only one processor, we first
 show  a comparison  with the  standard version  of GMRES  and our  algorithm. In
 Table~\ref{tab:01},  we  show  the  matrices  we  have used  and  some  of  them
 characteristics. For all  the matrices, the name, the field,  the number of rows
 show  a comparison  with the  standard version  of GMRES  and our  algorithm. In
 Table~\ref{tab:01},  we  show  the  matrices  we  have used  and  some  of  them
 characteristics. For all  the matrices, the name, the field,  the number of rows

-and the number of nonzero elements is given.
+and the number of nonzero elements are given.

 
 \begin{table}[htbp]
 \begin{center}
 
 \begin{table}[htbp]
 \begin{center}


@@ -780,7 +778,7 @@ torso3             & 2D/3D problem & 259,156 & 4,429,042 \\
 
 The following  parameters have been chosen  for our experiments.   As by default
 the restart  of GMRES is performed every  30 iterations, we have  chosen to stop
 
 The following  parameters have been chosen  for our experiments.   As by default
 the restart  of GMRES is performed every  30 iterations, we have  chosen to stop

-the GMRES every 30 iterations, $max\_iter_{kryl}=30$).  $s$ is set to 8. CGLS is
+the GMRES every 30 iterations (\emph{i.e.} $max\_iter_{kryl}=30$).  $s$ is set to 8. CGLS is

 chosen  to minimize  the least-squares  problem with  the  following parameters:
 $\epsilon_{ls}=1e-40$ and $max\_iter_{ls}=20$.  The external precision is set to
 $\epsilon_{tsirm}=1e-10$.  Those  experiments have been performed  on a Intel(R)
 chosen  to minimize  the least-squares  problem with  the  following parameters:
 $\epsilon_{ls}=1e-40$ and $max\_iter_{ls}=20$.  The external precision is set to
 $\epsilon_{tsirm}=1e-10$.  Those  experiments have been performed  on a Intel(R)


@@ -792,7 +790,7 @@ systems obtained with the previous matrices  with a GMRES variant and with out 2
 stage algorithm are  given. In the second column, it can  be noticed that either
 gmres or fgmres is used to  solve the linear system.  According to the matrices,
 different  preconditioner is used.   With TSIRM,  the same  solver and  the same
 stage algorithm are  given. In the second column, it can  be noticed that either
 gmres or fgmres is used to  solve the linear system.  According to the matrices,
 different  preconditioner is used.   With TSIRM,  the same  solver and  the same

-preconditionner is used.  This Table shows that TSIRM can drastically reduce the
+preconditionner are used.  This Table shows that TSIRM can drastically reduce the

 number of iterations to reach the  convergence when the number of iterations for
 the normal GMRES is more or less  greater than 500. In fact this also depends on
 tow  parameters: the  number  of iterations  to  stop GMRES  and  the number  of
 number of iterations to reach the  convergence when the number of iterations for
 the normal GMRES is more or less  greater than 500. In fact this also depends on
 tow  parameters: the  number  of iterations  to  stop GMRES  and  the number  of


@@ -826,12 +824,12 @@ torso3             & fgmres / sor  & 37.70 & 565 & 34.97 & 510 \\
 
 
 
 
 
 

-In order to perform larger  experiments, we have tested some example application
+In order to perform larger  experiments, we have tested some example applications

 of PETSc. Those  applications are available in the ksp part  which is suited for
 scalable linear equations solvers:
 \begin{itemize}
 \item ex15  is an example  which solves in  parallel an operator using  a finite
 of PETSc. Those  applications are available in the ksp part  which is suited for
 scalable linear equations solvers:
 \begin{itemize}
 \item ex15  is an example  which solves in  parallel an operator using  a finite

-  difference  scheme.   The  diagonal  is  equals to  4  and  4  extra-diagonals
+  difference  scheme.   The  diagonal  is  equal to  4  and  4  extra-diagonals

   representing the neighbors in each directions  is equal to -1. This example is
   used  in many  physical phenomena, for  example, heat  and fluid  flow, wave
   propagation...
   representing the neighbors in each directions  is equal to -1. This example is
   used  in many  physical phenomena, for  example, heat  and fluid  flow, wave
   propagation...


@@ -900,7 +898,7 @@ corresponds to 30*12, there are $max\_iter_{ls}$ which corresponds to 15.
 \begin{figure}[htbp]
 \centering
   \includegraphics[width=0.45\textwidth]{nb_iter_sec_ex15_juqueen}
 \begin{figure}[htbp]
 \centering
   \includegraphics[width=0.45\textwidth]{nb_iter_sec_ex15_juqueen}

-\caption{Number of iterations per second with ex15 and the same parameters than in Table~\ref{tab:03}}
+\caption{Number of iterations per second with ex15 and the same parameters than in Table~\ref{tab:03} (weak scaling)}

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


@@ -968,6 +966,13 @@ In Table~\ref{tab:04}, some experiments with example ex54 on the Curie architect
 \end{center}
 \end{table*}
 
 \end{center}
 \end{table*}
 

+\begin{figure}[htbp]
+\centering
+  \includegraphics[width=0.45\textwidth]{nb_iter_sec_ex54_curie}
+\caption{Number of iterations per second with ex54 and the same parameters than in Table~\ref{tab:05} (strong scaling)}
+\label{fig:02}
+\end{figure}
+

 %%%*********************************************************
 %%%*********************************************************
 
 %%%*********************************************************
 %%%*********************************************************