1 %%%%%%%%%%%%%%%%%%%%%%
\r
2 \documentclass{doublecol-new}
\r
3 %%%%%%%%%%%%%%%%%%%%%%
\r
5 \usepackage{natbib,stfloats}
\r
6 \usepackage{mathrsfs}
\r
7 \usepackage[utf8]{inputenc}
\r
8 \usepackage[T1]{fontenc}
\r
9 \usepackage{algorithm}
\r
10 \usepackage{algpseudocode}
\r
11 \usepackage{amsmath}
\r
12 \usepackage{amssymb}
\r
13 \usepackage{multirow}
\r
14 \usepackage{graphicx}
\r
19 \def\newblock{\hskip .11em plus .33em minus .07em}
\r
22 \newtheorem{lemma}{Lemma}
\r
23 \newtheorem{theorem}[lemma]{Theorem}
\r
24 \newtheorem{corrolary}[lemma]{Corrolary}
\r
25 \newtheorem{conjecture}[lemma]{Conjecture}
\r
26 \newtheorem{proposition}[lemma]{Proposition}
\r
27 \newtheorem{claim}[lemma]{Claim}
\r
28 \newtheorem{stheorem}[lemma]{Wrong Theorem}
\r
29 %\newtheorem{algorithm}{Algorithm}
\r
32 \theoremstyle{THrm}{
\r
33 \newtheorem{definition}{Definition}[section]
\r
34 \newtheorem{question}{Question}[section]
\r
35 \newtheorem{remark}{Remark}
\r
36 \newtheorem{scheme}{Scheme}
\r
39 \theoremstyle{THhit}{
\r
40 \newtheorem{case}{Case}[section]
\r
42 \algnewcommand\algorithmicinput{\textbf{Input:}}
\r
43 \algnewcommand\Input{\item[\algorithmicinput]}
\r
45 \algnewcommand\algorithmicoutput{\textbf{Output:}}
\r
46 \algnewcommand\Output{\item[\algorithmicoutput]}
\r
52 \def\theequation{\arabic{equation}}
\r
54 \JOURNALNAME{\TEN{\it International Journal of High Performance Computing and Networking}}
\r
58 %\thispagestyle{empty}%
\r
60 %\NINE\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lcr@{}}%
\r
62 %Copyright \copyright\ 2012 Inderscience Enterprises Ltd. & &%
\r
74 \setcounter{page}{1}
\r
76 \LRH{R. Couturier, L. Ziane Khodja and C. Guyeux}
\r
78 \RRH{TSIRM: A Two-Stage Iteration with least-squares Residual Minimization algorithm}
\r
92 \title{TSIRM: A Two-Stage Iteration with least-squares Residual Minimization algorithm to solve large sparse linear and non linear systems}
\r
95 \authorA{Rapha\"el Couturier}
\r
97 \affA{Femto-ST Institute, University of Bourgogne Franche-Comte, France\\
\r
98 E-mail: raphael.couturier@univ-fcomte.fr}
\r
101 \authorB{Lilia Ziane Khodja}
\r
102 \affB{LTAS-Mécanique numérique non linéaire, University of Liege, Belgium \\
\r
103 E-mail: l.zianekhodja@ulg.ac.be}
\r
105 \authorC{Christophe Guyeux}
\r
106 \affC{Femto-ST Institute, University of Bourgogne Franche-Comte, France\\
\r
107 E-mail: christophe.guyeux@univ-fcomte.fr}
\r
111 In this paper, a two-stage iterative algorithm is proposed to improve the
\r
112 convergence of Krylov based iterative methods, typically those of GMRES
\r
113 variants. The principle of the proposed approach is to build an external
\r
114 iteration over the Krylov method, and to frequently store its current residual
\r
115 (at each GMRES restart for instance). After a given number of outer iterations,
\r
116 a least-squares minimization step is applied on the matrix composed by the saved
\r
117 residuals, in order to compute a better solution and to make new iterations if
\r
118 required. It is proven that the proposal has the same convergence properties
\r
119 than the inner embedded method itself.
\r
121 Several experiments have been performed
\r
122 with the PETSc solver with linear and nonlinear problems. They show good
\r
123 speedups compared to GMRES with up to 16,394 cores with different
\r
130 \KEYWORD{Iterative Krylov methods; sparse linear and non linear systems; two stage iteration; least-squares residual minimization; PETSc.}
\r
132 %\REF{to this paper should be made as follows: Rodr\'{\i}guez
\r
133 %Bol\'{\i}var, M.P. and Sen\'{e}s Garc\'{\i}a, B. (xxxx) `The
\r
134 %corporate environmental disclosures on the internet: the case of
\r
135 %IBEX 35 Spanish companies', {\it International Journal of Metadata,
\r
136 %Semantics and Ontologies}, Vol. x, No. x, pp.xxx\textendash xxx.}
\r
139 Raphaël Couturier ....
\r
141 \noindent Lilia Ziane Khodja ...
\r
143 \noindent Christophe Guyeux ...
\r
150 \section{Introduction}
\r
152 Iterative methods have recently become more attractive than direct ones to solve
\r
153 very large sparse linear systems~\cite{Saad2003}. They are more efficient in a
\r
154 parallel context, supporting thousands of cores, and they require less memory
\r
155 and arithmetic operations than direct methods~\cite{bahicontascoutu}. This is
\r
156 why new iterative methods are frequently proposed or adapted by researchers, and
\r
157 the increasing need to solve very large sparse linear systems has triggered the
\r
158 development of such efficient iterative techniques suitable for parallel
\r
161 Most of the successful iterative methods currently available are based on
\r
162 so-called ``Krylov subspaces''. They consist in forming a basis of successive
\r
163 matrix powers multiplied by an initial vector, which can be for instance the
\r
164 residual. These methods use vectors orthogonality of the Krylov subspace basis
\r
165 in order to solve linear systems. The best known iterative Krylov subspace
\r
166 methods are conjugate gradient and GMRES ones (Generalized Minimal RESidual).
\r
169 However, iterative methods suffer from scalability problems on parallel
\r
170 computing platforms with many processors, due to their need of reduction
\r
171 operations, and to collective communications to achieve matrix-vector
\r
172 multiplications. The communications on large clusters with thousands of cores
\r
173 and large sizes of messages can significantly affect the performances of these
\r
174 iterative methods. As a consequence, Krylov subspace iteration methods are often
\r
175 used with preconditioners in practice, to increase their convergence and
\r
176 accelerate their performances. However, most of the good preconditioners are
\r
177 not scalable on large clusters.
\r
179 In this research work, a two-stage algorithm based on two nested iterations
\r
180 called inner-outer iterations is proposed. This algorithm consists in solving
\r
181 the sparse linear system iteratively with a small number of inner iterations,
\r
182 and restarting the outer step with a new solution minimizing some error
\r
183 functions over some previous residuals. For further information on two-stage
\r
184 iteration methods, interested readers are invited to
\r
185 consult~\cite{Nichols:1973:CTS}. Two-stage algorithms are easy to parallelize on
\r
186 large clusters. Furthermore, the least-squares minimization technique improves
\r
187 its convergence and performances.
\r
189 The present article is organized as follows. Related works are presented in
\r
190 Section~\ref{sec:02}. Section~\ref{sec:03} details the two-stage algorithm using
\r
191 a least-squares residual minimization, while Section~\ref{sec:04} provides
\r
192 convergence results regarding this method. Section~\ref{sec:05} shows some
\r
193 experimental results obtained on large clusters using routines of PETSc
\r
194 toolkit. This research work ends by a conclusion section, in which the proposal
\r
195 is summarized while intended perspectives are provided.
\r
199 %%%*********************************************************
\r
200 %%%*********************************************************
\r
204 %%%*********************************************************
\r
205 %%%*********************************************************
\r
206 \section{Related works}
\r
208 Krylov subspace iteration methods have increasingly become key
\r
209 techniques for solving linear and nonlinear systems, or eigenvalue problems,
\r
210 especially since the increasing development of
\r
211 preconditioners~\cite{Saad2003,Meijerink77}. One reason for the popularity of
\r
212 these methods is their generality, simplicity, and efficiency to solve systems of
\r
213 equations arising from very large and complex problems.
\r
215 GMRES is one of the most widely used Krylov iterative method for solving sparse
\r
216 and large linear systems. It has been developed by Saad \emph{et
\r
217 al.}~\cite{Saad86} as a generalized method to deal with unsymmetric and
\r
218 non-Hermitian problems, and indefinite symmetric problems too. In its original
\r
219 version called full GMRES, this algorithm minimizes the residual over the
\r
220 current Krylov subspace until convergence in at most $n$ iterations, where $n$
\r
221 is the size of the sparse matrix. Full GMRES is however too expensive in the
\r
222 case of large matrices, since the required orthogonalization process per
\r
223 iteration grows quadratically with the number of iterations. For that reason,
\r
224 GMRES is restarted in practice after each $m\ll n$ iterations, to avoid the
\r
225 storage of a large orthonormal basis. However, the convergence behavior of the
\r
226 restarted GMRES, called GMRES($m$), in many cases depends quite critically on
\r
227 the $m$ value~\cite{Huang89}. Therefore in most cases, a preconditioning
\r
228 technique is applied to the restarted GMRES method in order to improve its
\r
231 To enhance the robustness of Krylov iterative solvers, some techniques have been
\r
232 proposed allowing the use of different preconditioners, if necessary, within the
\r
233 iteration itself instead of restarting. Those techniques may lead to
\r
234 considerable savings in CPU time and memory requirements. Van der Vorst
\r
235 in~\cite{Vorst94} has for instance proposed variants of the GMRES algorithm in
\r
236 which a different preconditioner is applied in each iteration, leading to the
\r
237 so-called GMRESR family of nested methods. In fact, the GMRES method is
\r
238 effectively preconditioned with other iterative schemes (or GMRES itself), where
\r
239 the iterations of the GMRES method are called outer iterations while the
\r
240 iterations of the preconditioning process is referred to as inner iterations.
\r
241 Saad in~\cite{Saad:1993} has proposed Flexible GMRES (FGMRES) which is another
\r
242 variant of the GMRES algorithm using a variable preconditioner. In FGMRES the
\r
243 search directions are preconditioned whereas in GMRESR the residuals are
\r
244 preconditioned. However, in practice, good preconditioners are those based on
\r
245 direct methods, as ILU preconditioners, which are not easy to parallelize and
\r
246 suffer from the scalability problems on large clusters of thousands of cores.
\r
248 Recently, communication-avoiding methods have been developed to reduce the
\r
249 communication overheads in Krylov subspace iterative solvers. On modern computer
\r
250 architectures, communications between processors are much slower than
\r
251 floating-point arithmetic operations on a given
\r
252 processor. Communication-avoiding techniques reduce either communications
\r
253 between processors or data movements between levels of the memory hierarchy, by
\r
254 reformulating the communication-bound kernels (more frequently SpMV kernels) and
\r
255 the orthogonalization operations within the Krylov iterative solver. Different
\r
256 works have studied the communication-avoiding techniques for the GMRES method,
\r
257 so-called CA-GMRES, on multicore processors and multi-GPU
\r
258 machines~\cite{Mohiyuddin2009,Hoemmen2010,Yamazaki2014}.
\r
260 Compared to all these works and to all the other works on Krylov iterative
\r
261 methods, the originality of our work is to build a second iteration over a
\r
262 Krylov iterative method and to minimize the residuals with a least-squares
\r
263 method after a given number of outer iterations.
\r
265 %%%*********************************************************
\r
266 %%%*********************************************************
\r
270 %%%*********************************************************
\r
271 %%%*********************************************************
\r
272 \section{TSIRM: Two-stage iteration with least-squares residuals minimization algorithm}
\r
274 A two-stage algorithm is proposed to solve large sparse linear systems of the
\r
275 form $Ax=b$, where $A\in\mathbb{R}^{n\times n}$ is a sparse and square
\r
276 nonsingular matrix, $x\in\mathbb{R}^n$ is the solution vector, and
\r
277 $b\in\mathbb{R}^n$ is the right-hand side. As explained previously, the
\r
278 algorithm is implemented as an inner-outer iteration solver based on iterative
\r
279 Krylov methods. The main key-points of the proposed solver are given in
\r
280 Algorithm~\ref{algo:01}. It can be summarized as follows: the inner solver is a
\r
281 Krylov based one. In order to accelerate its convergence, the outer solver
\r
282 periodically applies a least-squares minimization on the residuals computed by
\r
285 At each outer iteration, the sparse linear system $Ax=b$ is partially solved
\r
286 using only $m$ iterations of an iterative method, this latter being initialized
\r
287 with the last obtained approximation. The GMRES method~\cite{Saad86}, or any of
\r
288 its variants, can potentially be used as inner solver. The current approximation
\r
289 of the Krylov method is then stored inside a $n \times s$ matrix $S$, which is
\r
290 composed by the $s$ last solutions that have been computed during the inner
\r
291 iterations phase. In the remainder, the $i$-th column vector of $S$ will be
\r
294 At each $s$ iterations, another kind of minimization step is applied in order to
\r
295 compute a new solution $x$. For that, the previous residuals of $Ax=b$ are
\r
296 computed by the inner iterations with $(b-AS)$. The minimization of the
\r
297 residuals is obtained by
\r
299 \underset{\alpha\in\mathbb{R}^{s}}{min}\|b-R\alpha\|_2
\r
302 with $R=AS$. The new solution $x$ is then computed with $x=S\alpha$.
\r
305 In practice, $R$ is a dense rectangular matrix belonging in $\mathbb{R}^{n\times
\r
306 s}$, with $s\ll n$. In order to minimize~\eqref{eq:01}, a least-squares
\r
307 method such as CGLS ~\cite{Hestenes52} or LSQR~\cite{Paige82} is used. Remark
\r
308 that these methods are more appropriate than a single direct method in a
\r
309 parallel context. CGLS has recently been used to improve the performance of multisplitting algorithms \cite{cz15:ij}.
\r
313 \begin{algorithm}[t]
\r
315 \begin{algorithmic}[1]
\r
316 \Input $A$ (sparse matrix), $b$ (right-hand side)
\r
317 \Output $x$ (solution vector)\vspace{0.2cm}
\r
318 \State Set the initial guess $x_0$
\r
319 \For {$k=1,2,3,\ldots$ until convergence ($error<\epsilon_{tsirm}$)} \label{algo:conv}
\r
320 \State $[x_k,error]=Solve(A,b,x_{k-1},max\_iter_{kryl})$ \label{algo:solve}
\r
321 \State $S_{k \mod s}=x_k$ \label{algo:store} \Comment{update column ($k \mod s$) of $S$}
\r
322 \If {$k \mod s=0$ {\bf and} $error>\epsilon_{kryl}$}
\r
323 \State $R=AS$ \Comment{compute dense matrix} \label{algo:matrix_mul}
\r
324 \State $\alpha=Least\_Squares(R,b,max\_iter_{ls})$ \label{algo:}
\r
325 \State $x_k=S\alpha$ \Comment{compute new solution}
\r
332 Algorithm~\ref{algo:01} summarizes the principle of the proposed method. The
\r
333 outer iteration is inside the \emph{for} loop. Line~\ref{algo:solve}, the Krylov
\r
334 method is called for a maximum of $max\_iter_{kryl}$ iterations. In practice,
\r
335 we suggest to set this parameter equal to the restart number in the GMRES-like
\r
336 method. Moreover, a tolerance threshold must be specified for the solver. In
\r
337 practice, this threshold must be much smaller than the convergence threshold of
\r
338 the TSIRM algorithm (\emph{i.e.}, $\epsilon_{tsirm}$). We also consider that
\r
339 after the call of the $Solve$ function, we obtain the vector $x_k$ and the
\r
340 $error$, which is defined by $||Ax_k-b||_2$.
\r
342 Line~\ref{algo:store}, $S_{k \mod s}=x_k$ consists in copying the solution
\r
343 $x_k$ into the column $k \mod s$ of $S$. After the minimization, the matrix
\r
344 $S$ is reused with the new values of the residuals. To solve the minimization
\r
345 problem, an iterative method is used. Two parameters are required for that:
\r
346 the maximum number of iterations ($max\_iter_{ls}$) and the threshold to stop
\r
347 the method ($\epsilon_{ls}$).
\r
349 Let us summarize the most important parameters of TSIRM:
\r
351 \item $\epsilon_{tsirm}$: the threshold that stops the TSIRM method;
\r
352 \item $max\_iter_{kryl}$: the maximum number of iterations for the Krylov method;
\r
353 \item $s$: the number of outer iterations before applying the minimization step;
\r
354 \item $max\_iter_{ls}$: the maximum number of iterations for the iterative least-squares method;
\r
355 \item $\epsilon_{ls}$: the threshold used to stop the least-squares method.
\r
359 The parallelization of TSIRM relies on the parallelization of all its
\r
360 parts. More precisely, except the least-squares step, all the other parts are
\r
361 obvious to achieve out in parallel. In order to develop a parallel version of
\r
362 our code, we have chosen to use PETSc~\cite{petsc-web-page}. In
\r
363 line~\ref{algo:matrix_mul}, the matrix-matrix multiplication is implemented and
\r
364 efficient since the matrix $A$ is sparse and the matrix $S$ contains few columns
\r
365 in practice. As explained previously, at least two methods seem to be
\r
366 interesting to solve the least-squares minimization, the CGLS and the LSQR
\r
369 %% In Algorithm~\ref{algo:02} we remind the CGLS algorithm. The LSQR method follows
\r
370 %% more or less the same principle but it takes more place, so we briefly explain
\r
371 %% the parallelization of CGLS which is similar to LSQR.
\r
373 %% \begin{algorithm}[t]
\r
375 %% \begin{algorithmic}[1]
\r
376 %% \Input $A$ (matrix), $b$ (right-hand side)
\r
377 %% \Output $x$ (solution vector)\vspace{0.2cm}
\r
378 %% \State Let $x_0$ be an initial approximation
\r
379 %% \State $r_0=b-Ax_0$
\r
380 %% \State $p_1=A^Tr_0$
\r
381 %% \State $s_0=p_1$
\r
382 %% \State $\gamma=||s_0||^2_2$
\r
383 %% \For {$k=1,2,3,\ldots$ until convergence ($\gamma<\epsilon_{ls}$)} \label{algo2:conv}
\r
384 %% \State $q_k=Ap_k$
\r
385 %% \State $\alpha_k=\gamma/||q_k||^2_2$
\r
386 %% \State $x_k=x_{k-1}+\alpha_kp_k$
\r
387 %% \State $r_k=r_{k-1}-\alpha_kq_k$
\r
388 %% \State $s_k=A^Tr_k$
\r
389 %% \State $\gamma_{old}=\gamma$
\r
390 %% \State $\gamma=||s_k||^2_2$
\r
391 %% \State $\beta_k=\gamma/\gamma_{old}$
\r
392 %% \State $p_{k+1}=s_k+\beta_kp_k$
\r
394 %% \end{algorithmic}
\r
400 The PETSc code we have develop is avaiable here: {\bf a mettre} and it will soon
\r
401 be integrated with the PETSc sources. TSIRM has been implemented as any solver
\r
402 for linear systems in PETSc. As it requires to use another solver, we have used
\r
403 a very interesting feature of PETSc that enables to use a preconditioner as a
\r
404 linear system with the function {\it PCKSPGetKSP}. As the LSQR function is
\r
405 already implemented in PETSc, we have used it. CGLS was not implemented yet, so
\r
406 we have implemented it and we plan to define it as a minimization solver in
\r
407 PETSc similarly to LSQR. Both CGLS and LSQR are not complex from the computation
\r
408 point of view. They involves matrix-vector multiplications and some classical
\r
409 operations: dot product, norm, multiplication, and addition on vectors. As
\r
410 presented in Section~\ref{sec:05} the minimization step is scalable.
\r
413 %%%*********************************************************
\r
414 %%%*********************************************************
\r
416 \section{Convergence results}
\r
420 We suppose in this section that GMRES($m$) is used as solver in the TSIRM algorithm applied on a complex matrix $A$.
\r
421 Let us denote $A^\ast$ the conjugate transpose of $A$, and let $\mathfrak{R}(A)=\dfrac{1}{2} \left( A + A^\ast\right)$, $\mathfrak{I}(A)=\dfrac{1}{2i} \left( A - A^\ast\right)$.
\r
423 \subsection{$\mathfrak{R}(A)$ is positive}
\r
425 \begin{proposition}
\r
426 \label{positiveConvergent}
\r
427 If $\mathfrak{R}(A)$ is positive, then the TSIRM algorithm is convergent.
\r
432 If $\mathfrak{R}(A)$ is positive, then even if $A$ is complex, it is possible to state that
\r
433 the GMRES algorithm is convergent, see, \emph{e.g.},~\cite{Huang89}. In particular, its residual norm
\r
436 At each iterate of the TSIRM algorithm, either a GMRES iteration is realized or a least square
\r
437 resolution (to find the minimum of $||b-Ax||_2$ is achieved on the linear span of the iterated approximation vectors
\r
438 $span\left(x_{k-s+1}, x_{k-s}+2, \hdots, x_{k} \right)$
\r
439 of the last GMRES stage,
\r
441 $\operatorname{span}(S) = \left \{ {\sum_{i=1}^k \lambda_i v_i \Big| k \in \mathbb{N}, v_i \in S, \lambda _i \in \mathbb{R}} \right \}$.
\r
443 Obviously, the minimum of $||b-Ax||_2$ on the set $span\left(x_{k-s+1}, x_{k-s}+2, \hdots, x_{k} \right)$
\r
444 is lower than or equal to $||b-Ax_k||_2$, which is the last obtained GMRES-residual norm. So we can
\r
445 conclude that the intermediate stage of least square resolution inserted into the GMRES algorithm
\r
446 does not break the decreasing to zero of the GMRES-residual norm.
\r
448 In other words, the TSIRM algorithm is convergent.
\r
452 Regarding the convergence speed, we can claim that,
\r
453 \begin{proposition}
\r
455 If $A$ is a positive matrix, then the convergence of the
\r
456 TSIRM algorithm is linear.
\r
458 Furthermore, let $r_k$ be the $k$-th residue of TSIRM, then
\r
459 we have the following boundaries:
\r
461 ||r_k|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{km}{2}} ||r_0|| ,
\r
463 where $M$ is the symmetric part of $A$, $\alpha = \lambda_{min}(M)^2$ and $\beta = \lambda_{max}(A^T A)$.
\r
467 Let us first recall that, when $A$ is a positive real matrix with symmetric part $M$, then the residual norm provided at the $m$-th step of GMRES satisfies:
\r
469 ||r_m|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{m}{2}} ||r_0|| ,
\r
471 where $\alpha$ and $\beta$ are defined as in Proposition~\ref{prop:saad}.
\r
472 These well-known results can be found, \emph{e.g.}, in~\cite{Saad86}.
\r
474 We will now prove by a mathematical induction that, for each $k \in \mathbb{N}^\ast$,
\r
475 $||r_k|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{mk}{2}} ||r_0||$ when $A$ is positive.
\r
477 The base case is obvious, as for $k=1$, the TSIRM algorithm simply consists in applying GMRES($m$) once, leading to a new residual $r_1$ that follows the inductive hypothesis due to the results recalled above.
\r
479 Suppose now that the claim holds for all $m=1, 2, \hdots, k-1$, that is, $\forall m \in \{1,2,\hdots, k-1\}$, $||r_m|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{km}{2}} ||r_0||$.
\r
480 We will show that the statement holds too for $r_k$. Two situations can occur:
\r
482 \item If $k \not\equiv 0 ~(\textrm{mod}\ m)$, then the TSIRM algorithm consists in executing GMRES once. In that case and by using the inductive hypothesis, we obtain either $||r_k|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{m}{2}} ||r_{k-1}||\leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{km}{2}} ||r_0||$.
\r
483 \item Else, the TSIRM algorithm consists in two stages: a first GMRES($m$) execution leads to a temporary $x_k$ whose residue satisfies:
\r
484 $$||r_k|| \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{m}{2}} ||r_{k-1}||\leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{km}{2}} ||r_0||$$
\r
485 and a least squares resolution.
\r
486 Let $\operatorname{span}(S) = \left \{ {\sum_{i=1}^k \lambda_i v_i \Big| k \in \mathbb{N}, v_i \in S, \lambda _i \in \mathbb{R}} \right \}$ be the linear span of a set of real vectors $S$. So,\\
\r
487 $\min_{\alpha \in \mathbb{R}^s} ||b-R\alpha ||_2 = \min_{\alpha \in \mathbb{R}^s} ||b-AS\alpha ||_2$
\r
490 & = \min_{x \in span\left(S_{k-s+1}, S_{k-s+2}, \hdots, S_{k} \right)} ||b-AS\alpha ||_2\\
\r
491 & = \min_{x \in span\left(x_{k-s+1}, x_{k-s}+2, \hdots, x_{k} \right)} ||b-AS\alpha ||_2\\
\r
492 & \leqslant \min_{x \in span\left( x_{k} \right)} ||b-Ax ||_2\\
\r
493 & \leqslant \min_{\lambda \in \mathbb{R}} ||b-\lambda Ax_{k} ||_2\\
\r
494 & \leqslant ||b-Ax_{k}||_2\\
\r
496 & \leqslant \left(1-\dfrac{\alpha}{\beta}\right)^{\frac{km}{2}} ||r_0||, \\
\r
499 which concludes the induction and the proof.
\r
504 \subsection{$\mathfrak{R}(A)$ is positive definite}
\r
506 \begin{proposition}
\r
508 Convergence of the TSIRM algorithm is at least linear when $\mathfrak{R}(A)$ is
\r
509 positive definite. Furthermore, the rate of convergence is lower
\r
510 than $$\min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{m}{2}};
\r
511 \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{m}{2}}\right) ,$$
\r
512 where ${\lambda_{min}^{X}}$ (resp. ${\lambda_{max}^{X}}$) is the lowest (resp. largest) eigenvalue of matrix $X$.
\r
517 If $\mathfrak{R}(A)$ is positive definite, then it is positive, and so the TSIRM algorithm
\r
518 is convergent due to Proposition~\ref{positiveConvergent}.
\r
520 Furthermore, as stated in the proof of Proposition~\ref{positiveConvergent}, the GMRES residue is under control
\r
521 when $\mathfrak{R}(A)$ is positive. More precisely, it has been proven in the literature that the residual norm
\r
522 provided at the $m$-th step of GMRES satisfies:
\r
524 \item $||r_m|| \leqslant \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{mk}{2}} ||r_0||$, see, \emph{e.g.},~\cite{citeulike:2951999},
\r
525 \item $||r_m|| \leqslant \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mk}{2}} ||r_0||$, see~\cite{ANU:137201},
\r
527 which proves the convergence of GMRES($m$) for all $m$ under such assumptions regarding $A$.
\r
529 We will now prove by a mathematical induction, and following the same canvas than in the proof of Prop.~\ref{positiveConvergent}, that: for each $k \in \mathbb{N}^\ast$, the TSIRM-residual norm satisfies
\r
533 ||r_k|| \leqslant & \min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{m}{2}}; \right. \\
\r
534 & \left. \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{m}{2}}\right) ||r_0||
\r
537 when $A$ is positive definite.
\r
540 The base case is obvious, as for $k=1$, the TSIRM algorithm simply consists in applying GMRES($m$) once, leading to a new residual $r_1$ that follows the inductive hypothesis due to the results recalled in the items listed above.
\r
542 Suppose now that the claim holds for all $u=1, 2, \hdots, k-1$, that is, $\forall u \in \{1,2,\hdots, k-1\}$, $||r_u|| \leqslant \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mu}{2}} ||r_0||$.
\r
543 We will show that the statement holds too for $r_k$. Two situations can occur:
\r
545 \item If $k \not\equiv 0 ~(\textrm{mod}\ m)$, then the TSIRM algorithm consists in executing GMRES once. In that case and by using the inductive hypothesis, we obtain
\r
546 $||r_k|| \leqslant \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{m}{2}} \leqslant \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{mk}{2}} ||r_0||$, due to~\cite{citeulike:2951999}. Furthermore, we have too that: $||r_k|| \leqslant \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{m}{2}} ||r_{k-1}|| \leqslant \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mk}{2}} ||r_0||$, as proven in~\cite{ANU:137201} and by using the inductive hypothesis. So we can conclude that
\r
547 $$\begin{array}{ll}||r_k|| \leqslant & \min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{mk}{2}}; \right. \\
\r
548 & \left. \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mk}{2}}\right) \times ||r_0||
\r
551 \item Else, the TSIRM algorithm consists in two stages: a first GMRES($m$) execution leads to a temporary $x_k$ whose residue satisfies, following the previous item:
\r
552 $$\begin{array}{ll}
\r
553 ||r_k|| & \leqslant \min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{m}{2}}; \right. \\
\r
554 & \left. \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{m}{2}}\right) \times ||r_{k-1}||\\
\r
555 & \leqslant \min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{mk}{2}}; \right. \\
\r
556 & \left. \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mk}{2}}\right) \times ||r_0||
\r
558 and the least squares resolution of $\min_{\alpha \in \mathbb{R}^s} ||b-R\alpha ||_2$.
\r
560 Let $\operatorname{span}(S) = \left \{ {\sum_{i=1}^k \lambda_i v_i \Big| k \in \mathbb{N}, v_i \in S, \lambda _i \in \mathbb{R}} \right \}$ be the linear span of a set of real vectors $S$, as defined previously. So,\\
\r
561 $\min_{\alpha \in \mathbb{R}^s} ||b-R\alpha ||_2 = \min_{\alpha \in \mathbb{R}^s} ||b-AS\alpha ||_2$
\r
564 & = \min_{x \in span\left(S_{k-s+1}, S_{k-s+2}, \hdots, S_{k} \right)} ||b-AS\alpha ||_2\\
\r
565 & = \min_{x \in span\left(x_{k-s+1}, x_{k-s}+2, \hdots, x_{k} \right)} ||b-AS\alpha ||_2\\
\r
566 & \leqslant \min_{x \in span\left( x_{k} \right)} ||b-Ax ||_2\\
\r
567 & \leqslant \min_{\lambda \in \mathbb{R}} ||b-\lambda Ax_{k} ||_2\\
\r
568 & \leqslant ||b-Ax_{k}||_2\\
\r
570 & \leqslant \min\left( \left(1- \dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{ \lambda_{min}^{\mathfrak{R}(A)} \lambda_{max}^{\mathfrak{R}(A)} + {\lambda_{max}^{\mathfrak{I}(A)}}^2}\right)^{\frac{mk}{2}}; \right. \\
\r
571 & \left. \left(1-\dfrac{{\lambda_{min}^{\mathfrak{R}(A)}}^2}{||A||^2}\right)^{\frac{mk}{2}}\right) \times ||r_0||
\r
574 due to the inductive hypothesis.
\r
575 So the statement of Equation~\eqref{induc} holds too for the $k$-th iterate, which concludes the induction and the proof.
\r
578 \subsection{A last linear convergence}
\r
581 \begin{proposition}
\r
582 Let us define the field of values of $A$ by
\r
583 $$\mathfrak{F}(A) = \left\{ \dfrac{x^\ast A x}{x^\ast x}, x \in \mathds{C}^n\setminus \{0\} \right\} .$$
\r
585 Then if $\mathfrak{F}(A)$ is included into a closed ball of radius $r$ and center $c$,
\r
586 which does not contain the origin, then the convergence of the TSIRM algorithm is at least linear.
\r
588 More precisely, the rate of convergence is lower
\r
589 than $2 \dfrac{r}{|c|}$.
\r
593 This inequality comes from the fact that, in the conditions of the proposition, the GMRES residue
\r
594 satisfies the inequality: $|r_k| \leqslant 2 \dfrac{r}{|c|}^k |r_0|$. An induction inspired by
\r
595 the proofs of Propositions~\ref{prop:saad} and~\ref{prop2} can transfer this inequality to the
\r
601 Remark that a similar proposition can be formulated at each time
\r
602 the given solver satisfies an inequality of the form $||r_n|| \leqslant \mu^n ||r_0||$,
\r
603 with $|\mu|<1$. Furthermore, it is \emph{a priori} possible in some particular cases
\r
605 that the proposed TSIRM converges while the GMRES($m$) does not.
\r
608 %%%*********************************************************
\r
609 %%%*********************************************************
\r
610 \section{Experiments using PETSc}
\r
613 In this section four kinds of experiments have been performed. First, some
\r
614 experiments on real matrices issued from the sparse matrix florida have been
\r
615 achieved out. Second, some experiments in parallel with some linear problems are
\r
616 reported and analyzed. Third, some experiments in parallèle with som nonlinear
\r
617 problems are illustrated. Finally some parameters of TSIRM are studied in order
\r
618 to understand their influences.
\r
621 \subsection{Real matrices}
\r
625 In order to see the behavior of our approach when considering only one processor,
\r
626 a first comparison with GMRES or FGMRES and the new algorithm detailed
\r
627 previously has been experimented. Matrices that have been used with their
\r
628 characteristics (names, fields, rows, and nonzero coefficients) are detailed in
\r
629 Table~\ref{tab:01}. These latter, which are real-world applications matrices,
\r
630 have been extracted from the Davis collection, University of
\r
631 Florida~\cite{Dav97}.
\r
633 \begin{table*}[htbp]
\r
635 \begin{tabular}{|c|c|r|r|r|}
\r
637 Matrix name & Field &\# Rows & \# Nonzeros \\\hline \hline
\r
638 crashbasis & Optimization & 160,000 & 1,750,416 \\
\r
639 parabolic\_fem & Comput. fluid dynamics & 525,825 & 2,100,225 \\
\r
640 epb3 & Thermal problem & 84,617 & 463,625 \\
\r
641 atmosmodj & Comput. fluid dynamics & 1,270,432 & 8,814,880 \\
\r
642 bfwa398 & Electromagnetics pb & 398 & 3,678 \\
\r
643 torso3 & 2D/3D problem & 259,156 & 4,429,042 \\
\r
647 \caption{Main characteristics of the sparse matrices chosen from the Davis collection}
\r
651 Chosen parameters are detailed below.
\r
652 We have stopped the GMRES every 30
\r
653 iterations (\emph{i.e.}, $max\_iter_{kryl}=30$), which is the default
\r
654 setting of GMRES restart parameter. The parameter $s$ has been set to 8. CGLS
\r
655 minimizes the least-squares problem with parameters
\r
656 $\epsilon_{ls}=1e-40$ and $max\_iter_{ls}=20$. The external precision is set to
\r
657 $\epsilon_{tsirm}=1e-10$. These experiments have been performed on an Intel(R)
\r
658 Core(TM) i7-3630QM CPU @ 2.40GHz with the 3.5.1 version of PETSc.
\r
661 Experiments comparing
\r
662 a GMRES variant with TSIRM in the resolution of linear systems are given in Table~\ref{tab:02}.
\r
663 The second column describes whether GMRES or FGMRES has been used for linear systems solving.
\r
664 Different preconditioners have been used according to the matrices. With TSIRM, the same
\r
665 solver and the same preconditioner are used. This table shows that TSIRM can
\r
666 drastically reduce the number of iterations needed to reach the convergence, when the
\r
667 number of iterations for the normal GMRES is more or less greater than 500. In
\r
668 fact this also depends on two parameters: the number of iterations before stopping GMRES
\r
669 and the number of iterations to perform the minimization.
\r
672 \begin{table*}[htbp]
\r
674 \begin{tabular}{|c|c|r|r|r|r|}
\r
677 \multirow{2}{*}{Matrix name} & Solver / & \multicolumn{2}{c|}{GMRES} & \multicolumn{2}{c|}{TSIRM CGLS} \\
\r
679 & precond & Time & \# Iter. & Time & \# Iter. \\\hline \hline
\r
681 crashbasis & gmres / none & 15.65 & 518 & 14.12 & 450 \\
\r
682 parabolic\_fem & gmres / ilu & 1009.94 & 7573 & 401.52 & 2970 \\
\r
683 epb3 & fgmres / sor & 8.67 & 600 & 8.21 & 540 \\
\r
684 atmosmodj & fgmres / sor & 104.23 & 451 & 88.97 & 366 \\
\r
685 bfwa398 & gmres / none & 1.42 & 9612 & 0.28 & 1650 \\
\r
686 torso3 & fgmres / sor & 37.70 & 565 & 34.97 & 510 \\
\r
690 \caption{Comparison between sequential standalone (F)GMRES and TSIRM with (F)GMRES (time in seconds).}
\r
697 \subsection{Parallel linear problems}
\r
700 In order to perform larger experiments, we have tested some example applications
\r
701 of PETSc. These applications are available in the \emph{ksp} part, which is
\r
702 suited for scalable linear equations solvers:
\r
704 \item ex15 is an example that solves in parallel an operator using a finite
\r
705 difference scheme. The diagonal is equal to 4 and 4 extra-diagonals
\r
706 representing the neighbors in each directions are equal to -1. This example is
\r
707 used in many physical phenomena, for example, heat and fluid flow, wave
\r
709 \item ex54 is another example based on a 2D problem discretized with quadrilateral
\r
710 finite elements. In this example, the user can define the scaling of material
\r
711 coefficient in embedded circle called $\alpha$.
\r
713 For more technical details on these applications, interested readers are invited
\r
714 to read the codes available in the PETSc sources. These problems have been
\r
715 chosen because they are scalable with many cores.
\r
717 In the following, larger experiments are described on two large scale
\r
718 architectures: Curie and Juqueen. Both these architectures are supercomputers
\r
719 respectively composed of 80,640 cores for Curie and 458,752 cores for
\r
720 Juqueen. Those machines are respectively hosted by GENCI in France and Jülich
\r
721 Supercomputing Center in Germany. They belong with other similar architectures
\r
722 to the PRACE initiative (Partnership for Advanced Computing in Europe), which
\r
723 aims at proposing high performance supercomputing architecture to enhance
\r
724 research in Europe. The Curie architecture is composed of Intel E5-2680
\r
725 processors at 2.7 GHz with 2Gb memory by core. The Juqueen architecture,
\r
727 composed by IBM PowerPC A2 at 1.6 GHz with 1Gb memory per core. Both those
\r
728 architectures are equipped with a dedicated high speed network.
\r
731 In many situations, using preconditioners is essential in order to find the
\r
732 solution of a linear system. There are many preconditioners available in PETSc.
\r
733 However, for parallel applications, all the preconditioners based on matrix factorization
\r
734 are not available. In our experiments, we have tested different kinds of
\r
735 preconditioners, but as it is not the subject of this paper, we will not
\r
736 present results with many preconditioners. In practice, we have chosen to use a
\r
737 multigrid (MG) and successive over-relaxation (SOR). For further details on the
\r
738 preconditioners in PETSc, readers are referred to~\cite{petsc-web-page}.
\r
742 \begin{table*}[htbp]
\r
744 \begin{tabular}{|r|r|r|r|r|r|r|r|r|}
\r
747 nb. cores & precond & \multicolumn{2}{c|}{FGMRES} & \multicolumn{2}{c|}{TSIRM CGLS} & \multicolumn{2}{c|}{TSIRM LSQR} & best gain \\
\r
749 & & Time & \# Iter. & Time & \# Iter. & Time & \# Iter. & \\\hline \hline
\r
750 2,048 & MG & 403.49 & 18,210 & 73.89 & 3,060 & 77.84 & 3,270 & 5.46 \\
\r
751 2,048 & SOR & 745.37 & 57,060 & 87.31 & 6,150 & 104.21 & 7,230 & 8.53 \\
\r
752 4,096 & MG & 562.25 & 25,170 & 97.23 & 3,990 & 89.71 & 3,630 & 6.27 \\
\r
753 4,096 & SOR & 912.12 & 70,194 & 145.57 & 9,750 & 168.97 & 10,980 & 6.26 \\
\r
754 8,192 & MG & 917.02 & 40,290 & 148.81 & 5,730 & 143.03 & 5,280 & 6.41 \\
\r
755 8,192 & SOR & 1,404.53 & 106,530 & 212.55 & 12,990 & 180.97 & 10,470 & 7.76 \\
\r
756 16,384 & MG & 1,430.56 & 63,930 & 237.17 & 8,310 & 244.26 & 7,950 & 6.03 \\
\r
757 16,384 & SOR & 2,852.14 & 216,240 & 418.46 & 21,690 & 505.26 & 23,970 & 6.82 \\
\r
761 \caption{Comparison of FGMRES and TSIRM with FGMRES for example ex15 of PETSc/KSP with two preconditioners (MG and SOR) having 25,000 components per core on Juqueen ($\epsilon_{tsirm}=1e-3$, $max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$, $\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
766 Table~\ref{tab:03} shows the execution times and the number of iterations of
\r
767 example ex15 of PETSc on the Juqueen architecture. Different numbers of cores
\r
768 are studied ranging from 2,048 up-to 16,383 with the two preconditioners {\it
\r
769 MG} and {\it SOR}. For those experiments, the number of components (or
\r
770 unknowns of the problems) per core is fixed at 25,000, also called weak
\r
771 scaling. This number can seem relatively small. In fact, for some applications
\r
772 that need a lot of memory, the number of components per processor requires
\r
773 sometimes to be small. Other parameters for this application are described in
\r
774 the legend of this table.
\r
778 In Table~\ref{tab:03}, we can notice that TSIRM is always faster than
\r
779 FGMRES. The last column shows the ratio between FGMRES and the best version of
\r
780 TSIRM according to the minimization procedure: CGLS or LSQR. Even if we have
\r
781 computed the worst case between CGLS and LSQR, it is clear that TSIRM is always
\r
782 faster than FGMRES. For this example, the multigrid preconditioner is faster
\r
783 than SOR. The gain between TSIRM and FGMRES is more or less similar for the two
\r
784 preconditioners. Looking at the number of iterations to reach the convergence,
\r
785 it is obvious that TSIRM allows the reduction of the number of iterations. It
\r
786 should be noticed that for TSIRM, in those experiments, only the iterations of
\r
787 the Krylov solver are taken into account. Iterations of CGLS or LSQR were not
\r
788 recorded but they are time-consuming. In general, each $max\_iter_{kryl}*s$
\r
789 iterations which corresponds to 30*12, there are $max\_iter_{ls}$ iterations for
\r
790 the least-squares method which corresponds to 15.
\r
792 \begin{figure}[htbp]
\r
794 \includegraphics[width=0.5\textwidth]{nb_iter_sec_ex15_juqueen}
\r
795 \caption{Number of iterations per second with ex15 and the same parameters as in Table~\ref{tab:03} (weak scaling)}
\r
800 In Figure~\ref{fig:01}, the number of iterations per second corresponding to
\r
801 Table~\ref{tab:03} is displayed. It can be noticed that the number of
\r
802 iterations per second of FMGRES is constant whereas it decreases with TSIRM with
\r
803 both preconditioners. This can be explained by the fact that when the number of
\r
804 cores increases, the time for the least-squares minimization step also increases
\r
805 but, generally, when the number of cores increases, the number of iterations to
\r
806 reach the threshold also increases, and, in that case, TSIRM is more efficient
\r
807 to reduce the number of iterations. So, the overall benefit of using TSIRM is
\r
815 \begin{table*}[htbp]
\r
817 \begin{tabular}{|r|r|r|r|r|r|r|r|r|}
\r
820 nb. cores & $\epsilon_{tsirm}$ & \multicolumn{2}{c|}{FGMRES} & \multicolumn{2}{c|}{TSIRM CGLS} & \multicolumn{2}{c|}{TSIRM LSQR} & best gain \\
\r
822 & & Time & \# Iter. & Time & \# Iter. & Time & \# Iter. & \\\hline \hline
\r
823 2,048 & 8e-5 & 108.88 & 16,560 & 23.06 & 3,630 & 22.79 & 3,630 & 4.77 \\
\r
824 2,048 & 6e-5 & 194.01 & 30,270 & 35.50 & 5,430 & 27.74 & 4,350 & 6.99 \\
\r
825 4,096 & 7e-5 & 160.59 & 22,530 & 35.15 & 5,130 & 29.21 & 4,350 & 5.49 \\
\r
826 4,096 & 6e-5 & 249.27 & 35,520 & 52.13 & 7,950 & 39.24 & 5,790 & 6.35 \\
\r
827 8,192 & 6e-5 & 149.54 & 17,280 & 28.68 & 3,810 & 29.05 & 3,990 & 5.21 \\
\r
828 8,192 & 5e-5 & 785.04 & 109,590 & 76.07 & 10,470 & 69.42 & 9,030 & 11.30 \\
\r
829 16,384 & 4e-5 & 718.61 & 86,400 & 98.98 & 10,830 & 131.86 & 14,790 & 7.26 \\
\r
833 \caption{Comparison of FGMRES and TSIRM with FGMRES algorithms for ex54 of PETSc/KSP (both with the MG preconditioner) with 25,000 components per core on Curie ($max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$, $\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
839 In Table~\ref{tab:04}, some experiments with example ex54 on the Curie
\r
840 architecture are reported. For this application, we fixed $\alpha=0.6$. As it
\r
841 can be seen in that table, the size of the problem has a strong influence on the
\r
842 number of iterations to reach the convergence. That is why we have preferred to
\r
843 change the threshold. If we set it to $1e-3$ as with the previous application,
\r
844 only one iteration is necessary to reach the convergence. So Table~\ref{tab:04}
\r
845 shows the results of different executions with different number of cores and
\r
846 different thresholds. As with the previous example, we can observe that TSIRM is
\r
847 faster than FGMRES. The ratio greatly depends on the number of iterations for
\r
848 FMGRES to reach the threshold. The greater the number of iterations to reach the
\r
849 convergence is, the better the ratio between our algorithm and FMGRES is. This
\r
850 experiment is also a weak scaling with approximately $25,000$ components per
\r
851 core. It can also be observed that the difference between CGLS and LSQR is not
\r
852 significant. Both can be good but it seems not possible to know in advance which
\r
853 one will be the best.
\r
855 Table~\ref{tab:05} shows a strong scaling experiment with example ex54 on the
\r
856 Curie architecture. So, in this case, the number of unknowns is fixed at
\r
857 $204,919,225$ and the number of cores ranges from $512$ to $8192$ with the power
\r
858 of two. The threshold is fixed at $5e-5$ and only the $mg$ preconditioner has
\r
859 been tested. Here again we can see that TSIRM is faster than FGMRES. The
\r
860 efficiency of each algorithm is reported. It can be noticed that the efficiency
\r
861 of FGMRES is better than the TSIRM one except with $8,192$ cores and that its
\r
862 efficiency is greater than one whereas the efficiency of TSIRM is lower than
\r
863 one. Nevertheless, the ratio of TSIRM with any version of the least-squares
\r
864 method is always faster. With $8,192$ cores when the number of iterations is
\r
865 far more important for FGMRES, we can see that it is only slightly more
\r
866 important for TSIRM.
\r
868 In Figure~\ref{fig:02} we report the number of iterations per second for the
\r
869 experiments reported in Table~\ref{tab:05}. This figure highlights that the
\r
870 number of iterations per second is more or less the same for FGMRES and TSIRM
\r
871 with a little advantage for FGMRES. It can be explained by the fact that, as we
\r
872 have previously explained, the iterations of the least-squares steps are not
\r
873 taken into account with TSIRM.
\r
875 \begin{table*}[htbp]
\r
877 \begin{tabular}{|r|r|r|r|r|r|r|r|r|r|r|}
\r
880 nb. cores & \multicolumn{2}{c|}{FGMRES} & \multicolumn{2}{c|}{TSIRM CGLS} & \multicolumn{2}{c|}{TSIRM LSQR} & best gain & \multicolumn{3}{c|}{efficiency} \\
\r
881 \cline{2-7} \cline{9-11}
\r
882 & Time & \# Iter. & Time & \# Iter. & Time & \# Iter. & & FGMRES & TS CGLS & TS LSQR\\\hline \hline
\r
883 512 & 3,969.69 & 33,120 & 709.57 & 5,790 & 622.76 & 5,070 & 6.37 & 1 & 1 & 1 \\
\r
884 1024 & 1,530.06 & 25,860 & 290.95 & 4,830 & 307.71 & 5,070 & 5.25 & 1.30 & 1.21 & 1.01 \\
\r
885 2048 & 919.62 & 31,470 & 237.52 & 8,040 & 194.22 & 6,510 & 4.73 & 1.08 & .75 & .80\\
\r
886 4096 & 405.60 & 28,380 & 111.67 & 7,590 & 91.72 & 6,510 & 4.42 & 1.22 & .79 & .84 \\
\r
887 8192 & 785.04 & 109,590 & 76.07 & 10,470 & 69.42 & 9,030 & 11.30 & .32 & .58 & .56 \\
\r
892 \caption{Comparison of FGMRES and TSIRM for ex54 of PETSc/KSP (both with the MG preconditioner) with 204,919,225 components on Curie with different number of cores ($\epsilon_{tsirm}=5e-5$, $max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$, $\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
897 \begin{figure}[htbp]
\r
899 \includegraphics[width=0.5\textwidth]{nb_iter_sec_ex54_curie}
\r
900 \caption{Number of iterations per second with ex54 and the same parameters as in Table~\ref{tab:05} (strong scaling)}
\r
905 \begin{figure}[htbp]
\r
907 \includegraphics[width=0.5\textwidth]{nb_iter_sec_ex45_curie}
\r
908 \caption{Number of iterations per second with ex45 and the same parameters as in Table~\ref{tab:06} (weak scaling)}
\r
914 It is well-known that preconditioners have a very strong influence on the
\r
915 convergence of linear systems. Previously, we have used some classical
\r
916 preconditioners provided by PETSc. HYPRE~\cite{Falgout06} is a very efficient
\r
917 preconditioner based on structured multigrid and element-based algebraic
\r
918 multigrid algorithms. In Table~\ref{tab:06} we report an experiment that show it
\r
919 reduces drastivally the number of iterations but sometimes it is very
\r
920 time-consuming compared to other simpler precondititioners. In this table, we
\r
921 can see that for $512$ and $2,048$ cores, HYPRE reduces drastically the number
\r
922 of iterations for FGMRES to reach the convergence. However, it is very
\r
923 time-consuming compared to TSIRM and FGMRES with the ASM preconditioner. For
\r
924 $4,096$ and $8,192$ cores, FGMRES with HYPRE did not converge in less than 1000s
\r
925 where FGMRES and TSIRM with the ASM converge very quickly. Finally, it can be
\r
926 noticed that TSIRM is also faster than FGMRES and it requires less iterations.
\r
929 \begin{table*}[htbp]
\r
931 \begin{tabular}{|r|r|r|r|r|r|r|r|}
\r
934 nb. cores & \multicolumn{2}{c|}{FGMRES/ASM} & \multicolumn{2}{c|}{TSIRM CGLS/ASM} & gain& \multicolumn{2}{c|}{FGMRES/HYPRE} \\
\r
935 \cline{2-5} \cline{7-8}
\r
936 & Time & \# Iter. & Time & \# Iter. & & Time & \# Iter. \\\hline \hline
\r
937 512 & 5.54 & 685 & 2.5 & 570 & 2.21 & 128.9 & 9 \\
\r
938 2048 & 14.95 & 1,560 & 4.32 & 746 & 3.48 & 335.7 & 9 \\
\r
939 4096 & 25.13 & 2,369 & 5.61 & 859 & 4.48 & >1000 & -- \\
\r
940 8192 & 44.35 & 3,197 & 7.6 & 1,083 & 5.84 & >1000 & -- \\
\r
945 \caption{Comparison of FGMRES and TSIRM for ex45 of PETSc/KSP with two preconditioner (ASM and HYPRE) having 5,000 components per core on Curie ($\epsilon_{tsirm}=1e-10$, $max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$,$\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
951 \subsection{Parallel nonlinear problems}
\r
953 With PETSc, linear solvers are used inside nonlinear solvers. The SNES
\r
954 (Scalable Nonlinear Equations Solvers) module in PETSc implements easy to use
\r
955 methods, like Newton-type, quasi-Newton or full approximation scheme (FAS)
\r
956 multigrid to solve systems of nonlinears equations. As SNES is based on the
\r
957 Krylov methods of PETSc, it is interesting to investigate if the TSIRM method is
\r
958 also efficient and scalable with non linear problems. In PETSc, some examples
\r
959 are provided. An important criteria is the scalability of the initial code with
\r
960 classical solvers. Consequently, we have chosen two of these examples: ex14 and
\r
961 ex20. In ex14, the code solves the Bratu (SFI - solid fuel ignition) nonlinear
\r
962 partial difference equations in 3 dimension. In ex20, the code solves a 3
\r
963 dimension radiative transport test problem. For more details on these examples,
\r
964 interested readers are invited to see the code in the PETSc examples. For both
\r
965 these examples, a weak scaling case is chosen where processors have
\r
966 approximately a number of components equals to 100,000.
\r
968 In Table~\ref{tab:07} we report the result of our experiments for the example
\r
969 ex14 with the block Jacobi preconditioner. For TSIRM the CGLS algorithm is used
\r
970 to solve the minimization step. In this table, we can see that the number of
\r
971 iterations used by the linear solver is smaller with TSIRM compared with FGMRES.
\r
972 Consequently the execution times are smaller with TSIRM. The gain between TSIRM
\r
973 and FGMRES is around 6 and 7. The parameters of TSIRM are expressed in the
\r
974 caption of the table.
\r
976 \begin{table*}[htbp]
\r
978 \begin{tabular}{|r|r|r|r|r|r|}
\r
981 nb. cores & \multicolumn{2}{c|}{FGMRES/BJAC} & \multicolumn{2}{c|}{TSIRM CGLS/BJAC} & gain \\
\r
983 & Time & \# Iter. & Time & \# Iter. & \\\hline \hline
\r
984 1,024 & 159.52 & 11,584 & 26.34 & 1,563 & 6.06 \\
\r
985 2,048 & 226.24 & 16,459 & 37.23 & 2,248 & 6.08\\
\r
986 4,096 & 391.21 & 27,794 & 50.93 & 2,911 & 7.69\\
\r
987 8,192 & 543.23 & 37,770 & 79.21 & 4,324 & 6.86 \\
\r
992 \caption{Comparison of FGMRES and TSIRM for ex14 of PETSc/SNES with a Block Jacobi preconditioner having 100,000 components per core on Curie ($\epsilon_{tsirm}=1e-10$, $max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$, $\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
997 In Table~\ref{tab:08}, the results of the experiments with the example ex20 are
\r
998 reported. The block Jacobi preconditioner has also been used and CGLS to solve
\r
999 the minimization step for TSIRM. For this example, we can observ that the number
\r
1000 of iterations for FMGRES increase drastically when the number of cores
\r
1001 increases. With TSIRM, we can see that the number of iterations is initially
\r
1002 very small compared to the FGMRES ones and when the number of cores increase,
\r
1003 the number of iterations increases slighther with TSIRM than with FGMRES. For
\r
1004 this example, the gain between TSIRM and FGMRES ranges between 8 with 1,024
\r
1005 cores to more than 16 with 8,192 cores.
\r
1007 \begin{table*}[htbp]
\r
1009 \begin{tabular}{|r|r|r|r|r|r|}
\r
1012 nb. cores & \multicolumn{2}{c|}{FGMRES/BJAC} & \multicolumn{2}{c|}{TSIRM CGLS/BJAC} & gain \\
\r
1014 & Time & \# Iter. & Time & \# Iter. & \\\hline \hline
\r
1015 1,024 & 667.92 & 48,732 & 81.65 & 5,087 & 8.18 \\
\r
1016 2,048 & 966.87 & 77,177 & 90.34 & 5,716 & 10.70\\
\r
1017 4,096 & 1,742.31 & 124,411 & 119.21 & 6,905 & 14.61\\
\r
1018 8,192 & 2,739.21 & 187,626 & 168.9 & 9,000 & 16.22\\
\r
1023 \caption{Comparison of FGMRES and TSIRM for ex20 of PETSc/SNES with a Block Jacobi preconditioner having 100,000 components per core on Curie ($\epsilon_{tsirm}=1e-10$, $max\_iter_{kryl}=30$, $s=12$, $max\_iter_{ls}=15$, $\epsilon_{ls}=1e-40$), time is expressed in seconds.}
\r
1035 \subsection{Influence of parameters for TSIRM}
\r
1036 In this section we present some experimental results in order to study the influence of some parameters on the TSIRM algorithm. We conducted experiments on $16$ cores to solve 3D problems of size $200,000$ components per core. We solved nonlinear problems token from examples of PETSc. We fixed some parameters of the TSIRM algorithm as follows: the nonlinear systems are solved with a precision of $10^{-8}$, block Jacobi preconditioner is used, the tolerance threshold $\epsilon_{tsirm}$ is $10^{-8}$ , the maximum number of iterations $max\_iter_{tsirm}$ is set to $10,000$ iterations, the FGMRES method is used as the inner solver with a tolerance threshold $\epsilon_{kryl}=10^{-10}$ and the least-squares problem is solved with a precision $\epsilon_{ls}=10^{-40}$ in the minimization process.
\r
1038 %time mpirun ../ex34 -da_grid_x 147 -da_grid_y 147 -da_grid_z 147 -ksp_type tsirm -ksp_pc_type asm -pc_type ksp -ksp_tsirm_tol 1e-10 -ksp_tsirm_maxiter 10000 -ksp_ksp_type fgmres -ksp_tsirm_max_inner_iter 30 -ksp_tsirm_inner_tol 1e-10 -ksp_tsirm_cgls 0 -ksp_tsirm_tol_ls 1.e-40 -ksp_tsirm_maxiter_ls 20 -ksp_tsirm_size_ls 10
\r
1039 \begin{figure}[htbp]
\r
1041 \includegraphics[width=0.5\textwidth]{ksp_tsirm_cgls}
\r
1042 \caption{Number of total iterations using two different methods for the minimization: CGLS and LSQR.}
\r
1046 \begin{figure}[htbp]
\r
1048 \includegraphics[width=0.5\textwidth]{snes_ex14}
\r
1049 \caption{Total number of iterations in example {\it snes ex14} of PETSc by varyin the number of inner iterations and the size of the least-squares problem.}
\r
1050 \label{fig:snes_ex14}
\r
1062 \subsection{Experiments conclusions }
\r
1066 Concerning the experiments some other remarks are interesting.
\r
1068 \item We have tested other examples of PETSc/KSP (ex29, ex45, ex49). For all these
\r
1069 examples, we have also obtained similar gains between GMRES and TSIRM but
\r
1070 those examples are not scalable with many cores. In general, we had some
\r
1071 problems with more than $4,096$ cores.
\r
1072 \item We have tested many iterative solvers available in PETSc. In fact, it is
\r
1073 possible to use most of them with TSIRM. From our point of view, the condition
\r
1074 to use a solver inside TSIRM is that the solver must have a restart
\r
1075 feature. More precisely, the solver must support to be stopped and restarted
\r
1076 without decreasing its convergence. That is why with GMRES we stop it when it
\r
1077 is naturally restarted (\emph{i.e.} with $m$ the restart parameter). The
\r
1078 Conjugate Gradient (CG) and all its variants do not have ``restarted'' version
\r
1079 in PETSc, so they are not efficient. They will converge with TSIRM but not
\r
1080 quickly because if we compare a normal CG with a CG which is stopped and
\r
1081 restarted every 16 iterations (for example), the normal CG will be far more
\r
1082 efficient. Some restarted CG or CG variant versions exist and may be
\r
1083 interesting to study in future works.
\r
1089 %%%*********************************************************
\r
1090 %%%*********************************************************
\r
1091 \section{Conclusion}
\r
1093 %The conclusion goes here. this is more of the conclusion
\r
1094 %%%*********************************************************
\r
1095 %%%*********************************************************
\r
1098 In this paper a new two-stage algorithm TSIRM has been described. This method allows us to improve the convergence of Krylov iterative methods. It is based
\r
1099 on a least-squares minimization step which uses the Krylov residuals.
\r
1102 We have implemented our code in PETSc in order to show that it is efficient and scalable. Some experiments with classical examples of PETSc for linear and nonlinear problems have been performed. We observed that TSIRM outperforms GMRES variants when the number of iterations is large. TSIRM is also scalable since we made some experiments with up to 16,394 cores.
\r
1104 We also observed that TSIRM is efficient with different preconditioners. The hypre preconditioner that is globally very efficient for many problems is also very time consuming. Consequently, sometimes using a less performent preconditioners may be a better solution. In that case, TSIRM is also more efficient than traditional Krylov methods.
\r
1106 {\bf A CHECKER !!}
\r
1107 The influence of some important parameters of TSIRM have been studied. It can be noticed that they have a strong influence on the convergence speed
\r
1109 In future works, we plan to study other problems coming from different research areas. Other efficient Krylov optimisation methods as communication avoiding technique may be interesting to be investigated
\r
1114 % use section* for acknowledgement
\r
1115 %%%*********************************************************
\r
1116 %%%*********************************************************
\r
1117 \section*{Acknowledgment}
\r
1118 This paper is partially funded by the Labex ACTION program (contract
\r
1119 ANR-11-LABX-01-01). We acknowledge PRACE for awarding us access to resources
\r
1120 Curie and Juqueen respectively based in France and Germany.
\r
1126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\r
1128 \bibliography{biblio}
\r
1129 \bibliographystyle{plain}
\r
1130 %\bibliographystyle{alpha}
\r
