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 \usepackage[american]{babel}
25 % Extension pour les liens intra-documents (tagged PDF)
26 % et l'affichage correct des URL (commande \url{http://example.com})
27 %\usepackage{hyperref}
30 \DeclareUrlCommand\email{\urlstyle{same}}
32 \usepackage[autolanguage,np]{numprint}
34 \renewcommand*\npunitcommand[1]{\text{#1}}
35 \npthousandthpartsep{}}
38 \usepackage[textsize=footnotesize]{todonotes}
40 \newcommand{\AG}[2][inline]{%
41 \todo[color=green!50,#1]{\sffamily\textbf{AG:} #2}\xspace}
42 \newcommand{\RC}[2][inline]{%
43 \todo[color=red!10,#1]{\sffamily\textbf{RC:} #2}\xspace}
44 \newcommand{\LZK}[2][inline]{%
45 \todo[color=blue!10,#1]{\sffamily\textbf{LZK:} #2}\xspace}
46 \newcommand{\RCE}[2][inline]{%
47 \todo[color=yellow!10,#1]{\sffamily\textbf{RCE:} #2}\xspace}
49 \algnewcommand\algorithmicinput{\textbf{Input:}}
50 \algnewcommand\Input{\item[\algorithmicinput]}
52 \algnewcommand\algorithmicoutput{\textbf{Output:}}
53 \algnewcommand\Output{\item[\algorithmicoutput]}
55 \newcommand{\TOLG}{\mathit{tol_{gmres}}}
56 \newcommand{\MIG}{\mathit{maxit_{gmres}}}
57 \newcommand{\TOLM}{\mathit{tol_{multi}}}
58 \newcommand{\MIM}{\mathit{maxit_{multi}}}
61 \usepackage{color, colortbl}
62 \newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
63 \newcolumntype{Z}[1]{>{\raggedleft}m{#1}}
65 \newcolumntype{g}{>{\columncolor{Gray}}c}
66 \definecolor{Gray}{gray}{0.9}
71 \RCE{Titre a confirmer.}
72 \title{Comparative performance analysis of simulated grid-enabled numerical iterative algorithms}
73 %\itshape{\journalnamelc}\footnotemark[2]}
75 \author{ Charles Emile Ramamonjisoa and
78 Lilia Ziane Khodja and
84 Femto-ST Institute - DISC Department\\
85 Université de Franche-Comté\\
87 Email: \email{{raphael.couturier,arnaud.giersch,david.laiymani,charles.ramamonjisoa}@univ-fcomte.fr}
90 %% Lilia Ziane Khodja: Department of Aerospace \& Mechanical Engineering\\ Non Linear Computational Mechanics\\ University of Liege\\ Liege, Belgium. Email: l.zianekhodja@ulg.ac.be
96 \keywords{Algorithm; distributed; iterative; asynchronous; simulation; simgrid; performance}
100 \section{Introduction}
102 \section{The asynchronous iteration model}
106 %%%%%%%%%%%%%%%%%%%%%%%%%
107 %%%%%%%%%%%%%%%%%%%%%%%%%
109 \section{Two-stage splitting methods}
111 \subsection{Multisplitting methods for sparse linear systems}
113 Let us consider the following sparse linear system of $n$ equations in $\mathbb{R}$
118 where $A$ is a sparse square and nonsingular matrix, $b$ is the right-hand side and $x$ is the solution of the system. The multisplitting methods solve the linear system~(\ref{eq:01}) iteratively as follows
120 x^{k+1}=\displaystyle\sum^L_{\ell=1} E_\ell M^{-1}_\ell (N_\ell x^k + b),~k=1,2,3,\ldots
123 where a collection of $L$ triplets $(M_\ell, N_\ell, E_\ell)$ defines the multisplitting of matrix $A$~\cite{O'leary85,White86}, such that: the different splittings are defined as $A=M_\ell-N_\ell$ where $M_\ell$ are nonsingular matrices, and $\sum_\ell{E_\ell=I}$ are diagonal nonnegative weighting matrices and $I$ is the identity matrix. The iterations of the multisplitting methods can naturally be computed in parallel such that each processor or cluster of processors is responsible for solving one splitting as a linear sub-system
125 M_\ell y_\ell = c_\ell^k,\mbox{~such that~} c_\ell^k = N_\ell x^k + b,
128 then the weighting matrices $E_\ell$ are used to compute the solution of the global system~(\ref{eq:01})
130 x^{k+1}=\displaystyle\sum^L_{\ell=1} E_\ell y_\ell.
133 The convergence of the multisplitting methods, based on synchronous or asynchronous iterations, is studied by many authors for example~\cite{O'leary85,bahi97,Bai99,bahi07}. %It is dependent on the condition
135 %\rho(\displaystyle\sum_{\ell=1}^L E_\ell M^{-1}_\ell N_\ell) < 1,
138 %where $\rho$ is the spectral radius of the square matrix.
139 The multisplitting methods are convergent:
141 \item if $A^{-1}>0$ and the splittings of matrix $A$ are weak regular (i.e. $M^{-1}\geq 0$ and $M^{-1}N\geq 0$) when the iterations are synchronous, or
142 \item if $A$ is M-matrix and its splittings are regular (i.e. $M^{-1}\geq 0$ and $N\geq 0$) when the iterations are asynchronous.
144 The solutions of the different linear sub-systems~(\ref{eq:03}) arising from the multisplitting of matrix $A$ can be either computed exactly with a direct method or approximated with an iterative method. In the latter case, the multisplitting methods are called {\it inner-outer iterative methods} or {\it two-stage multisplitting methods}. This kind of methods uses two nested iterations: the outer iteration and the inner iteration (that of the iterative method).
146 In this paper we are focused on two-stage multisplitting methods, in their both versions synchronous and asynchronous, where the well-known iterative method GMRES ({\it Generalized Minimal RESidual})~\cite{saad86} is used as an inner iteration. Furthermore, 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. weighting matrices $E_\ell$ have only zero and one factors). In this case, the iteration of the multisplitting method presented by (\ref{eq:03}) and~(\ref{eq:04}) can be rewritten in the following form
148 A_{\ell\ell} x_\ell^{k+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
151 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. In each outer iteration $k$ until the convergence, each sub-system arising from the block Jacobi multisplitting
153 A_{\ell\ell} x_\ell = c_\ell,
156 is solved iteratively using GMRES method and independently from other sub-systems by a cluster of processors. 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$. Algorithm~\ref{alg:01} shows the main key points of our block Jacobi two-stage method executed by a cluster of processors. In line~\ref{solve}, the linear sub-system~(\ref{eq:06}) is solved in parallel using GMRES method where $\MIG$ and $\TOLG$ are the maximum number of inner iterations and the tolerance threshold of GMRES respectively.
159 \caption{Block Jacobi two-stage multisplitting method}
160 \begin{algorithmic}[1]
161 \Input $A_\ell$ (sparse matrix), $b_\ell$ (right-hand side)
162 \Output $x_\ell$ (solution vector)\vspace{0.2cm}
163 \State Set the initial guess $x^0$
164 \For {$k=1,2,3,\ldots$ until convergence}
165 \State $c_\ell=b_\ell-\sum_{m\neq\ell}A_{\ell m}x_m^{k-1}$
166 \State $x^k_\ell=Solve(A_{\ell\ell},c_\ell,x^{k-1}_\ell,\MIG,\TOLG)$\label{solve}
167 \State Send $x_\ell^k$ to neighboring clusters\label{send}
168 \State Receive $\{x_m^k\}_{m\neq\ell}$ from neighboring clusters\label{recv}
174 Multisplitting methods are more advantageous for large distributed computing platforms composed of hundreds or even thousands of processors interconnected by high latency networks. In this context, the parallel asynchronous model is preferred to the synchronous one to reduce overall execution times of the algorithms, even if it generally requires more iterations to converge. The asynchronous model allows the communications to be overlapped by computations which suppresses the idle times resulting from the synchronizations. So in 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 Algorithm~\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
176 k\geq\MIM\mbox{~or~}\|x_\ell^{k+1}-x_\ell^k\|_{\infty }\leq\TOLM,
179 where $\MIM$ is the maximum number of outer iterations and $\TOLM$ is the tolerance threshold of the two-stage algorithm. The procedure of the convergence detection is implemented as follows. All clusters are interconnected by a virtual unidirectional ring network around which a Boolean token circulates from a cluster to another.
187 \subsection{Simulation of two-stage methods using SimGrid framework}
189 %%%%%%%%%%%%%%%%%%%%%%%%%
190 %%%%%%%%%%%%%%%%%%%%%%%%%
192 \section{Experimental, Results and Comments}
195 \textbf{V.1. Setup study and Methodology}
197 To conduct our study, we have put in place the following methodology
198 which can be reused with any grid-enabled applications.
200 \textbf{Step 1} : Choose with the end users the class of algorithms or
201 the application to be tested. Numerical parallel iterative algorithms
202 have been chosen for the study in the paper.
204 \textbf{Step 2} : Collect the software materials needed for the
205 experimentation. In our case, we have three variants algorithms for the
206 resolution of three 3D-Poisson problem: (1) using the classical GMRES alias Algo-1 in this
207 paper, (2) using the multisplitting method alias Algo-2 and (3) an
208 enhanced version of the multisplitting method as Algo-3. In addition,
209 SIMGRID simulator has been chosen to simulate the behaviors of the
210 distributed applications. SIMGRID is running on the Mesocentre
211 datacenter in Franche-Comte University $[$10$]$ but also in a virtual
214 \textbf{Step 3} : Fix the criteria which will be used for the future
215 results comparison and analysis. In the scope of this study, we retain
216 in one hand the algorithm execution mode (synchronous and asynchronous)
217 and in the other hand the execution time and the number of iterations of
218 the application before obtaining the convergence.
220 \textbf{Step 4 }: Setup up the different grid testbeds environment
221 which will be simulated in the simulator tool to run the program. The
222 following architecture has been configured in Simgrid : 2x16 - that is a
223 grid containing 2 clusters with 16 hosts (processors/cores) each -, 4x8,
224 4x16, 8x8 and 2x50. The network has been designed to operate with a
225 bandwidth equals to 10Gbits (resp. 1Gbits/s) and a latency of 8E-6
226 microseconds (resp. 5E-5) for the intra-clusters links (resp.
227 inter-clusters backbone links).
229 \textbf{Step 5}: Process an extensive and comprehensive testings
230 within these configurations in varying the key parameters, especially
231 the CPU power capacity, the network parameters and also the size of the
232 input matrix. Note that some parameters should be invariant to allow the
233 comparison like some program input arguments.
235 \textbf{Step 6} : Collect and analyze the output results.
237 \textbf{ V.2. Factors impacting distributed applications performance in
240 From our previous experience on running distributed application in a
241 computational grid, many factors are identified to have an impact on the
242 program behavior and performance on this specific environment. Mainly,
243 first of all, the architecture of the grid itself can obviously
244 influence the performance results of the program. The performance gain
245 might be important theoretically when the number of clusters and/or the
246 number of nodes (processors/cores) in each individual cluster increase.
248 Another important factor impacting the overall performance of the
249 application is the network configuration. Two main network parameters
250 can modify drastically the program output results : (i) the network
251 bandwidth (bw=bits/s) also known as "the data-carrying capacity"
252 $[$13$]$ of the network is defined as the maximum of data that can pass
253 from one point to another in a unit of time. (ii) the network latency
254 (lat : microsecond) defined as the delay from the start time to send the
255 data from a source and the final time the destination have finished to
256 receive it. Upon the network characteristics, another impacting factor
257 is the application dependent volume of data exchanged between the nodes
258 in the cluster and between distant clusters. Large volume of data can be
259 transferred in transit between the clusters and nodes during the code
262 In a grid environment, it is common to distinguish in one hand, the
263 "\,intra-network" which refers to the links between nodes within a
264 cluster and in the other hand, the "\,inter-network" which is the
265 backbone link between clusters. By design, these two networks perform
266 with different speed. The intra-network generally works like a high
267 speed local network with a high bandwith and very low latency. In
268 opposite, the inter-network connects clusters sometime via heterogeneous
269 networks components thru internet with a lower speed. The network
270 between distant clusters might be a bottleneck for the global
271 performance of the application.
273 \textbf{V.3 Comparing GMRES and Multisplitting algorithms in
276 In the scope of this paper, our first objective is to demonstrate the
277 Algo-2 (Multisplitting method) shows a better performance in grid
278 architecture compared with Algo-1 (Classical GMRES) both running in
279 \textbf{\textit{synchronous mode}}. Better algorithm performance
280 should mean a less number of iterations output and a less execution time
281 before reaching the convergence. For a systematic study, the experiments
282 should figure out that, for various grid parameters values, the
283 simulator will confirm the targeted outcomes, particularly for poor and
284 slow networks, focusing on the impact on the communication performance
285 on the chosen class of algorithm $[$12$]$.
287 The following paragraphs present the test conditions, the output results
291 \textit{3.a Executing the algorithms on various computational grid
292 architecture scaling up the input matrix size}
297 \begin{tabular}{r c }
299 Grid & 2x16, 4x8, 4x16 and 8x8\\ %\hline
300 Network & N2 : bw=1Gbs-lat=5E-05 \\ %\hline
301 Input matrix size & N$_{x}$ =150 x 150 x 150 and\\ %\hline
302 - & N$_{x}$ =170 x 170 x 170 \\ \hline
307 Table 1 : Clusters x Nodes with NX=150 or NX=170
309 \RCE{J'ai voulu mettre les tableaux des données mais je pense que c'est inutile et ça va surcharger}
312 The results in figure 1 show the non-variation of the number of
313 iterations of classical GMRES for a given input matrix size; it is not
314 the case for the multisplitting method.
316 %\begin{wrapfigure}{l}{60mm}
319 \includegraphics[width=60mm]{cluster_x_nodes_nx_150_and_nx_170.pdf}
320 \caption{Cluster x Nodes NX=150 and NX=170}
325 Unless the 8x8 cluster, the time
326 execution difference between the two algorithms is important when
327 comparing between different grid architectures, even with the same number of
328 processors (like 2x16 and 4x8 = 32 processors for example). The
329 experiment concludes the low sensitivity of the multisplitting method
330 (compared with the classical GMRES) when scaling up to higher input
333 \textit{3.b Running on various computational grid architecture}
337 \begin{tabular}{r c }
339 Grid & 2x16, 4x8\\ %\hline
340 Network & N1 : bw=10Gbs-lat=8E-06 \\ %\hline
341 - & N2 : bw=1Gbs-lat=5E-05 \\
342 Input matrix size & N$_{x}$ =150 x 150 x 150\\ \hline \\
346 %Table 2 : Clusters x Nodes - Networks N1 x N2
347 %\RCE{idem pour tous les tableaux de donnees}
350 %\begin{wrapfigure}{l}{60mm}
353 \includegraphics[width=60mm]{cluster_x_nodes_n1_x_n2.pdf}
354 \caption{Cluster x Nodes N1 x N2}
359 The experiments compare the behavior of the algorithms running first on
360 speed inter- cluster network (N1) and a less performant network (N2).
361 The figure 2 shows that end users will gain to reduce the execution time
362 for both algorithms in using a grid architecture like 4x16 or 8x8: the
363 performance was increased in a factor of 2. The results depict also that
364 when the network speed drops down, the difference between the execution
365 times can reach more than 25\%.
367 \textit{\\\\\\\\\\\\\\\\\\3.c Network latency impacts on performance}
371 \begin{tabular}{r c }
373 Grid & 2x16\\ %\hline
374 Network & N1 : bw=1Gbs \\ %\hline
375 Input matrix size & N$_{x}$ =150 x 150 x 150\\ \hline\\
379 Table 3 : Network latency impact
384 \includegraphics[width=60mm]{network_latency_impact_on_execution_time.pdf}
385 \caption{Network latency impact on execution time}
390 According the results in table and figure 3, degradation of the network
391 latency from 8.10$^{-6}$ to 6.10$^{-5}$ implies an absolute time
392 increase more than 75\% (resp. 82\%) of the execution for the classical
393 GMRES (resp. multisplitting) algorithm. In addition, it appears that the
394 multisplitting method tolerates more the network latency variation with
395 a less rate increase. Consequently, in the worst case (lat=6.10$^{-5
396 }$), the execution time for GMRES is almost the double of the time for
397 the multisplitting, even though, the performance was on the same order
398 of magnitude with a latency of 8.10$^{-6}$.
400 \textit{3.d Network bandwidth impacts on performance}
404 \begin{tabular}{r c }
406 Grid & 2x16\\ %\hline
407 Network & N1 : bw=1Gbs - lat=5E-05 \\ %\hline
408 Input matrix size & N$_{x}$ =150 x 150 x 150\\ \hline
412 Table 4 : Network bandwidth impact
416 \includegraphics[width=60mm]{network_bandwith_impact_on_execution_time.pdf}
417 \caption{Network bandwith impact on execution time}
423 The results of increasing the network bandwidth depict the improvement
424 of the performance by reducing the execution time for both of the two
425 algorithms. However, and again in this case, the multisplitting method
426 presents a better performance in the considered bandwidth interval with
427 a gain of 40\% which is only around 24\% for classical GMRES.
429 \textit{3.e Input matrix size impacts on performance}
433 \begin{tabular}{r c }
436 Network & N2 : bw=1Gbs - lat=5E-05 \\ %\hline
437 Input matrix size & N$_{x}$ = From 40 to 200\\ \hline
441 Table 5 : Input matrix size impact
445 \includegraphics[width=60mm]{pb_size_impact_on_execution_time.pdf}
446 \caption{Pb size impact on execution time}
450 In this experimentation, the input matrix size has been set from
451 Nx=Ny=Nz=40 to 200 side elements that is from 40$^{3}$ = 64.000 to
452 200$^{3}$ = 8.000.000 points. Obviously, as shown in the figure 5,
453 the execution time for the algorithms convergence increases with the
454 input matrix size. But the interesting result here direct on (i) the
455 drastic increase (300 times) of the number of iterations needed before
456 the convergence for the classical GMRES algorithm when the matrix size
457 go beyond Nx=150; (ii) the classical GMRES execution time also almost
458 the double from Nx=140 compared with the convergence time of the
459 multisplitting method. These findings may help a lot end users to setup
460 the best and the optimal targeted environment for the application
461 deployment when focusing on the problem size scale up. Note that the
462 same test has been done with the grid 2x16 getting the same conclusion.
464 \textit{3.f CPU Power impact on performance}
468 \begin{tabular}{r c }
470 Grid & 2x16\\ %\hline
471 Network & N2 : bw=1Gbs - lat=5E-05 \\ %\hline
472 Input matrix size & N$_{x}$ = 150 x 150 x 150\\ \hline
476 Table 6 : CPU Power impact
480 \includegraphics[width=60mm]{cpu_power_impact_on_execution_time.pdf}
481 \caption{CPU Power impact on execution time}
485 Using the SIMGRID simulator flexibility, we have tried to determine the
486 impact on the algorithms performance in varying the CPU power of the
487 clusters nodes from 1 to 19 GFlops. The outputs depicted in the figure 6
488 confirm the performance gain, around 95\% for both of the two methods,
489 after adding more powerful CPU. Note that the execution time axis in the
490 figure is in logarithmic scale.
492 \textbf{V.4 Comparing GMRES in native synchronous mode and
493 Multisplitting algorithms in asynchronous mode}
495 The previous paragraphs put in evidence the interests to simulate the
496 behavior of the application before any deployment in a real environment.
497 We have focused the study on analyzing the performance in varying the
498 key factors impacting the results. In the same line, the study compares
499 the performance of the two proposed methods in \textbf{synchronous mode
500 }. In this section, with the same previous methodology, the goal is to
501 demonstrate the efficiency of the multisplitting method in \textbf{
502 asynchronous mode} compare with the classical GMRES staying in the
505 Note that the interest of using the asynchronous mode for data exchange
506 is mainly, in opposite of the synchronous mode, the non-wait aspects of
507 the current computation after a communication operation like sending
508 some data between nodes. Each processor can continue their local
509 calculation without waiting for the end of the communication. Thus, the
510 asynchronous may theoretically reduce the overall execution time and can
511 improve the algorithm performance.
513 As stated supra, SIMGRID simulator tool has been used to prove the
514 efficiency of the multisplitting in asynchronous mode and to find the
515 best combination of the grid resources (CPU, Network, input matrix size,
516 \ldots ) to get the highest "\,relative gain" in comparison with the
517 classical GMRES time.
520 The test conditions are summarized in the table below :
524 \begin{tabular}{r c }
526 Grid & 2x50 totaling 100 processors\\ %\hline
527 Processors & 1 GFlops to 1.5 GFlops\\
528 Intra-Network & bw=1.25 Gbits - lat=5E-05 \\ %\hline
529 Inter-Network & bw=5 Mbits - lat=2E-02\\
530 Input matrix size & N$_{x}$ = From 62 to 150\\ %\hline
531 Residual error precision: 10$^{-5}$ to 10$^{-9}$\\ \hline
535 Again, comprehensive and extensive tests have been conducted varying the
536 CPU power and the network parameters (bandwidth and latency) in the
537 simulator tool with different problem size. The relative gains greater
538 than 1 between the two algorithms have been captured after each step of
539 the test. Table I below has recorded the best grid configurations
540 allowing a multiplitting method time more than 2.5 times lower than
541 classical GMRES execution and convergence time. The finding thru this
542 experimentation is the tolerance of the multisplitting method under a
543 low speed network that we encounter usually with distant clusters thru the
546 % use the same column width for the following three tables
547 \newlength{\mytablew}\settowidth{\mytablew}{\footnotesize\np{E-11}}
548 \newenvironment{mytable}[1]{% #1: number of columns for data
549 \renewcommand{\arraystretch}{1.3}%
550 \begin{tabular}{|>{\bfseries}r%
551 |*{#1}{>{\centering\arraybackslash}p{\mytablew}|}}}{%
556 \caption{Relative gain of the multisplitting algorithm compared with
558 \label{tab.cluster.2x50}
563 & 5 & 5 & 5 & 5 & 5 & 50 \\
566 & 0.02 & 0.02 & 0.02 & 0.02 & 0.02 & 0.02 \\
569 & 1 & 1 & 1 & 1.5 & 1.5 & 1.5 \\
572 & 62 & 62 & 62 & 100 & 100 & 110 \\
575 & \np{E-5} & \np{E-8} & \np{E-9} & \np{E-11} & \np{E-11} & \np{E-11} \\
578 & 0.396 & 0.392 & 0.396 & 0.391 & 0.393 & 0.395 \\
587 & 50 & 50 & 50 & 50 & 10 & 10 \\
590 & 0.02 & 0.02 & 0.02 & 0.02 & 0.03 & 0.01 \\
593 & 1.5 & 1.5 & 1.5 & 1.5 & 1 & 1.5 \\
596 & 120 & 130 & 140 & 150 & 171 & 171 \\
599 & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-11} & \np{E-5} & \np{E-5} \\
602 & 0.398 & 0.388 & 0.393 & 0.394 & 0.63 & 0.778 \\
611 \section*{Acknowledgment}
614 The authors would like to thank\dots{}
617 \bibliographystyle{wileyj}
618 \bibliography{biblio}
626 %%% ispell-local-dictionary: "american"