1 \documentclass[times]{cpeauth}
5 %\usepackage[dvips,colorlinks,bookmarksopen,bookmarksnumbered,citecolor=red,urlcolor=red]{hyperref}
7 %\newcommand\BibTeX{{\rmfamily B\kern-.05em \textsc{i\kern-.025em b}\kern-.08em
8 %T\kern-.1667em\lower.7ex\hbox{E}\kern-.125emX}}
16 \usepackage[T1]{fontenc}
17 \usepackage[utf8]{inputenc}
18 \usepackage{amsfonts,amssymb}
20 \usepackage{algorithm}
21 \usepackage{algpseudocode}
24 % Extension pour les liens intra-documents (tagged PDF)
25 % et l'affichage correct des URL (commande \url{http://example.com})
26 %\usepackage{hyperref}
31 \DeclareUrlCommand\email{\urlstyle{same}}
33 \usepackage[autolanguage,np]{numprint}
35 \renewcommand*\npunitcommand[1]{\text{#1}}
36 \npthousandthpartsep{}}
39 \usepackage[textsize=footnotesize]{todonotes}
41 \newcommand{\AG}[2][inline]{%
42 \todo[color=green!50,#1]{\sffamily\textbf{AG:} #2}\xspace}
43 \newcommand{\RC}[2][inline]{%
44 \todo[color=red!10,#1]{\sffamily\textbf{RC:} #2}\xspace}
45 \newcommand{\LZK}[2][inline]{%
46 \todo[color=blue!10,#1]{\sffamily\textbf{LZK:} #2}\xspace}
47 \newcommand{\RCE}[2][inline]{%
48 \todo[color=yellow!10,#1]{\sffamily\textbf{RCE:} #2}\xspace}
49 \newcommand{\DL}[2][inline]{%
50 \todo[color=pink!10,#1]{\sffamily\textbf{DL:} #2}\xspace}
52 \algnewcommand\algorithmicinput{\textbf{Input:}}
53 \algnewcommand\Input{\item[\algorithmicinput]}
55 \algnewcommand\algorithmicoutput{\textbf{Output:}}
56 \algnewcommand\Output{\item[\algorithmicoutput]}
58 \newcommand{\TOLG}{\mathit{tol_{gmres}}}
59 \newcommand{\MIG}{\mathit{maxit_{gmres}}}
60 \newcommand{\TOLM}{\mathit{tol_{multi}}}
61 \newcommand{\MIM}{\mathit{maxit_{multi}}}
62 \newcommand{\TOLC}{\mathit{tol_{cgls}}}
63 \newcommand{\MIC}{\mathit{maxit_{cgls}}}
66 \usepackage{color, colortbl}
67 \newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
68 \newcolumntype{Z}[1]{>{\raggedleft}m{#1}}
70 \newcolumntype{g}{>{\columncolor{Gray}}c}
71 \definecolor{Gray}{gray}{0.9}
76 \title{Grid-enabled simulation of large-scale linear iterative solvers}
77 %\itshape{\journalnamelc}\footnotemark[2]}
79 \author{Charles Emile Ramamonjisoa\affil{1},
80 Lilia Ziane Khodja\affil{2},
81 David Laiymani\affil{1},
82 Raphaël Couturier\affil{1} and
83 Arnaud Giersch\affil{1}
88 Femto-ST Institute, DISC Department,
89 University of Franche-Comté,
91 Email:~\email{{charles.ramamonjisoa,david.laiymani,raphael.couturier,arnaud.giersch}@univ-fcomte.fr}\break
93 Department of Aerospace \& Mechanical Engineering,
94 Non Linear Computational Mechanics,
95 University of Liege, Liege, Belgium.
96 Email:~\email{l.zianekhodja@ulg.ac.be}
99 \begin{abstract} %% The behavior of multi-core applications is always a challenge
100 %% to predict, especially with a new architecture for which no experiment has been
101 %% performed. With some applications, it is difficult, if not impossible, to build
102 %% accurate performance models. That is why another solution is to use a simulation
103 %% tool which allows us to change many parameters of the architecture (network
104 %% bandwidth, latency, number of processors) and to simulate the execution of such
105 %% applications. The main contribution of this paper is to show that the use of a
106 %% simulation tool (here we have decided to use the SimGrid toolkit) can really
107 %% help developers to better tune their applications for a given multi-core
110 %% In this paper we focus our attention on the simulation of iterative algorithms to solve sparse linear systems on large clusters. We study the behavior of the widely used GMRES algorithm and two different variants of the Multisplitting algorithms: one using synchronous iterations and another one with asynchronous iterations.
111 %% For each algorithm we have simulated
112 %% different architecture parameters to evaluate their influence on the overall
114 %% The simulations confirm the real results previously obtained on different real multi-core architectures and also confirm the efficiency of the asynchronous Multisplitting algorithm on distant clusters compared to the synchronous GMRES algorithm.
116 The behavior of multi-core applications always proves quite challenging to predict, especially with a new architecture for which no experiment has yet been performed. With some applications, it is difficult, if not impossible, to build accurate performance models. That is why another solution is to use a simulation tool which allows us to change many parameters of the architecture (network bandwidth, latency, number of processors) and to simulate the execution of such applications.
118 In this paper we focus on the simulation of iterative algorithms to solve sparse linear systems. We study the behavior of the GMRES algorithm and two different variants of the multisplitting algorithms: using synchronous or asynchronous iterations. For each algorithm we have simulated different architecture parameters to evaluate their influence on the overall execution time. The simulations confirm the real results previously obtained on different real multi-core architectures and also confirm the efficiency of the asynchronous multisplitting algorithm on distant clusters compared to the GMRES algorithm.
122 %\keywords{Algorithm; distributed; iterative; asynchronous; simulation; simgrid;
124 \keywords{ Performance evaluation, Simulation, SimGrid, Synchronous and asynchronous iterations, Multisplitting algorithms}
128 \section{Introduction} The use of multi-core architectures to solve large
129 scientific problems seems to become imperative in many situations.
130 Whatever the scale of these architectures (distributed clusters, computational
131 grids, embedded multi-core,~\ldots) they are generally well adapted to execute
132 complex parallel applications operating on a large amount of data.
133 Unfortunately, users (industrials or scientists), who need such computational
134 resources, may not have an easy access to such efficient architectures. The cost
135 of using the platform and/or the cost of testing and deploying an application
136 are often very important. So, in this context, it is difficult to optimize a
137 given application for a given architecture. In this way and in order to reduce
138 the access cost to these computing resources it seems very interesting to use a
139 simulation environment. The advantages are numerous: life cycle development,
140 code debugging, ability to obtain results quickly\dots{} In return, the simulation results need to be consistent with the real ones.
142 In this paper we focus on a class of highly efficient parallel algorithms called
143 \emph{iterative algorithms}. The parallel scheme of iterative methods is quite
144 simple. It generally involves the division of the problem into several
145 \emph{blocks} that will be solved in parallel on multiple processing
146 units. Each processing unit has to compute an iteration to send/receive some
147 data dependencies to/from its neighbors and to iterate this process until the
148 convergence of the method. Several well-known studies demonstrate the
149 convergence of these algorithms~\cite{BT89,bahi07}. In this processing mode a
150 task cannot begin a new iteration while it has not received data dependencies
151 from its neighbors. The iteration computation is said to follow a
152 \textit{synchronous} scheme. In the asynchronous scheme a task can compute a new
153 iteration without having to wait for the data dependencies coming from its
154 neighbors. Both communications and computations are \textit{asynchronous}
155 inducing that there is no more idle time, due to synchronizations, between two
156 iterations~\cite{bcvc06:ij}. This model presents some advantages and drawbacks
157 that we detail in Section~\ref{sec:asynchro}. Even if the number of
158 iterations required to converge is generally greater than for the synchronous
159 case, it appears that the asynchronous iterative scheme can significantly
160 reduce overall execution times by suppressing idle times due to
161 synchronizations~(see~\cite{bahi07} for more details).
163 Nevertheless, in both cases (synchronous or asynchronous) it is extremely time
164 consuming to find optimal configurations and deployment requirements for a given
165 application on a given multi-core architecture. Finding good resource
166 allocations policies under varying CPU power, network speeds and loads is very
167 challenging and labor intensive~\cite{Calheiros:2011:CTM:1951445.1951450}. This
168 problematic is even more difficult for the asynchronous scheme where a small
169 parameter variation of the execution platform and of the application data can
170 lead to very different numbers of iterations to reach the convergence and consequently to
171 very different execution times. In this challenging context we think that the
172 use of a simulation tool can greatly leverage the possibility of testing various
175 The {\bf main contribution of this paper} is to show that the use of a
176 simulation tool (i.e. the SimGrid toolkit~\cite{SimGrid}) in the context of real
177 parallel applications (i.e. large linear system solvers) can help developers to
178 better tune their applications for a given multi-core architecture. To show the
179 validity of this approach we first compare the simulated execution of the Krylov
180 multisplitting algorithm with the GMRES (Generalized Minimal RESidual)
181 solver~\cite{saad86} in synchronous mode. The simulation results allow us to
182 determine which method to choose for a given multi-core architecture.
183 Moreover, the obtained results on different simulated multi-core architectures
184 confirm the real results previously obtained on real physical architectures.
185 More precisely the simulated results are in accordance (i.e. with the same order
186 of magnitude) with the works presented in~\cite{couturier15}, which show that
187 the synchronous Krylov multisplitting method is more efficient than GMRES for large
188 scale clusters. Simulated results also confirm the efficiency of the
189 asynchronous multisplitting algorithm compared to the synchronous GMRES
190 especially in case of geographically distant clusters.
192 Thus, with a simple computing architecture (a laptop) SimGrid allows us
193 to run a test campaign of real parallel iterative applications on
194 different simulated multi-core architectures. To our knowledge, there is no
195 related work on the large-scale multi-core simulation of a real synchronous and
196 asynchronous iterative application.
198 This paper is organized as follows. Section~\ref{sec:asynchro} presents the
199 iteration model we use and more particularly the asynchronous scheme. In
200 Section~\ref{sec:simgrid} the SimGrid simulation toolkit is presented.
201 Section~\ref{sec:04} details the different solvers that we use. Finally our
202 experimental results are presented in Section~\ref{sec:expe} followed by some
203 concluding remarks and perspectives.
206 \section{The asynchronous iteration model and the motivations of our work}
209 Asynchronous iterative methods have been studied for many years both theoretically and
210 practically. Many methods have been considered and convergence results have been
211 proved. These methods can be used to solve, in parallel, fixed point problems
212 (i.e. problems for which the solution is $x^\star =f(x^\star)$). In practice,
213 asynchronous iteration methods can be used to solve, for example, linear and
214 non-linear systems of equations or optimization problems. Interested readers are
215 invited to read~\cite{BT89,bahi07}.
217 Before using an asynchronous iterative method, the convergence must be
218 studied. Otherwise, there is no guarantee that the application will reach the convergence. An
219 algorithm that supports both the synchronous or the asynchronous iteration model
220 requires very few modifications to be able to be executed in both variants. In
221 practice, only the communications management and the convergence detection are different. In
222 the synchronous mode, iterations are synchronized, whereas, in the asynchronous
223 one, they are not. It should be noticed that non-blocking communications can be
224 used in both modes. Concerning the convergence detection, synchronous variants
225 can use a global convergence procedure which acts as a global synchronization
226 point. In the asynchronous model, the convergence detection is more tricky as
227 it must not synchronize all the processors. Interested readers can
228 consult~\cite{myBCCV05c,bahi07,ccl09:ij}.
230 The number of iterations required to reach the convergence is generally greater
231 for the asynchronous scheme (this number depends on the delay of the
232 messages). Note that, it is not the case in the synchronous mode where the
233 number of iterations is the same as in the sequential mode. Thus, the
234 set of the parameters of the platform (number of nodes, power of nodes,
235 inter and intra clusters bandwidth and latency,~\ldots) and of the
236 application can drastically change the number of iterations required to get the
237 convergence. It follows that asynchronous iterative algorithms are difficult to
238 optimize since the financial and deployment costs on large scale multi-core
239 architectures are often very important. So, prior to deployment and tests it
240 seems very promising to be able to simulate the behavior of asynchronous
241 iterative algorithms. The problematic is then to show that the results produced
242 by simulation are in accordance with reality (i.e. of the same order of
243 magnitude). To our knowledge, there is no study on this problematic.
248 In the scope of this paper, we have chosen the SimGrid
249 toolkit~\cite{SimGrid,casanova+legrand+quinson.2008.simgrid,casanova+giersch+legrand+al.2014.versatile}
250 to simulate the behavior of parallel iterative linear solvers on different
251 computational grid configurations. In opposite to most of the simulators which
252 are stayed very application-oriented, the SimGrid framework is designed to study
253 the behavior of many large-scale distributed computing platforms as Grids,
254 Peer-to-Peer systems, Clouds or High Performance Computation systems. It is
255 still actively developed by the scientific community and distributed as an open
258 SimGrid provides four user interfaces which can be convenient for different
259 distributed applications. In this paper we are interested on the SMPI
260 (Simulated MPI) user interface which implements about \np[\%]{80} of the MPI 2.0
261 standard~\cite{bedaride+degomme+genaud+al.2013.toward}, and allows minor
262 modifications of the initial code (see Section~\ref{sec:04.02}). SMPI enables
263 the direct simulation of the execution, as in the real life, of an unmodified
264 MPI distributed application, and gets accurate results with the detailed
265 resources consumption.
267 SimGrid simulator uses an XML input file describing the computational grid
268 resources: the number of clusters in the grid, the number of processors/cores in
269 each cluster, the detailed description of the intra and inter networks and the
270 list of the hosts in each cluster (see the details in
271 Section~\ref{sec:expe}). SimGrid employs a fluid model to simulate the use of
272 these resources along the program execution. This model produces accurate
273 results while still running relatively
274 fast~\cite{bedaride+degomme+genaud+al.2013.toward,velho+schnorr+casanova+al.2013.validity}.
275 During the simulation, the computations are really executed, but the communications
276 are intercepted and their execution time evaluated according to the parameters
277 of the simulated platform. It is also possible for SimGrid/SMPI to only keep the
278 duration of large computations by skipping them. Moreover, when applicable, the
279 application can be run by sharing some in-memory structures between the
280 simulated processes and thus allowing the use of very large-scale data.
282 The choice of SimGrid/SMPI as a simulator tool in this study has been emphasized
283 by the results obtained by several studies to validate, in the real
284 environments, the behavior of different network models simulated in
285 SimGrid~\cite{velho+schnorr+casanova+al.2013.validity}. Other studies underline
286 the comparison between the real MPI application executions and the SimGrid/SMPI
287 ones~\cite{guermouche+renard.2010.first,clauss+stillwell+genaud+al.2011.single,bedaride+degomme+genaud+al.2013.toward}. These
288 works show the accuracy of SimGrid simulations compared to the executions on
289 real physical architectures.
291 %% In the scope of this paper, the SimGrid toolkit~\cite{SimGrid,casanova+legrand+quinson.2008.simgrid,casanova+giersch+legrand+al.2014.versatile},
292 %% an open source framework actively developed by its scientific community, has been chosen to simulate the behavior of iterative linear solvers in different computational grid configurations. SimGrid pretends to be non-specialized in opposite to some other simulators which stayed to be very specific oriented-application. One of the well-known SimGrid advantage is its SMPI (Simulated MPI) user interface. SMPI purpose is to execute by simulation in a similar way as in real life, an MPI distributed application and to get accurate results with the detailed resources
293 %% consumption.Several studies have demonstrated the accuracy of the simulation
294 %% compared with execution on real physical architectures. In addition of SMPI,
295 %% Simgrid provides other API which can be convienent for different distrbuted
296 %% applications: computational grid applications, High Performance Computing (HPC),
297 %% P2P but also clouds applications. In this paper we use the SMPI API. It
298 %% implements about \np[\%]{80} of the MPI 2.0 standard and allows minor
299 %% modifications of the initial code~\cite{bedaride+degomme+genaud+al.2013.toward}
300 %% (see Section~\ref{sec:04.02}).
303 %% Provided as an input to the simulator, at least $3$ XML files describe the
304 %% computational grid resources: number of clusters in the grid, number of
305 %% processors/cores in each cluster, detailed description of the intra and inter
306 %% networks and the list of the hosts in each cluster (see the details in Section~\ref{sec:expe}). Simgrid uses a fluid model to simulate the program execution.
307 %% This gives several simulation modes which produce accurate
308 %% results~\cite{bedaride+degomme+genaud+al.2013.toward,
309 %% velho+schnorr+casanova+al.2013.validity}. For instance, the "in vivo" mode
310 %% really executes the computation but "intercepts" the communications (running
311 %% time is then evaluated according to the parameters of the simulated platform).
312 %% It is also possible for SimGrid/SMPI to only keep duration of large
313 %% computations by skipping them. Moreover the application can be run "in vitro"
314 %% by sharing some in-memory structures between the simulated processes and
315 %% thus allowing the use of very large data scale.
318 %% The choice of Simgrid/SMPI as a simulator tool in this study has been emphasized
319 %% by the results obtained by several studies to validate, in real environments,
320 %% the behavior of different network models simulated in
321 %% Simgrid~\cite{velho+schnorr+casanova+al.2013.validity}. Other studies underline
322 %% the comparison between real MPI executions and SimGrid/SMPI
323 %% ones\cite{guermouche+renard.2010.first, clauss+stillwell+genaud+al.2011.single,
324 %% bedaride+degomme+genaud+al.2013.toward}. These works show the accuracy of
325 %% SimGrid simulations.
332 % SimGrid~\cite{SimGrid,casanova+legrand+quinson.2008.simgrid,casanova+giersch+legrand+al.2014.versatile} is a discrete event simulation framework to study the behavior of large-scale distributed computing platforms as Grids, Peer-to-Peer systems, Clouds and High Performance Computation systems. It is widely used to simulate and evaluate heuristics, prototype applications or even assess legacy MPI applications. It is still actively developed by the scientific community and distributed as an open source software.
334 % %%%%%%%%%%%%%%%%%%%%%%%%%
335 % % SimGrid~\cite{SimGrid,casanova+legrand+quinson.2008.simgrid,casanova+giersch+legrand+al.2014.versatile}
336 % % is a simulation framework to study the behavior of large-scale distributed
337 % % systems. As its name suggests, it emanates from the grid computing community,
338 % % but is nowadays used to study grids, clouds, HPC or peer-to-peer systems. The
339 % % early versions of SimGrid date back from 1999, but it is still actively
340 % % developed and distributed as an open source software. Today, it is one of the
341 % % major generic tools in the field of simulation for large-scale distributed
344 % SimGrid provides several programming interfaces: MSG to simulate Concurrent
345 % Sequential Processes, SimDAG to simulate DAGs of (parallel) tasks, and SMPI to
346 % run real applications written in MPI~\cite{MPI}. Apart from the native C
347 % interface, SimGrid provides bindings for the C++, Java, Lua and Ruby programming
348 % languages. SMPI is the interface that has been used for the work described in
349 % this paper. The SMPI interface implements about \np[\%]{80} of the MPI 2.0
350 % standard~\cite{bedaride+degomme+genaud+al.2013.toward}, and supports
351 % applications written in C or Fortran, with little or no modifications (cf Section IV - paragraph B).
353 % Within SimGrid, the execution of a distributed application is simulated by a
354 % single process. The application code is really executed, but some operations,
355 % like communications, are intercepted, and their running time is computed
356 % according to the characteristics of the simulated execution platform. The
357 % description of this target platform is given as an input for the execution, by
358 % means of an XML file. It describes the properties of the platform, such as
359 % the computing nodes with their computing power, the interconnection links with
360 % their bandwidth and latency, and the routing strategy. The scheduling of the
361 % simulated processes, as well as the simulated running time of the application
362 % are computed according to these properties.
364 % To compute the durations of the operations in the simulated world, and to take
365 % into account resource sharing (e.g. bandwidth sharing between competing
366 % communications), SimGrid uses a fluid model. This allows users to run relatively fast
367 % simulations, while still keeping accurate
368 % results~\cite{bedaride+degomme+genaud+al.2013.toward,
369 % velho+schnorr+casanova+al.2013.validity}. Moreover, depending on the
370 % simulated application, SimGrid/SMPI allows to skip long lasting computations and
371 % to only take their duration into account. When the real computations cannot be
372 % skipped, but the results are unimportant for the simulation results, it is
373 % also possible to share dynamically allocated data structures between
374 % several simulated processes, and thus to reduce the whole memory consumption.
375 % These two techniques can help to run simulations on a very large scale.
377 % The validity of simulations with SimGrid has been asserted by several studies.
378 % See, for example, \cite{velho+schnorr+casanova+al.2013.validity} and articles
379 % referenced therein for the validity of the network models. Comparisons between
380 % real execution of MPI applications on the one hand, and their simulation with
381 % SMPI on the other hand, are presented in~\cite{guermouche+renard.2010.first,
382 % clauss+stillwell+genaud+al.2011.single,
383 % bedaride+degomme+genaud+al.2013.toward}. All these works conclude that
384 % SimGrid is able to simulate pretty accurately the real behavior of the
386 %%%%%%%%%%%%%%%%%%%%%%%%%
388 \section{Two-stage multisplitting methods}
390 \subsection{Synchronous and asynchronous two-stage methods for sparse linear systems}
392 In this paper we focus on two-stage multisplitting methods in both their versions (synchronous and asynchronous)~\cite{Frommer92,Szyld92,Bru95}. These iterative methods are based on multisplitting methods~\cite{O'leary85,White86,Alefeld97} and use two nested iterations: the outer iteration and the inner iteration. Let us consider the following sparse linear system of $n$ equations in $\mathbb{R}$:
397 where $A$ is a sparse square and nonsingular matrix, $b$ is the right-hand side and $x$ is the solution of the system. Our work in this paper is restricted to the block Jacobi splitting method. This approach of multisplitting consists in partitioning the matrix $A$ into $L$ horizontal band matrices of order $\frac{n}{L}\times n$ without overlapping (i.e. sub-vectors $\{x_\ell\}_{1\leq\ell\leq L}$ are disjoint). Two-stage multisplitting methods solve the linear system~(\ref{eq:01}) iteratively as follows:
399 x_\ell^{k+1} = A_{\ell\ell}^{-1}(b_\ell - \displaystyle\sum^{L}_{\substack{m=1\\m\neq\ell}}{A_{\ell m}x^k_m}),\mbox{~for~}\ell=1,\ldots,L\mbox{~and~}k=1,2,3,\ldots
402 where $x_\ell$ are sub-vectors of the solution $x$, $b_\ell$ are the sub-vectors of the right-hand side $b$, and $A_{\ell\ell}$ and $A_{\ell m}$ are diagonal and off-diagonal blocks of matrix $A$ respectively. The iterations of these methods can naturally be computed in parallel so that each processor or cluster of processors is responsible for solving one splitting as a linear sub-system:
404 A_{\ell\ell} x_\ell = c_\ell,\mbox{~for~}\ell=1,\ldots,L,
407 where the right-hand sides $c_\ell=b_\ell-\sum_{m\neq\ell}A_{\ell m}x_m$ are computed using the shared vectors $x_m$. In this paper, we use the well-known iterative method GMRES~\cite{saad86} as an inner iteration to approximate the solutions of the different splittings arising from the block Jacobi multisplitting of matrix $A$. The algorithm in Figure~\ref{alg:01} shows the main key points of our block Jacobi two-stage method executed by a cluster of processors. Line~\ref{solve}, the linear sub-system~(\ref{eq:03}) is solved in parallel using the GMRES method where $\MIG$ and $\TOLG$ are the maximum number of inner iterations and the tolerance threshold for GMRES respectively. The convergence of the two-stage multisplitting methods, based on synchronous or asynchronous iterations, has been studied by many authors for example~\cite{Bru95,bahi07}.
410 %\begin{algorithm}[t]
411 %\caption{Block Jacobi two-stage multisplitting method}
412 \begin{algorithmic}[1]
413 \Input $A_\ell$ (sparse matrix), $b_\ell$ (right-hand side)
414 \Output $x_\ell$ (solution vector)\vspace{0.2cm}
415 \State Set the initial guess $x^0$
416 \For {$k=1,2,3,\ldots$ until convergence}
417 \State $c_\ell=b_\ell-\sum_{m\neq\ell}A_{\ell m}x_m^{k-1}$
418 \State $x^k_\ell=Solve_{gmres}(A_{\ell\ell},c_\ell,x^{k-1}_\ell,\MIG,\TOLG)$\label{solve}
419 \State Send $x_\ell^k$ to neighboring clusters\label{send}
420 \State Receive $\{x_m^k\}_{m\neq\ell}$ from neighboring clusters\label{recv}
423 \caption{Block Jacobi two-stage multisplitting method}
428 In this paper, we propose two algorithms of two-stage multisplitting methods. The first algorithm is based on the asynchronous scheme which allows communications to be overlapped by computations and reduces the idle times resulting from the synchronizations. So in the asynchronous mode, our two-stage algorithm uses asynchronous outer iterations and asynchronous communications between clusters. The communications (i.e. lines~\ref{send} and~\ref{recv} in Figure~\ref{alg:01}) are performed by message passing using MPI non-blocking communication routines. The convergence of the asynchronous iterations is detected when all clusters have locally converged:
430 k\geq\MIM\mbox{~or~}\|x_\ell^{k+1}-x_\ell^k\|_{\infty }\leq\TOLM,
433 where $\MIM$ is the maximum number of outer iterations and $\TOLM$ is the tolerance threshold for the two-stage algorithm.
435 The second two-stage algorithm is based on synchronous outer iterations. We propose to use the Krylov iteration based on residual minimization to improve the slow convergence of the multisplitting methods. In this case, a $n\times s$ matrix $S$ is set using solutions issued from the inner iteration:
437 S=[x^1,x^2,\ldots,x^s],~s\ll n.
440 At each $s$ outer iterations, the algorithm computes a new approximation $\tilde{x}=S\alpha$ which minimizes the residual:
442 \min_{\alpha\in\mathbb{R}^s}{\|b-AS\alpha\|_2}.
445 The algorithm in Figure~\ref{alg:02} includes the procedure of the residual minimization and the outer iteration is restarted with a new approximation $\tilde{x}$ at every $s$ iterations. The least-squares problem~(\ref{eq:06}) is solved in parallel by all clusters using the CGLS method~\cite{Hestenes52} sosuch that $\MIC$ is the maximum number of iterations and $\TOLC$ is the tolerance threshold for this method (line~\ref{cgls} in Figure~\ref{alg:02}).
448 %\begin{algorithm}[t]
449 %\caption{Krylov two-stage method using block Jacobi multisplitting}
450 \begin{algorithmic}[1]
451 \Input $A_\ell$ (sparse matrix), $b_\ell$ (right-hand side)
452 \Output $x_\ell$ (solution vector)\vspace{0.2cm}
453 \State Set the initial guess $x^0$
454 \For {$k=1,2,3,\ldots$ until convergence}
455 \State $c_\ell=b_\ell-\sum_{m\neq\ell}A_{\ell m}x_m^{k-1}$
456 \State $x^k_\ell=Solve_{gmres}(A_{\ell\ell},c_\ell,x^{k-1}_\ell,\MIG,\TOLG)$
457 \State $S_{\ell,k\mod s}=x_\ell^k$
459 \State $\alpha = Solve_{cgls}(AS,b,\MIC,\TOLC)$\label{cgls}
460 \State $\tilde{x_\ell}=S_\ell\alpha$
461 \State Send $\tilde{x_\ell}$ to neighboring clusters
463 \State Send $x_\ell^k$ to neighboring clusters
465 \State Receive $\{x_m^k\}_{m\neq\ell}$ from neighboring clusters
468 \caption{Krylov two-stage method using block Jacobi multisplitting}
473 \subsection{Simulation of the two-stage methods using SimGrid toolkit}
476 One of our objectives when simulating the application in SimGrid is, as in real
477 life, to get accurate results (solutions of the problem) but also to ensure the
478 test reproducibility under similar conditions. According to our experience,
479 very few modifications are required to adapt a MPI program for the SimGrid
480 simulator using SMPI (Simulated MPI). The first modification is to include SMPI
481 libraries and related header files (\verb+smpi.h+). The second modification is to
482 suppress all global variables by replacing them with local variables or using a
483 SimGrid selector called "runtime automatic switching"
484 (smpi/privatize\_global\_variables). Indeed, global variables can generate side
485 effects on runtime between the threads running in the same process and generated by
486 SimGrid to simulate the grid environment.
488 \paragraph{Simulation parameters for SimGrid}
489 \ \\ \noindent Before running a SimGrid benchmark, many parameters for the
490 computation platform must be defined. For our experiments, we consider platforms
491 in which several clusters are geographically distant, so that there are intra and
492 inter-cluster communications. In the following, these parameters are described:
495 \item hostfile: hosts description file,
496 \item platform: file describing the platform architecture: clusters (CPU power,
497 \dots{}), intra cluster network description, inter cluster network (bandwidth $bw$,
498 latency $lat$, \dots{}),
499 \item archi : grid computational description (number of clusters, number of
500 nodes/processors in each cluster).
503 In addition, the following arguments are given to the programs at runtime:
506 \item maximum number of inner iterations $\MIG$ and outer iterations $\MIM$,
507 \item inner precision $\TOLG$ and outer precision $\TOLM$,
508 \item matrix sizes of the problem: N$_{x}$, N$_{y}$ and N$_{z}$ on axis $x$, $y$ and $z$ respectively (in our experiments, we solve 3D problem, see Section~\ref{3dpoisson}),
509 \item matrix diagonal value is fixed to $6.0$ for synchronous experiments and $6.2$ for asynchronous ones,
510 \item matrix off-diagonal value is fixed to $-1.0$,
511 \item number of vectors in matrix $S$ (i.e. value of $s$),
512 \item maximum number of iterations $\MIC$ and precision $\TOLC$ for CGLS method,
513 \item maximum number of iterations and precision for the classical GMRES method,
514 \item maximum number of restarts for the Arnorldi process in GMRES method,
515 \item execution mode: synchronous or asynchronous.
518 It should also be noticed that both solvers have been executed with the SimGrid selector \texttt{-cfg=smpi/running\_power} which determines the computational power (here 19GFlops) of the simulator host machine.
520 %%%%%%%%%%%%%%%%%%%%%%%%%
521 %%%%%%%%%%%%%%%%%%%%%%%%%
523 \section{Experimental results}
526 In this section, experiments for both multisplitting algorithms are reported. First the 3D Poisson problem used in our experiments is described.
528 \subsection{The 3D Poisson problem}
530 We use our two-stage algorithms to solve the well-known Poisson problem $\nabla^2\phi=f$~\cite{Polyanin01}. In three-dimensional Cartesian coordinates in $\mathbb{R}^3$, the problem takes the following form:
532 \frac{\partial^2}{\partial x^2}\phi(x,y,z)+\frac{\partial^2}{\partial y^2}\phi(x,y,z)+\frac{\partial^2}{\partial z^2}\phi(x,y,z)=f(x,y,z)\mbox{~in the domain~}\Omega
537 \phi(x,y,z)=0\mbox{~on the boundary~}\partial\Omega
539 where the real-valued function $\phi(x,y,z)$ is the solution sought, $f(x,y,z)$ is a known function and $\Omega=[0,1]^3$. The 3D discretization of the Laplace operator $\nabla^2$ with the finite difference scheme includes 7 points stencil on the computational grid. The numerical approximation of the Poisson problem on three-dimensional grid is repeatedly computed as $\phi=\phi^\star$ such that:
542 \phi^\star(x,y,z)=&\frac{1}{6}(\phi(x-h,y,z)+\phi(x,y-h,z)+\phi(x,y,z-h)\\&+\phi(x+h,y,z)+\phi(x,y+h,z)+\phi(x,y,z+h)\\&-h^2f(x,y,z))
546 until convergence where $h$ is the grid spacing between two adjacent elements in the 3D computational grid.
548 In the parallel context, the 3D Poisson problem is partitioned into $L\times p$
549 sub-problems such that $L$ is the number of clusters and $p$ is the number of
550 processors in each cluster. We apply the three-dimensional partitioning instead
551 of the row-by-row one in order to reduce the size of the data shared at the
552 sub-problems boundaries. In this case, each processor is in charge of
553 parallelepipedic block of the problem and has at most six neighbors in the same
554 cluster or in distant clusters with which it shares data at boundaries.
556 \subsection{Study setup and simulation methodology}
558 First, to conduct our study, we propose the following methodology
559 which can be reused for any grid-enabled applications.\\
561 \textbf{Step 1}: Choose with the end users the class of algorithms or
562 the application to be tested. Numerical parallel iterative algorithms
563 have been chosen for the study in this paper. \\
565 \textbf{Step 2}: Collect the software materials needed for the experimentation.
566 In our case, we have two variants for the resolution of the
567 3D-Poisson problem: (1) using the classical GMRES; (2) using the multisplitting
568 method. In addition, the SimGrid simulator has been chosen to simulate the
569 behaviors of the distributed applications. SimGrid is running in a virtual
570 machine on a simple laptop. \\
572 \textbf{Step 3}: Fix the criteria which will be used for the future
573 results comparison and analysis. In the scope of this study, we retain
574 on the one hand the algorithm execution mode (synchronous and asynchronous)
575 and on the other hand the execution time and the number of iterations to reach the convergence. \\
577 \textbf{Step 4}: Set up the different grid testbed environments that will be
578 simulated in the simulator tool to run the program. The following architectures
579 have been configured in SimGrid : 2$\times$16, 4$\times$8, 4$\times$16, 8$\times$8 and 2$\times$50. The first number
580 represents the number of clusters in the grid and the second number represents
581 the number of hosts (processors/cores) in each cluster. \\
583 \textbf{Step 5}: Conduct extensive and comprehensive testings
584 within these configurations by varying the key parameters, especially
585 the CPU power capacity, the network parameters and also the size of the
588 \textbf{Step 6} : Collect and analyze the output results.
590 \subsection{Factors impacting distributed applications performance in a grid environment}
592 When running a distributed application in a computational grid, many factors may
593 have a strong impact on the performance. First of all, the architecture of the
594 grid itself can obviously influence the performance results of the program. The
595 performance gain might be important theoretically when the number of clusters
596 and/or the number of nodes (processors/cores) in each individual cluster
599 Another important factor impacting the overall performance of the application
600 is the network configuration. Two main network parameters can modify drastically
601 the program output results:
603 \item the network bandwidth ($bw$ in Gbits/s) also known as "the data-carrying
604 capacity" of the network is defined as the maximum amount of data that can transit
605 from one point to another in a unit of time.
606 \item the network latency ($lat$ in microseconds) defined as the delay from the
607 starting time to send a simple data from a source to a destination.
609 Among the network characteristics, another impacting factor is the volume of data exchanged between the nodes in the cluster
610 and between distant clusters. This parameter is application dependent.
612 In a grid environment, it is common to distinguish, on the one hand, the
613 \textit{intra-network} which refers to the links between nodes within a
614 cluster and on the other hand, the \textit{inter-network} which is the
615 backbone link between clusters. In practice, these two networks have
616 different speeds. The intra-network generally works like a high speed
617 local network with a high bandwidth and very low latency. On the contrary, the
618 inter-network connects clusters sometimes via heterogeneous networks components the through internet with a lower speed. The network between distant clusters
619 might be a bottleneck for the global performance of the application.
622 \subsection{Comparison between GMRES and two-stage multisplitting algorithms in
624 In the scope of this paper, our first objective is to analyze
625 when the synchronous Krylov two-stage method has better performances than the
626 classical GMRES method. With a synchronous iterative method, better performances
627 mean a smaller number of iterations and execution time before reaching the
630 Table~\ref{tab:01} summarizes the parameters used in the different simulations:
631 the grid architectures (i.e. the number of clusters and the number of nodes per
632 cluster), the network of inter-clusters backbone links and the matrix sizes of
633 the 3D Poisson problem. However, for all simulations we fix the network
634 parameters of the intra-clusters links: the bandwidth $bw$=10Gbit/s and the latency
635 $lat=8\mu$s. In what follows, we will present the test conditions, the output
636 results and our comments.
642 Grid architecture & 2$\times$16, 4$\times$8, 4$\times$16 and 8$\times$8\\
643 \multirow{2}{*}{Network inter-clusters} & $N1$: $bw$=10Gbit/s, $lat=8\mu$s \\
644 & $N2$: $bw$=1Gbit/s, $lat=50\mu$s \\
645 \multirow{2}{*}{Matrix size} & $Mat1$: N$_{x}\times$N$_{y}\times$N$_{z}$=150$\times$150$\times$150\\
646 & $Mat2$: N$_{x}\times$N$_{y}\times$N$_{z}$=170$\times$170$\times$170 \\ \hline
648 \caption{Parameters for the different simulations}
653 \subsubsection{Simulations for various grid architectures and scaling-up matrix sizes\\}
655 In this section, we analyze the simulations conducted on various grid
656 configurations and for different sizes of the 3D Poisson problem. The parameters
657 of the network between clusters is fixed to $N2$ (see
658 Table~\ref{tab:01}). Figure~\ref{fig:01} shows, for all grid configurations and
659 a given matrix size of 170$^3$ elements, a non-variation in the number of
660 iterations for the classical GMRES algorithm, which is not the case of the
661 Krylov two-stage algorithm. In fact, with multisplitting algorithms, the number
662 of splittings (in our case, it is equal to the number of clusters) influences on the
663 convergence speed. The higher the number of splittings is, the slower the
664 convergence of the algorithm is (see the output results obtained from
665 configurations 2$\times$16 vs. 4$\times$8 and configurations 4$\times$16 vs.
668 The execution times between both algorithms are significant with different grid
669 architectures. The synchronous Krylov two-stage algorithm presents better
670 performances than the GMRES algorithm, even for a high number of clusters (it is about
671 $32\%$ more efficient on a grid of 8$\times$8 than the GMRES). In addition, we can
672 observe a better sensitivity of the Krylov two-stage algorithm (compared to the
673 GMRES one) when scaling up the number of the processors in the computational
674 grid: the Krylov two-stage algorithm is about $48\%$ and the GMRES algorithm is
675 about $40\%$ better on $64$ processors (grid of 8$\times$8) than $32$ processors
676 (grid of 2$\times$16).
680 \includegraphics[width=100mm]{cluster_x_nodes_nx_150_and_nx_170.pdf}
682 \caption{Various grid configurations with two matrix sizes: $150^3$ and $170^3$}
686 \subsubsection{Simulations for two different inter-clusters network speeds\\}
687 In Figure~\ref{fig:02} we present the execution times of both algorithms to
688 solve a 3D Poisson problem of size $150^3$ on two different simulated network
689 $N1$ and $N2$ (see Table~\ref{tab:01}). As previously mentioned, we can see from
690 this figure that the Krylov two-stage algorithm is sensitive to the number of
691 clusters (i.e. it is better to have a small number of clusters). However, we can
692 notice an interesting behavior of the Krylov two-stage algorithm. It is less
693 sensitive to bad network bandwidth and latency for the inter-clusters links than
694 the GMRES algorithms. This means that the multisplitting methods are more
695 efficient for distributed systems with high latency networks.
699 \includegraphics[width=100mm]{cluster_x_nodes_n1_x_n2.pdf}
700 \caption{Various grid configurations with two networks parameters: $N1$ vs. $N2$}
701 %\LZK{CE, remplacer les ``,'' des décimales par un ``.''}
706 \subsubsection{Network latency impacts on performances\\}
707 Figure~\ref{fig:03} shows the impact of the network latency on the performances of both algorithms. The simulation is conducted on a computational grid of 2 clusters of 16 processors each (i.e. configuration 2$\times$16) interconnected by a network of bandwidth $bw$=1Gbit/s to solve a 3D Poisson problem of size $150^3$. According to the results, a degradation of the network latency from $8\mu$s to $60\mu$s implies an absolute execution time increase for both algorithms, but not with the same rate of degradation. The GMRES algorithm is more sensitive to the latency degradation than the Krylov two-stage algorithm.
711 \includegraphics[width=100mm]{network_latency_impact_on_execution_time.pdf}
712 \caption{Network latency impacts on performances}
716 \subsubsection{Network bandwidth impacts on performances\\}
718 Figure~\ref{fig:04} reports the results obtained for the simulation of a grid of
719 $2\times16$ processors interconnected by a network of latency $lat=50\mu$s to
720 solve a 3D Poisson problem of size $150^3$. Increasing the
721 network bandwidth from $1$Gbit/s to $10$Gbit/s results in improving the performances of both algorithms by reducing the execution times. However, the Krylov two-stage
722 algorithm presents a better performance gain in the considered bandwidth
723 interval with a gain of $40\%$ compared to only about $24\%$ for the classical
728 \includegraphics[width=100mm]{network_bandwith_impact_on_execution_time.pdf}
729 \caption{Network bandwith impacts on performances}
733 \subsubsection{Matrix size impacts on performances\\}
735 In these experiments, the matrix size of the 3D Poisson problem varies from
736 $50^3$ to $190^3$ elements. The simulated computational grid is composed of $4$
737 clusters of $8$ processors each interconnected by the network $N2$ (see
738 Table~\ref{tab:01}). As shown in Figure~\ref{fig:05}, the execution
739 times for both algorithms increase with increased matrix sizes. For all problem
740 sizes, the GMRES algorithm is always slower than the Krylov two-stage algorithm.
741 Moreover, for this benchmark, it seems that the greater the problem size is, the
742 bigger the ratio between execution times of both algorithms is. We can also
743 observe that for some problem sizes, the convergence (and thus the execution
744 time) of the Krylov two-stage algorithm varies quite a lot.
745 %This is due to the 3D partitioning of the 3D matrix of the Poisson problem.
746 These findings may greatly help a lot end users to setup the best and the optimal targeted environment for the application deployment when focusing on the problem size scale up.
750 \includegraphics[width=100mm]{pb_size_impact_on_execution_time.pdf}
751 \caption{Problem size impacts on performances}
755 \subsubsection{CPU power impacts on performances\\}
757 Using the SimGrid simulator flexibility, we have tried to determine the impact
758 of the CPU power of the processors in the different clusters on performances of
759 both algorithms. We have varied the CPU power from $1$GFlops to $19$GFlops. The
760 simulation is conducted on a grid of $2\times16$ processors interconnected by
761 the network $N2$ (see Table~\ref{tab:01}) to solve a 3D Poisson problem of size
762 $150^3$. The results depicted in Figure~\ref{fig:06} confirm the performance
763 gain, about $95\%$ for both algorithms, after improving the CPU power of
768 \includegraphics[width=100mm]{cpu_power_impact_on_execution_time.pdf}
769 \caption{CPU Power impacts on performances}
774 To conclude these series of experiments, with SimGrid we have been able to make
775 many simulations with many parameters variations. Doing all these experiments
776 with a real platform is most of the time impossible or very costly. Moreover
777 the behavior of both GMRES and Krylov two-stage algorithms is in accordance
778 with larger real executions on large scale supercomputers~\cite{couturier15}.
781 \subsection{Comparison between synchronous GMRES and asynchronous two-stage multisplitting algorithms}
783 The previous paragraphs put in evidence the interest to simulate the behavior
784 of the application before any deployment in a real environment. In this
785 section, following the same previous methodology, our goal is to compare the
786 efficiency of the multisplitting method in \textit{ asynchronous mode} compared with the
787 classical GMRES in \textit{synchronous mode}.
789 The interest of using an asynchronous algorithm is that there is no more
790 synchronization. With geographically distant clusters, this may be essential.
791 In this case, each processor can compute its iterations freely without any
792 synchronization with the other processors. Thus, an asynchronous algorithm may
793 theoretically reduce the overall execution time and can also improve the algorithm
796 In this section, the SimGrid simulator is used to compare the behavior of the
797 two-stage algorithm in asynchronous mode with GMRES in synchronous mode.
798 Several benchmarks have been performed with various combinations of the grid
799 resources (CPU, Network, matrix size, \ldots). The test conditions are
800 summarized in Table~\ref{tab:02}.
804 %\LZK{Quelle table repporte les gains relatifs?? Sûrement pas Table II !!}
805 %\RCE{Table III avec la nouvelle numerotation}
812 Grid architecture & 2$\times$50 totaling 100 processors\\
813 Processors Power & 1 GFlops to 1.5 GFlops \\
814 \multirow{2}{*}{Network inter-clusters} & $bw$: 5 Mbits to 50 Mbits\\
816 Matrix size & from $62^3$ to $150^3$\\
817 Residual error precision & $10^{-5}$ to $10^{-11}$\\ \hline \\
819 \caption{Test conditions: GMRES in synchronous mode vs. two-stage multisplitting in asynchronous mode}
824 % use the same column width for the following three tables
825 \newlength{\mytablew}\settowidth{\mytablew}{\footnotesize\np{E-11}}
826 \newenvironment{mytable}[1]{% #1: number of columns for data
827 \renewcommand{\arraystretch}{1.3}%
828 \begin{tabular}{|>{\bfseries}r%
829 |*{#1}{>{\centering\arraybackslash}p{\mytablew}|}}}{%
836 % \caption{Relative gain of the multisplitting algorithm compared with the classical GMRES}
841 & 5 & 5 & 5 & 5 & 5 & 50 & 50 & 50 & 50 & 50 \\
844 & 20 & 20 & 20 & 20 & 20 & 20 & 20 & 20 & 20 & 20 \\
847 & 1 & 1 & 1 & 1.5 & 1.5 & 1.5 & 1.5 & 1 & 1.5 & 1.5 \\
850 & 62 & 62 & 62 & 100 & 100 & 110 & 120 & 130 & 140 & 150 \\
853 & \np{E-5} & \np{E-8} & \np{E-9} & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-11}\\
856 & 2.52 & 2.55 & 2.52 & 2.57 & 2.54 & 2.53 & 2.51 & 2.58 & 2.55 & 2.54 \\
860 \caption{Relative gains of the asynchronous two-stage multisplitting algorithm compared to the classical synchronous GMRES algorithm}
865 Table~\ref{tab:03} reports the relative gains between both algorithms. It is
866 defined by the ratio between the execution time of GMRES and the execution time
867 of the multisplitting. The ratio is greater than one because the asynchronous
868 multisplitting version is faster than GMRES. On average, the two-stage
869 multisplitting algorithm is more than $2.5$ times faster than the classical
870 GMRES. These experiments also show the relative tolerance of the multisplitting
871 algorithm when using a low speed network as usually observed with geographically
872 distant clusters through the internet.
876 In this paper we have presented the simulation of the execution of three
877 different parallel solvers on some multi-core architectures. We have shown that
878 the SimGrid toolkit is an interesting simulation tool that has allowed us to
879 determine which method to choose given a specified multi-core architecture.
880 Moreover the simulated results are in accordance (i.e. with the same order of
881 magnitude) with the works presented in~\cite{couturier15}. Simulated results
882 also confirm the efficiency of the asynchronous multisplitting
883 algorithm compared to the synchronous GMRES especially in case of
884 geographically distant clusters.
886 These results are important since it is very time consuming to find optimal
887 configuration and deployment requirements for a given application on a given
888 multi-core architecture. Finding good resource allocation policies under
889 varying CPU power, network speeds and loads is very challenging and labor
890 intensive. This problematic is even more difficult for the asynchronous
891 scheme where a small parameter variation of the execution platform and of the
892 application data can lead to very different numbers of iterations to reach the
893 convergeence and consequently to very different execution times.
896 In future works, we plan to investigate how to simulate the behavior of really
897 large scale applications. For example, if we are interested in simulating the
898 execution of the solvers of this paper with thousand or even dozens of thousands
899 of cores, it is not possible to do that with SimGrid. In fact, this tool will
900 make the real computation. That is why, we plan to focus our research on that problematic.
904 %\section*{Acknowledgment}
906 This work is partially funded by the Labex ACTION program (contract ANR-11-LABX-01-01).
908 \bibliographystyle{wileyj}
909 \bibliography{biblio}
918 %%% ispell-local-dictionary: "american"
921 % LocalWords: Ramamonjisoa Ziane Khodja Laiymani Raphaël Arnaud Giersch Femto
922 % LocalWords: Franche Comté Belfort GMRES multisplitting SimGrid Krylov SMPI