Asynchronous Time-Parallel Method based on Laplace Transform

by   Frederic Magoules, et al.

Laplace transform method has proved to be very efficient and easy to parallelize for the solution of time-dependent problems. However, the synchronization delay among processors implies an upper bound on the expectable acceleration factor, which leads to a lot of wasted time. In this paper, we propose an original asynchronous Laplace transform method formalized for quasilinear problems based on the well-known Gaver-Stehfest algorithm. Parallel experiments show the convergence of our new method, as well as several interesting properties compared with the classical algorithms.


page 1

page 2

page 3

page 4


Asynchronous Communications Library for the Parallel-in-Time Solution of Black-Scholes Equation

The advent of asynchronous iterative scheme gives high efficiency to num...

Time-Dependent Wave-Structure Interaction Revisited: Thermo-piezoelectric Scatterers

In this paper, we are concerned with a time-dependent transmission probl...

Convergence analysis of the PML method for time-domain electromagnetic scattering problems

In this paper, a perfectly matched layer (PML) method is proposed to sol...

Asynchronous Parareal Algorithm Applied to European Option Pricing

Asynchronous iterations arise naturally in parallel computing if one wan...

Optimized numerical inverse Laplace transformation

Among the numerical inverse Laplace transformation (NILT) methods, those...

M/D/1 Queues with LIFO and SIRO Policies

While symbolics for the equilibrium M/D/1-LIFO waiting time density are ...

Solving Laplace problems with the AAA algorithm

We present a novel application of the recently developed AAA algorithm t...

1 Introduction

Laplace transforms are powerful tools employed to derive the analytical solutions of differential equations. However, it is too difficult to find or evaluate the analytical inverse transform in closed form. Researchers try to use numerical Laplace transform methods that convert the time-dependent equations to parameter-dependent problems, which could be parallelized easily in frequency domain. There exist a great deal of approaches that could be used to obtain the inverse transform over the past five decades. Davies and Martin in 1979 wrote a survey [13] investigating many promising methods, in which 14 specific algorithms were tested and compared. Duffy [15] concentrated on three methods developed later, in which two methods are based on existing techniques but considerably improved by other researchers. Cohen wrote a comprehensive review book [10] that shows all aspects of Laplace transform inversion. Finally, Kuhlman surveyed five different methods in [23] and their implementations.

There exist many time-parallel methods, including shooting-type methods, space-time domain decomposition methods, space-time multi-grid methods, and direct methods [20]. Obviously, Laplace transform belongs to the direct methods where [12] and [36] are cited as pioneering works in this category. Our paper is based on the algorithm described in [12] that is called Gaver-Stehfest algorithm, which has a good accuracy on a fairly wide range of functions as illustrated in [13], thereafter they developed an iterative scheme in [26, 25] through linearization of implied volatility, leading to a series of iterations throughout the time-stepping process.

Our purpose in this paper is to employ asynchronous iterative scheme in modeling a Laplace-type solver, which might be extremely flexible in data transmission and exploitation throughout the processing phase. Such iterative scheme was first theoretically proposed in [9] for relaxation algorithms, formalized further in [34, 4, 17] with norm-based contraction model and in [5, 6] with nested set model. Another well-known time-parallel method, called Parareal, has been successfully generalized to the asynchronous iterative scheme [31, 29]. In the next section, we recall and present the classical Gaver-Stehfest algorithm. In Section 3, we introduce the iterative Laplace transform method for an option pricing model with variational volatility. In Section 4, we propose a new Laplace method based on asynchronous iterative scheme with some remarks about the its convergence. Finally, some numerical experiments are given in Section 5 and concluding remarks are presented in Section 6.

2 The Gaver-Stehfest algorithm

Given a second-order linear elliptic operator , consider the initial value problem


where the boundary conditions are supposed defined to have a well-posed problem. In the following we note if the space variable does not hamper the formulation and computation. To solve this equation, we employ the Laplace transform defined as


A general contour integral formula for the inverse Laplace transform is given as


where is the Bromwich contour that must be further determined, thus yielding a lot of numerical applications in recovering the original function.

We are interested in the Gaver-Stehfest algorithm, which aims to approximate Equation (3) by a sequence of functions


where are defined as follows


where is an even number denoting the number of processors. This numerical solver could be deduced from Equation (3) by Cauchy integral formula with specific parameters. Notice that

will change sign during iterations. Such behavior is remarkable and indeed affects the numerical performance of Gaver-Stehfest algorithm in some cases, which will be shown in the following sections, thus ignored for the moment.

The implementation of direct Laplace transform method is quite simple. Given the time-dependent problem (1), making Laplace transform yields

where is the Laplace transform of depending on , and obviously

Finally, making inverse Laplace transform yields

which leads to Algorithm 1.

1:Compute using (5)
2:Compute by Laplace transform
3:Reduce and to process 1
4:if rank == 1 then
5:     Compute using (4)
6:end if
Algorithm 1 Direct Laplace transform method for linear equation.

We notice that one should operate a reduction operation after Laplace transform, which is a commonly used term borrowed from message-passing specification, where a designated root node receives data from all nodes and executes an arithmetic or logical operation. Here the specific operation is summation as indicated above. Furthermore, we illustrate the corresponding state diagram in Figure 1.

Figure 1: State diagram of direct Laplace transform method.

Note that Laplace transform method operates in the frequency domain, whereas most of other methods are proceeding in the original space, in which time-dependent terms must be solved by appropriate temporal methods, such as Euler methods, Runge-Kutta methods, and multistep methods.

3 Iterative Laplace transform method for quasilinear problem

It is not desirable to solve quasilinear problems directly and thus we deal with them in a different direction. Throughout this paper, we take Black-Scholes equation as an example to investigate the behavior of iterative methods (see, e.g., [26]).

We note that option pricing is a crucial target in financial decision-making. The breakthrough came with the appearance of the Black-Scholes equation [7], which has a huge influence to the financial market and drives an unexpected prosperity in the trading of derivatives. This equation is often called the Black-Scholes-Merton (BSM) equation since it was further generalized by Merton in important ways [32, 33]. Consider the following European call option pricing problem

where is the option price, depending on stock price and time . As before, we note in the following if possible. Volatility and risk-free interest rate are the constant parameters, with the final and boundary conditions given by


where is the maturity of the option and is the strike price. If there exist transaction costs, then the price prediction may become much more complex. The volatility can be treated in different ways, for instance using a modified volatility function . Here we assume that the price of option is a parameter of , thus leading to the BSM equation with implied volatility


There are many ways to choose an appropriate . Here we adopt the function described in [26] as following


where denotes constant historical volatility.

We should apply Laplace transform method to Equation (7), but this one is too complex to be employed directly. Moreover, from (6) we notice that (7) is a backward equation. Hence, we perform the following variable transformation

substituting into Equation (7) gives


with the corresponding conditions


Implied volatility defined in (8) becomes


Notice that (9) is forward parabolic and the fractional term is flexible for both linear and quasilinear model. If we import a constant volatility, the fractional term will vanish; if we use implied volatility like the one defined (11), the constant coefficient therein will also disappear.

This equation, however, depends on a quasilinear coefficient that should be tackled before getting forward. We share the idea of [26] by introducing a one-step retard into , which linearizes the quasilinear equation (9) in the form


We note that can be seen as a “frozen variable” which is assigned the previous value of and thus does not depend on . Now we prepare to apply the Laplace transform method to Equation (12). Setting



More specifically, this yields an iterative scheme illustrated in Algorithm 2. The same algorithm has been investigated in [26] under the name of “iterative coefficient–inverse Laplace transform”.

2:Compute using (5)
5:     Compute by Laplace transform
6:     Reduce and to process 1
7:     if rank == 1 then
8:         Compute using (4)
9:     end if
10:     Broadcast to other processes
Algorithm 2 Synchronous Laplace transform method for quasilinear equation.

Since all processors need the values of at the beginning of each iteration, we should broadcast such buffer to all neighbors. Note that “broadcast” is a collective communication operation which involves a specific root node communicating its data to all other nodes. We observe that is computed outside because it depends only on the rank and the number of processors . Similarly, the state diagram is shown in Figure 2.

Figure 2: State diagram of synchronous Laplace transform method.

We mention here that there are large numbers of models dealing with the transaction cost, such as those in [8] and [3], some of which are nonlinear differential equations that have different forms and properties. There exist meanwhile many other types of equations related to the option pricing problem, such as the martingale pricing theory [21], the binomial options pricing model [11], and the stochastic volatility option models [22, 14]. They are completely different from the Black-Scholes model and thus should be addressed separately. In general, one might employ Algorithm 2 to other types of equations generalized from Equation (1), such as

Substituting the frozen variable and applying Laplace transform yields

where linear operator can be solved by any appropriate discretization methods.

4 Asynchronous Laplace transform method

4.1 Formalization

Asynchronous iterative method releases the restriction of strict data dependency whereby iterations are executed by several processors in arbitrary order without any synchronization during computation, first formalized in [9] for linear systems. Obviously, it must be subject to some constraints that each processor still needs to continuously obtain the latest information. Since there is no synchronization imposed to such scheme, asynchronous iterative algorithm may exhibit incredible flexibility and efficiency even with the expansion of computational nodes, which contributes to overcome the fault tolerance and the load balancing problems.

Consider a vector space



Practically, denotes the number of processors. Hence

leading to the classical parallel iterative scheme


which is presented in Figure 3.

Figure 3: Synchronous iterations with asynchronous communications.

In order to compute th iteration’s in processor , one needs to collect all data of iteration from other processors, which imposes a synchronization point on the computational framework at the end of each iteration.

Now we define the integer subsets such that

and let be nonnegative integer-valued functions satisfying


with . We introduce here an iterative scheme that generates the following sequence


which is actually the mathematical model of asynchronous iterative scheme and illustrated in Figure 4.

Figure 4: Asynchronous iterations with asynchronous communications.

Generally, for the sake of convergence analysis, Assumption 1 and Assumption 2 are armed therewith, which ensure that processors read eventually the latest information for each element from essential neighbors and no processor stops updating during iterations.

Assumption 1.


Assumption 2.


Now we apply asynchronous iterative scheme to the iterative Laplace transform method depicted in Algorithm 2. Since newer value depends on the previous value , thereafter we use and to denote these two values during iterations. let


and for a fixed time span , we define


Notice that

which leads to a handy operator


Combining (17), (18), and (19) gives




which follows exactly the classical parallel iterative scheme (14). Hence, iterative Laplace method can be generalized naturally to the asynchronous iterative scheme subject to (15) and assumed to satisfy Assumption 1 and Assumption 2


then leading to Algorithm 3.

2:Compute using (5)
7:     for all   do
8:         Request receiving from process
9:     end for
10:     Compute by Laplace transform
11:     for all   do
12:         Request sending to process
13:     end for
14:     Compute using (4)
Algorithm 3 Asynchronous Laplace transform method for quasilinear equation.

The most remarkable notation therein is that exhibits the chaotic behavior of asynchronous iterative scheme, where intrinsically retard term and execution set play important roles during iterations. Meanwhile, another hinge consists in the sending and receiving operations without synchronization point, which yields the asynchronous performance, as opposed to the aforementioned blocking reduction and broadcast operations. Note that one may also employ nonblocking reduction and broadcast operations to implement Algorithm 3. The state diagram is depicted in Figure 5.

Figure 5: State diagram of asynchronous Laplace transform method.

4.2 Remarks on convergence

To the best of our knowledge, unfortunately, there is no available theory that can be used to prove the convergence of asynchronous Laplace algorithm. According to the norm-based convergence theory (see [17], which is based on the results in [34, 4]), one should prove a relationship in the form

where and denotes the weighted maximum norm with , in order to establish the convergence result. This imposes a contraction condition to the operator . We recall that in our case there exists a summation operator , as seen in (20), and thus is not an assembly but a summation of subvectors. It is clear that the summation operator could not lead to desired partial order under contraction. A more general but similar idea proposed in [5] (see also [6]) is based on the nested set theory, in which a technique known as “box condition” allows to obtain convergence result without using any norm. This can be written as

where for all . This condition involves a superposition of infinite subsets and thus suffers from the same problem as the previous approach. In addition, two-stage and nonstationary models asynchronous models were discussed in [18, 19] and some papers made use of the partial ordering such as [16, 35]. These models have gained much popularity in some problems but do not correspond to our case.

The expansion of summation operator is indeed compensated by the frequency decomposition, which may contribute to the convergence. In what follows we will not pursue the convergence issue further; instead, we will focus on the numerical behavior of our method.

5 Numerical experiments

5.1 Environmental setting

We have conducted several experiments for the direct Laplace transform method, using the Gaver-Stehfest algorithm as illustrated in Algorithm 1, as well as the iterative variants based on Algorithms 2 and 3. All tests are executed using (13) with in the linear case, whereas (11) is adopted in the quasilinear case.

All tests are launched on an SGI ICE X cluster connected with InfiniBand, involving MPI library to run parallel applications, which is supported by SGI-MPT 2.14. Mathematical operations and linear systems solvers are implemented by the Alinea programming library [27] and iterative algorithms are programmed by JACK [28, 30] for both synchronous and asynchronous variants.

5.2 Direct scheme for linear equation

According to [37], number of processors must be even. We first illustrate the numerical results of Algorithm 1 applied to linear equation, shown in Figure 6.

Figure 6: Direct Laplace transform method for linear equation: absolute error (left), relative error (right).

Here we run programs with different numbers of processes and observe that error behaves like a convex function for each specific time span, where denoting the maturity of the specific option contract. We choose as convergence threshold, which is adequate and conspicuous for the financial data. Table 1 gives the convergence interval of direct method.

convergence interval of remarks
1 8, 10, 12, 14, 16, 18 : inaccurate; : inaccurate
5 10, 12, 14, 16, 18 : inaccurate; : inaccurate
10 14, 16, 18 : inaccurate; : inaccurate
20 inaccurate, oscillating
Table 1: Convergence zone of direct Laplace transform method with large maturities.

Here, convergence interval and convergence zone are defined as the set of that leads to convergent result. Notice that the larger we assign for time span, the narrower convergence interval we get. When , there is no convergence interval observed throughout experiments. In this case, we should compute results step-by-step, as described in [26]. Here “oscillating” indicates that we could not observe a contracting behavior when increases monotonously. On the contrary, results fall into an inconsistent solution range and oscillate divergently around initial point. We mention here that large , denoting the maturity with the unit of year, is not a normal test case in option pricing problem, which is shown only for the experimental purpose. In the following tests, we will use primarily small to illustrate the real behaviors.

Davies and Martin [13] mentioned that the Laplace transform method based on the Gaver-Stehfest algorithm gives good accuracy on a fairly wide range of functions. Figure 6 and Table 1 show that a well-posed option pricing problem leads to good accuracy with a fairly wide range of convergence interval (e.g., lead to a convergence interval ). Theoretically, the more processors we use, the more precision we obtain through parallel processing [24]. In practice, however, we must consider the arithmetical precision of the test machine, namely, we decrease the truncation error thanks to the large number of processors and meanwhile increase the rounding error due to digit width limit. Unfortunately, the side effect is dramatically large due to the rapid growth of with sign alternating. Thus we could not exploit the power of large scale parallel computing, which acts as an intrinsic defect of the Gaver-Stehfest algorithm.

5.3 Iterative scheme for quasilinear equation

Now we focus on the numerical behavior of iterative Laplace transform methods illustrated in Algorithms 2 and 3. We continue to use the convergence threshold as the termination condition, which may lead to an infinite loop in our program. Experimental results are shown in Tables 2 and 3.

convergence interval of remarks
0.01 6, 8, 10, 12 : inaccurate; : divergent
0.1 4, 6, 8, 10, 12 : inaccurate; : divergent
1 4, 6, 8, 10, 12 : inaccurate; : divergent
Table 2: Convergence zone of synchronous Laplace transform method.
convergence interval of remarks
0.01 6 : inaccurate; : divergent
0.1 6 : inaccurate; : divergent
1 4, 6 : inaccurate; : divergent
Table 3: Convergence zone of asynchronous Laplace transform method.

We notice that the convergence interval of the asynchronous Laplace transform method is much narrower than the synchronous method. This is a reasonable compensation because asynchronous methods require less waiting time in communication but more computational iterations than classical parallel scheme, which amplify the unstable issue of with rapid growth and sign alternating.

On the other hand, for a large , both synchronous and asynchronous algorithms become divergent. This means that we could not obtain a result underneath the given threshold since the uncertainty propagates throughout iterations. Hence, such behavior is caused by an intrinsic issue of the Gaver-Stehfest algorithm, which is the same as that in the direct method.

We illustrate in Table 4 the results of direct Laplace method with the same interval of as in Tables 2 and 3.

convergence interval of remarks
0.01 4, 6, 8, 10, 12, 14, 16, 18, 20, 22 : inaccurate; : inaccurate
0.1 6, 8, 10, 12, 14, 16, 18, 20 : inaccurate; : inaccurate
1 6, 8, 10, 12, 14, 16, 18 : inaccurate; : inaccurate
Table 4: Convergence zone of direct Laplace transform method with small maturities.

One should keep in mind that the iterative results are obtained from the BSM equation with implied volatility, whereas the direct method, as seen in Algorithm 2, is applied to the linear equation where is a constant. Thus, comparing the former with the latter could not lead to an obvious conclusion for the performance of our method. We show these results only for clarity and completeness.

Figure 7 illustrates another experiment that divides temporal interval into several smaller parts, as described in [26].

Figure 7: An example of successive asynchronous Laplace iteration with , , : first time interval (left), last time interval (right).

As mentioned before, we address primarily the real situation in option pricing problem where is small. For a small we need no more subdivision in the time domain. Thus, we conduct this test only for evaluating the performance of the algorithm illustrated in [26], where we execute asynchronous Laplace transform method step-by-step along time domain. Similar to the classical temporal integration schemes like Euler methods, we perform computations in each step using always the latest data, but employ a rather wide subinterval to fit the characteristic of the Laplace transform method. The number of processors, the step size, and the number of steps are chosen as , , , respectively. We can see that the execution process in the first time interval is highly structured, except for some unordered data due to the asynchronous scheme. Nonetheless, the last time interval gives different behavior that seems somewhat affected by the unfinished chores due to the asynchronous iterative scheme. This produces larger retards and converges more quickly when receiving the latest data in some steps.

Finally, we test the accuracy of the asynchronous Laplace method in comparison with the synchronous version. We fix the volatility , risk-free rate , maturity , and number of processor . Recall that

The constant volatility and the strike price in the right hand side are fixed, whereas the volatility changes with . Then we change the stoke price and strike price to compare the two methods. Results are shown in Table 5.

Table 5: Synchronous and asynchronous Laplace results with synchronous option prices , asynchronous option prices , absolute errors , and relative errors , given , , , and .

Here, and are the option prices obtained from synchronous and asynchronous iterations, respectively. and act as absolute and relative errors such that

We notice that the asynchronous option prices are close enough to the synchronous cases. This observation achieves our expectations, as we speculated that the chaotic iterative process does not affect dramatically the final result. Since the synchronous Laplace transform method has proved to be accurate in predicting the option price values using the BSM equation with implied volatility, we can draw a conclusion that our method yields accurate results as well.

5.4 Further discussion

We have seen that Laplace transform methods are highly influenced by the number of digits of machine precision that are prone to round-off error propagation, which leads to an unstable behavior. Notice that (5) increases quickly with commutative plus-minus sign. The accuracy mounts in the beginning as the number of processes increases and then decreases rapidly.

One could avoid such effect for Algorithm 1 by simply adopting the recommended precision in a multi-precision computational environment as discussed in [1]. Meanwhile, the fixed Talbot method [38] requires the precision , which performs a less restrictive behavior and seems a better alternative to the Gaver-Stehfest algorithm according to the experiments therein. Further investigation gives a unified framework for various types of inverse Laplace transform algorithms [2]. Recall that the contour integral is

Let , we have

Then taking


Finally, using Cauchy integral formula, this becomes



where is an even integer, and notice that , we get again the Gaver-Stehfest formula (4). We could further rewrite (22) in the form


Combining (17) and (19), we could obtain the asynchronous scheme (21) for the unified Laplace transform model, whereby the methods discussed in [2], such as Fourier series methods and Talbot method, might also be executed in asynchronous mode, although no rigorous proofs are available. The contour integral approach with quadrature is another choice and has proved to be very efficient for parabolic problems with time-independent coefficients [36]. One may expect to combine these strategies to tackle the quasilinear equations in a more effective way.

These issues are beyond the scope of this paper. We leave this work, as well as a rigorous proof of the asynchronous convergence, as a direction of future investigation.

6 Concluding remarks

In this paper, we investigated the Laplace transform and proposed a new method based on the asynchronous iterative scheme. We formalized the mathematical model and gave some remarks on its convergence. Numerical experiments showed that asynchronous Laplace transform method converges chaotically to the correct solution sets. The intrinsic numerical approximation problem affects the accuracy of Gaver-Stehfest algorithm and led to a narrow convergence interval. There exists much remains to be done to prove the convergence of this method, which might require a complete extension of the existing asynchronous theory. Further investigation could focus on the asynchronous formalization for other efficient inversion formulas and nonlinear solvers.


  • [1] J. Abate and P. P. Valkó. Multi-precision Laplace transform inversion. Int. J. Numer. Methods Eng., 60(5):979–993, 2004.
  • [2] J. Abate and W. Whitt. A unified framework for numerically inverting Laplace transforms. INFORMS J. Comput., 18(4):408–421, 2006.
  • [3] G. Barles and H. M. Soner. Option pricing with transaction costs and a nonlinear Black-Scholes equation. Financ. Stoch., 2(4):369–397, 1998.
  • [4] G. M. Baudet. Asynchronous iterative methods for multiprocessors. J. ACM, 25(2):226–244, 1978.
  • [5] D. P. Bertsekas. Distributed asynchronous computation of fixed points. Math. Program., 27(1):107–120, 1983.
  • [6] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed Computation: Numerical Methods. Prentice-Hall, 1989.
  • [7] F. Black and M. Scholes. The pricing of options and corporate liabilities. J. Political Econ., 81(3):637–654, 1973.
  • [8] P. P. Boyle and T. Vorst. Option replication in discrete time with transaction costs. J. Financ., 47(1):271–293, 1992.
  • [9] D. Chazan and W. Miranker. Chaotic relaxation. Linear Algebra Appl., 2(2):199–222, 1969.
  • [10] A. M. Cohen. Numerical Methods for Laplace Transform Inversion. Springer, 2007.
  • [11] J. C. Cox, S. A. Ross, and M. Rubinstein. Option pricing: A simplified approach. J. Financ. Econ., 7(3):229–263, 1979.
  • [12] D. Crann, A. J. Davies, C.-H. Lai, and S. H. Leong. Time domain decomposition for European options in financial modelling. Contemp. Math., 218:486–491, 1998.
  • [13] B. Davies and B. Martin. Numerical inversion of the Laplace transform: a survey and comparison of methods. J. Comput. Phys., 33(1):1–32, 1979.
  • [14] J.-C. Duan. The GARCH option pricing model. Math. Financ., 5(1):13–32, 1995.
  • [15] D. G. Duffy. On the numerical inversion of Laplace transforms: Comparison of three new methods on characteristic problems from applications. ACM Trans. Math. Softw., 19(3):333–359, 1993.
  • [16] D. El Baz, P. Spitéri, J.-C. Miellou, and D. Gazen. Asynchronous iterative algorithms with flexible communication for nonlinear network flow problems. J. Parallel Distrib. Comput., 38(1):1–15, 1996.
  • [17] M. N. El Tarazi. Some convergence results for asynchronous algorithms. Numer. Math., 39(3):325–340, 1982. (in French).
  • [18] A. Frommer and D. B. Szyld. Asynchronous two-stage iterative methods. Numer. Math., 69(2):141–153, 1994.
  • [19] A. Frommer and D. B. Szyld. On asynchronous iterations. J. Comput. Appl. Math., 123(1-2):201–216, 2000.
  • [20] M. J. Gander. 50 years of time parallel time integration. In T. Carraro, M. Geiger, S. Körkel, and R. Rannacher, editors, Multiple Shooting and Time Domain Decomposition Methods, pages 69–113. Springer, Cham, 2015.
  • [21] J. M. Harrison and D. M. Kreps. Martingales and arbitrage in multiperiod securities markets. J. Econ. Theory, 20(3):381–408, 1979.
  • [22] J. Hull and A. White. The pricing of options on assets with stochastic volatilities. J. Financ., 42(2):281–300, 1987.
  • [23] K. L. Kuhlman. Review of inverse Laplace transform algorithms for Laplace-space numerical approaches. Numer. Algorithms, 63(2):339–355, 2013.
  • [24] A. Kuznetsov. On the convergence of the Gaver–Stehfest algorithm. SIAM J. Numer. Anal., 51(6):2984–2998, 2013.
  • [25] C.-H. Lai. On transformation methods and the induced parallel properties for the temporal domain. In F. Magoulès, editor, Substructuring Techniques and Domain Decomposition Methods, chapter 3, pages 45–70. Saxe-Coburg Publications, Stirlingshire, UK, 2010.
  • [26] C.-H. Lai, A. K. Parrott, S. Rout, and M. E. Honnor. A distributed algorithm for European options with nonlinear volatility. Comput. Math. Appl., 49(5):885–894, 2005.
  • [27] F. Magoulès and A.-K. Cheik Ahamed. Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units. Int. J. High Perform. Comput. Appl., 29(3):284–310, 2015.
  • [28] F. Magoulès and G. Gbikpi-Benissan. JACK: An asynchronous communication kernel library for iterative algorithms. J. Supercomput., 73(8):3468–3487, 2017.
  • [29] F. Magoulès and G. Gbikpi-Benissan.

    Asynchronous Parareal time discretization for partial differential equations.

    SIAM J. Sci. Comput., 40(6):C704–C725, 2018.
  • [30] F. Magoulès and G. Gbikpi-Benissan. JACK2: An MPI-based communication library with non-blocking synchronization for asynchronous iterations. Adv. Eng. Softw., 119:116–133, 2018.
  • [31] F. Magoulès, G. Gbikpi-Benissan, and Q. Zou. Asynchronous iterations of Parareal algorithm for option pricing models. Mathematics, 6(4):1–18, 2018.
  • [32] R. C. Merton. Theory of rational option pricing. Bell J. Econ., 4(1):141–183, 1973.
  • [33] R. C. Merton. On the pricing of corporate debt: the risk structure of interest rates. J. Financ., 29(2):449–470, 1974.
  • [34] J.-C. Miellou. Algorithmes de relaxation chaotique à retards. ESAIM: M2AN, 9(R1):55–82, 1975. (in French).
  • [35] J.-C. Miellou, D. El Baz, and P. Spitéri. A new class of asynchronous iterative algorithms with order intervals. Math. Comput., 67(221):237–255, 1998.
  • [36] D. Sheen, I. H. Sloan, and V. Thomée. A parallel method for time-discretization of parabolic problems based on contour integral representation and quadrature. Math. Comput., 69(229):177–195, 2000.
  • [37] H. Stehfest. Numerical inversion of Laplace transforms. Commun. ACM, 13(1):47–49, 1970.
  • [38] A. Talbot. The accurate numerical inversion of Laplace transforms. IMA J. Appl. Math., 23(1):97–120, 1979.