The Waveform Relaxation method for systems of differential/algebraic equations

This paper reports efforts towards establishing a parallel numerical algorithm known as Waveform Relaxation (WR) for simulating large systems of differential/algebraic equations. The WR algorithm was established as a relaxation based iterative method for the numerical integration of systems of ODEs over a finite time interval. In the WR approach, the system is broken into subsystems which are solved independently, with each subsystem using the previous iterate waveform as “guesses” about the behavior of the state variables in other subsystems. Waveforms are then exchanged between subsystems, and the subsystems are then resolved repeatedly with this improved information about the other subsystems until convergence is achieved.

In this paper, a WR algorithm is introduced for the simulation of generalized high-index DAE systems. As with ODEs, DAE systems often exhibit a multirate behavior in which the states vary as differing speeds. This can be exploited by partitioning the system into subsystems as in the WR for ODEs. One additional benefit of partitioning the DAE system into subsystems is that some of the resulting subsystems may be of lower index and, therefore, do not suffer from the numerical complications that high-index systems do. These lower index subsystems may therefore be solved by less specialized simulations. This increases the efficiency of the simulation since only a portion of the problem must be solved with specially tailored code. In addition, this paper established solvability requirements and convergence theorems for varying index DAE systems for WR simulation.     

Parallel direct methods for solving the system of linear equations with pipelining on a multicore using OpenMP

Recent developments in high performance computer architecture have a significant effect on all fields of scientific computing. Linear algebra and especially the solution of linear systems of equations lie at the heart of many applications in scientific computing. This paper describes and analyzes three parallel versions of the dense direct methods such as the Gaussian elimination method and the LU form of Gaussian elimination that are used in linear system solving on a multicore using an OpenMP interface. More specifically, we present two naive parallel algorithms based on row block and row cyclic data distribution and we put special emphasis on presenting a third parallel algorithm based on the pipeline technique. Further, we propose an implementation of the pipelining technique in OpenMP. Experimental results on a multicore CPU show that the proposed OpenMP pipeline implementation achieves good overall performance compared to the other two naive parallel methods. Finally, in this work we propose a simple, fast and reasonably analytical model to predict the performance of the direct methods with the pipelining technique.     

Block s‐step Krylov iterative methods

Block (including s‐step) iterative methods for (non)symmetric linear systems have been studied and implemented in the past. In this article we present a (combined) block s‐step Krylov iterative method for nonsymmetric linear systems. We then consider the problem of applying any block iterative method to solve a linear system with one right‐hand side using many linearly independent initial residual vectors. We present a new algorithm which combines the many solutions obtained (by any block iterative method) into a single solution to the linear system. This approach of using block methods in order to increase the parallelism of Krylov methods is very useful in parallel systems. We implemented the new method on a parallel computer and we ran tests to validate the accuracy and the performance of the proposed methods. It is expected that the block s‐step methods performance will scale well on other parallel systems because of their efficient use of memory hierarchies and their reduction of the number of global communication operations over the standard methods. Copyright © 2009 John Wiley & Sons, Ltd.     

Parallel iterative methods using factorized preconditioning matrices for solving elliptic equations on triangular grids

Parallel analogs of the variants of the incomplete Cholesky-conjugate gradient method and the modified incomplete Cholesky-conjugate gradient method for solving elliptic equations on uniform triangular and unstructured triangular grids on parallel computer systems with the MIMD architecture are considered. The construction of parallel methods is based on the use of various variants of ordering the grid points depending on the decomposition of the computation domain. Results of the theoretic and experimental studies of the convergence rate of these methods are presented. The solution of model problems on a moderate number processors is used to examine the efficiency of the proposed parallel methods.     

Waveform Relaxation Methods for Lie-Group Equations

In this paper, we derive and analyse waveform relaxation (WR) methods for solving differential equations evolving on a Lie-group. We present both continuous-time and discrete-time WR methods and study their convergence properties. In the discrete-time case, the novel methods are constructed by combining WR methods with Runge-Kutta-Munthe-Kaas (RK-MK) methods. The obtained methods have both advantages of WR methods and RK-MK methods, which simplify the computation by decoupling strategy and preserve the numerical solution of Lie-group equations on a manifold. Three numerical experiments are given to illustrate the feasibility of the new WR methods.     

Computing periodic solutions of linear differential-algebraic equations by waveform relaxation

We propose an algorithm, which is based on the waveform relaxation (WR) approach, to compute the periodic solutions of a linear system described by differential-algebraic equations. For this kind of two-point boundary problems, we derive an analytic expression of the spectral set for the periodic WR operator. We show that the periodic WR algorithm is convergent if the supremum value of the spectral radii for a series of matrices derived from the system is less than 1. Numerical examples, where discrete waveforms are computed with a backward-difference formula, further illustrate the correctness of the theoretical work in this paper.

     

Waveform relaxation methods for implicit differential equations

We apply a Runge-Kutta-based waveform relaxation method to initial-value problems for implicit differential equations. In the implementation of such methods, a sequence of nonlinear systems has to be solved iteratively in each step of the integration process. The size of these systems increases linearly with the number of stages of the underlying Runge-Kutta method, resulting in high linear algebra costs in the iterative process for high-order Runge-Kutta methods. In our earlier investigations of iterative solvers for implicit initial-value problems, we designed an iteration method in which the linear algebra costs are almost independent of the number of stages when implemented on a parallel computer system. In this paper, we use this parallel iteration process in the Runge-Kutta waveform relaxation method. In particular, we analyse the convergence of the method. The theoretical results are illustrated by a few numerical examples.     

Almost block diagonal linear systems: sequential and parallel solution techniques,and applications

Almost block diagonal (ABD) linear systems arise in a variety of contexts, specifically in numerical methods for two‐point boundary value problems for ordinary differential equations and in related partial differential equation problems. The stable, efficient sequential solution of ABDs has received much attention over the last fifteen years and the parallel solution more recently. We survey the fields of application with emphasis on how ABDs and bordered ABDs (BABDs) arise. We outline most known direct solution techniques, both sequential and parallel, and discuss the comparative efficiency of the parallel methods. Finally, we examine parallel iterative methods for solving BABD systems. Copyright © 2000 John Wiley & Sons, Ltd.     

Parallel iteration of the extended backward differentiation formulas

The extended backward differentiation formulas (EBDFs) and theirmodified form (MEBDF) were proposed by Cash in the 1980s forsolving initial value problems (IVPs) for stiff systems of ordinarydifferential equations (ODEs). In a recent performance evaluationof various IVP solvers, including a variable-step-variable-orderimplementation of the MEBDF method by Cash, it turned out thatthe MEBDF code often performs more efficiently than codes likeRADAU5, DASSL and VODE. This motivated us to look at possibleparallel implementations of the MEBDF method. Each MEBDF stepessentially consists of successively solving three non-linearsystems by means of modified Newton iteration using the sameJacobian matrix. In a direct implementation of the MEBDF methodon a parallel computer system, the only scope for (coarse grain)parallelism consists of a number of parallel vector updates.However, all forward–backward substitutions and all right-hand-sideevaluations have to be done in sequence. In this paper, ourstarting point is the original (unmodified) EBDF method. Asa consequence, two different Jacobian matrices are involvedin the modified Newton method, but on a parallel computer system,the effective Jacobian-evaluation and the LU decomposition costsare not increased. Furthermore, we consider the simultaneoussolution, rather than the successive solution, of the threenon-linear systems, so that in each iteration the forward–backwardsubstitutions and the right-hand-side evaluations can be doneconcurrently. A mutual comparison of the performance of theparallel EBDF approach and the MEBDF approach shows that wecan expect a speed-up factor of about 2 on three processors.     

A parareal waveform relaxation algorithm for semi-linear parabolic partial differential equations

We report a new parallel iterative algorithm for semi-linear parabolic partial differential equations (PDEs) by combining a kind of waveform relaxation (WR) techniques into the classical parareal algorithm. The parallelism can be simultaneously exploited by WR and parareal in different directions. We provide sharp error estimations for the new algorithm on bounded time domain and on unbounded time domain, respectively. The iterations of the parareal and the WR are balanced to optimize the performance of the algorithm. Furthermore, the speedup and the parallel efficiency of the new approach are analyzed. Numerical experiments are carried out to verify the effectiveness of the theoretic work.     

A family of explicit parallel Runge-Kutta-Nyström methods

A procedure for the construction of high-order explicit parallel Runge-Kutta-Nyström (RKN) methods for solving second-order nonstiff initial value problems (IVPs) is analyzed. The analysis reveals that starting the procedure with a reference symmetric RKN method it is possible to construct high-order RKN schemes which can be implemented in parallel on a small number of processors. These schemes are defined by means of a convex combination of k disjoint s_i-stage explicit RKN methods which are constructed by connecting s_i steps of a reference explicit symmetric method. Based on the reference second-order Störmer-Verlet methods we derive a family of high-order explicit parallel schemes which can be implemented in variable-step codes without additional cost. The numerical experiments carried out show that the new parallel schemes are more efficient than some sequential and parallel codes proposed in the scientific literature for solving second-order nonstiff IVPs.     

Parker-Sochacki method for solving systems of ordinary differential equations using graphics processors

The Parker-Sochacki method, which is used for solving systems of ordinary differential equations, and implementation of this method on graphics processors are described. The solution to the classical N-body problem is considered as a test. The algorithm makes it possible to effectively use massive parallel graphics processors and provides acceptable accuracy with multiple time reduction, as compared to processors of a conventional architecture.     

Simulation and solution of stochastic systems with a component-based software design

Often uncertainties occur in the numerical simulation of real-world problems. A useful method is stochastic modelling. In this context Monte Carlo (MC) methods and the stochastic Galerkin method are well known and frequently used. In the field of component-based software systems a component can be seen as an independent software or as a part of a software system. A corresponding architecture allows clearness, flexibility and reusability. The Component Template Library (CTL) is a so called middleware to realise component-based software systems, and fits to the requirements of scientific computing. In this paper we present the design of a software system for the stochastic simulation using MC and stochastic Galerkin methods based on the CTL. For each stochastic method a number of components are realised to allow parallel features. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

WAVEFORM RELAXATION METHODS AND ACCURACY INCREASE

ＷＡＶＥＦＯＲＭ　ＲＥＬＡＸＡＴＩＯＮ　ＭＥＴＨＯＤＳ　ＡＮＤ　ＡＣＣＵＲＡＣＹ　ＩＮＣＲＥＡＳＥＷＡＶＥＦＯＲＭＲＥＬＡＸＡＴＩＯＮＭＥＴＨＯＤＳＡＮＤＡＣＣＵＲＡＣＹＩＮＣＲＥＡＳＥ￥ＳｏｎｇＹｏｎｇｚｈｏｎｇ（ＮａｎｊｉｎｇＮｏｒｍａｌＵｎｉｖ...     

Parallel linear system solvers for Runge-Kutta methods

If the nonlinear systems arising in implicit Runge-Kutta methods like the Radau IIA methods are iterated by (modified) Newton, then we have to solve linear systems whose matrix of coefficients is of the form I-A hJ with A the Runge-Kutta matrix and J an approximation to the Jacobian of the righthand side function of the system of differential equations. For larger systems of differential equations, the solution of these linear systems by a direct linear solver is very costly, mainly because of the LU-decompositions. We try to reduce these costs by solving the linear systems by a second (inner) iteration process. This inner iteration process is such that each inner iteration again requires the solution of a linear system. However, the matrix of coefficients in these new linear systems is of the form I - B hJ where B is similar to a diagonal matrix with positive diagonal entries. Hence, after performing a similarity transformation, the linear systems are decoupled into s subsystems, so that the costs of the LU-decomposition are reduced to the costs of s LU-decompositions of dimension d. Since these LU-decompositions can be computed in parallel, the effective LU-costs on a parallel computer system are reduced by a factor s ³ . It will be shown that matrices B can be constructed such that the inner iterations converge whenever A and J have their eigenvalues in the positive and nonpositive halfplane, respectively. The theoretical results will be illustrated by a few numerical examples. A parallel implementation on the four-processor Cray-C98/4256 shows a speed-up ranging from at least 2.4 until at least 3.1 with respect to RADAU5 applied in one-processor mode.     

Nondeterministic dynamic programming on a parallel coprocessing system

We describe how nondeterministic dynamic programming (DP) algorithms can be designed for a new class of parallel coprocessing systems using “functional memory”, an architecture based upon dataflow computer principles. We also show how Petri nets can be used to model and express such parallel DP algorithms. Finally, we discuss architectural improvements that would facilitate the processing of Petri net models of nondeterministic DP algorithms on functional memory computers (FMC).     

Indefinitely preconditioned conjugate gradient method for large sparse equality and inequality constrained quadratic problems

This paper is concerned with the numerical solution of a symmetric indefinite system which is a generalization of the Karush–Kuhn–Tucker system. Following the recent approach of Luk?an and Vl?ek, we propose to solve this system by a preconditioned conjugate gradient (PCG) algorithm and we devise two indefinite preconditioners with good theoretical properties. In particular, for one of these preconditioners, the finite termination property of the PCG method is stated. The PCG method combined with a parallel version of these preconditioners is used as inner solver within an inexact Interior‐Point (IP) method for the solution of large and sparse quadratic programs. The numerical results obtained by a parallel code implementing the IP method on distributed memory multiprocessor systems enable us to confirm the effectiveness of the proposed approach for problems with special structure in the constraint matrix and in the objective function. Copyright © 2002 John Wiley & Sons, Ltd.     

Parallel direct methods for solving the system of linear equations with pipelining on a multicore using OpenMP

Recent developments in high performance computer architecture have a significant effect on all fields of scientific computing. Linear algebra and especially the solution of linear systems of equations lie at the heart of many applications in scientific computing. This paper describes and analyzes three parallel versions of the dense direct methods such as the Gaussian elimination method and the LU form of Gaussian elimination that are used in linear system solving on a multicore using an OpenMP interface. More specifically, we present two naive parallel algorithms based on row block and row cyclic data distribution and we put special emphasis on presenting a third parallel algorithm based on the pipeline technique. Further, we propose an implementation of the pipelining technique in OpenMP. Experimental results on a multicore CPU show that the proposed OpenMP pipeline implementation achieves good overall performance compared to the other two naive parallel methods. Finally, in this work we propose a simple, fast and reasonably analytical model to predict the performance of the direct methods with the pipelining technique.     

Efficient implementation of Jacobi's diagonalization method on the DAP

Summary The DAP architecture brings into consideration the Jacobi method where several non-interacting rotations can be performed in parallel. However the design of the algorithm is crucial in a parallel environment. In this paper we shall consider two techniques specifically designed to reduce the organisation and the arithmetic components of the parallel Jacobi method.     

Parallel Newton two‐stage methods based on ILU factorizations for nonlinear systems

Parallel iterative algorithms based on the Newton method and on two of its variants, the Shamanskii method and the Chord method, for solving nonlinear systems are proposed. These algorithms are based on two‐stage multisplitting methods where incomplete LU factorizations are considered as a mean of constructing the inner splittings. Convergence properties of these parallel methods are studied for H‐matrices. Computational results of these methods on two parallel computing systems are discussed. The reported experiments show the effectiveness of these methods. Copyright © 2006 John Wiley & Sons, Ltd.