Efficient implementation of Gauss collocation and Hamiltonian boundary value methods

In this paper we define an efficient implementation for the family of low-rank energy-conserving Runge-Kutta methods named Hamiltonian Boundary Value Methods (HBVMs), recently defined in the last years. The proposed implementation relies on the particular structure of the Butcher matrix defining such methods, for which we can derive an efficient splitting procedure. The very same procedure turns out to be automatically suited for the efficient implementation of Gauss-Legendre collocation methods, since these methods are a special instance of HBVMs. The linear convergence analysis of the splitting procedure exhibits excellent properties, which are confirmed by a few numerical tests.     

Efficient implementation of Radau collocation methods

In this paper we define an efficient implementation of Runge–Kutta methods of Radau IIA type, which are commonly used when solving stiff ODE-IVPs problems. The proposed implementation relies on an alternative low-rank formulation of the methods, for which a splitting procedure is easily defined. The linear convergence analysis of this splitting procedure exhibits excellent properties, which are confirmed by its performance on a few numerical tests.     

An efficient time-step-based self-adaptive algorithm for predictor-corrector methods of Runge-Kutta type

Finding an efficient implementation variant for the numerical solution of problems from computational science and engineering involves many implementation decisions that are strongly influenced by the specific hardware architecture. The complexity of these architectures makes it difficult to find the best implementation variant by manual tuning. For numerical solution methods from linear algebra, auto-tuning techniques based on a global search engine as they are used for ATLAS or FFTW can be used successfully. These techniques generate different implementation variants at installation time and select one of these implementation variants either at installation time or at runtime, before the computation starts. For some numerical methods, auto-tuning at installation time cannot be applied directly, since the best implementation variant may strongly depend on the specific numerical problem to be solved. An example is solution methods for initial value problems (IVPs) of ordinary differential equations (ODEs), where the coupling structure of the ODE system to be solved has a large influence on the efficient use of the memory hierarchy of the hardware architecture. In this context, it is important to use auto-tuning techniques at runtime, which is possible because of the time-stepping nature of ODE solvers.In this article, we present a sequential self-adaptive ODE solver that selects the best implementation variant from a candidate pool at runtime during the first time steps, i.e., the auto-tuning phase already contributes to the progress of the computation. The implementation variants differ in the loop structure and the data structures used to realize the numerical algorithm, a predictor-corrector (PC) iteration scheme with Runge-Kutta (RK) corrector considered here as an example. For those implementation variants in the candidate pool that use loop tiling to exploit the memory hierarchy of a given hardware platform we investigate the selection of tile sizes. The self-adaptive ODE solver combines empirical search with a model-based approach in order to reduce the search space of possible tile sizes. Runtime experiments demonstrate the efficiency of the self-adaptive solver for different IVPs across a range of problem sizes and on different hardware architectures.     

On systolic arrays for updating the Cholesky factorization

We have considered the systolic implementation of several methods for updating the Cholesky factorization. For positive rank-k changes there are simple one-pass arrays that implement algorithms based on elimination and plane rotations. In the case of negative rank-one changes, we do not feel that the standard algorithm [2] has a practical implementation. We have introduced a new algorithm for the case of a negative rank-k change and provided an attractive two-pass systolic implementation.     

An efficient time-step-based self-adaptive algorithm for predictor–corrector methods of Runge–Kutta type

Finding an efficient implementation variant for the numerical solution of problems from computational science and engineering involves many implementation decisions that are strongly influenced by the specific hardware architecture. The complexity of these architectures makes it difficult to find the best implementation variant by manual tuning. For numerical solution methods from linear algebra, auto-tuning techniques based on a global search engine as they are used for ATLAS or FFTW can be used successfully. These techniques generate different implementation variants at installation time and select one of these implementation variants either at installation time or at runtime, before the computation starts. For some numerical methods, auto-tuning at installation time cannot be applied directly, since the best implementation variant may strongly depend on the specific numerical problem to be solved. An example is solution methods for initial value problems (IVPs) of ordinary differential equations (ODEs), where the coupling structure of the ODE system to be solved has a large influence on the efficient use of the memory hierarchy of the hardware architecture. In this context, it is important to use auto-tuning techniques at runtime, which is possible because of the time-stepping nature of ODE solvers.In this article, we present a sequential self-adaptive ODE solver that selects the best implementation variant from a candidate pool at runtime during the first time steps, i.e., the auto-tuning phase already contributes to the progress of the computation. The implementation variants differ in the loop structure and the data structures used to realize the numerical algorithm, a predictor–corrector (PC) iteration scheme with Runge–Kutta (RK) corrector considered here as an example. For those implementation variants in the candidate pool that use loop tiling to exploit the memory hierarchy of a given hardware platform we investigate the selection of tile sizes. The self-adaptive ODE solver combines empirical search with a model-based approach in order to reduce the search space of possible tile sizes. Runtime experiments demonstrate the efficiency of the self-adaptive solver for different IVPs across a range of problem sizes and on different hardware architectures.     

On implementing a primal-dual interior-point method for conic quadratic optimization

 Based on the work of the Nesterov and Todd on self-scaled cones an implementation of a primal-dual interior-point method for solving large-scale sparse conic quadratic optimization problems is presented. The main features of the implementation are it is based on a homogeneous and self-dual model, it handles rotated quadratic cones directly, it employs a Mehrotra type predictor-corrector extension and sparse linear algebra to improve the computational efficiency. Finally, the implementation exploits fixed variables which naturally occurs in many conic quadratic optimization problems. This is a novel feature for our implementation. Computational results are also presented to document that the implementation can solve very large problems robustly and efficiently. Received: November 18, 2000 / Accepted: January 18, 2001 Published online: September 27, 2002 Key Words. conic optimization – interior-point methods – large-scale implementation     

A domain-specific embedded language in C++ for lowest-order discretizations of diffusive problems on general meshes

In this work we propose an original implementation of a large family of lowest-order methods for diffusive problems including standard and hybrid finite volume methods, mimetic finite difference-type schemes, and cell centered Galerkin methods. The key idea is to regard the method at hand as a (Petrov–)Galerkin scheme based on possibly incomplete, broken affine spaces defined from a gradient reconstruction and a point value. The resulting unified framework serves as a basis for the development of a FreeFEM-like domain specific language targeted at defining discrete linear and bilinear forms. Both the back-end and the front-end of the language are extensively discussed, and several examples of applications are provided. The overhead of the language is evaluated with respect to a more traditional implementation. A benchmark including the comparison with more classical finite element methods on standard meshes is also proposed.     

Tensor-Krylov methods for large nonlinear equations

In this paper, we describe tensor methods for large systems of nonlinear equations based on Krylov subspace techniques for approximately solving the linear systems that are required in each tensor iteration. We refer to a method in this class as a tensor-Krylov algorithm. We describe comparative testing for a tensor-Krylov implementation versus an analogous implementation based on a Newton-Krylov method. The test results show that tensor-Krylov methods are much more efficient and robust than Newton-Krylov methods on hard nonlinear equations problems.Part of this work was performed while the author was research associate at CERFACS (Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique).Research supported in part by the Office of Scientific Computing, U.S. Department of Energy, under Contract W-31-109-Eng-38.     

Symplectic phase flow approximation for the numerical integration of canonical systems

Summary New methods are presented for the numerical integration of ordinary differential equations of the important family of Hamiltonian dynamical systems. These methods preserve the Poincaré invariants and, therefore, mimic relevant qualitative properties of the exact solutions. The methods are based on a Runge-Kutta-type ansatz for the generating function to realize the integration steps by canonical transformations. A fourth-order method is given and its implementation is discussed. Numerical results are presented for the Hénon-Heiles system, which describes the motion of a star in an axisymmetric galaxy.     

Solution of Differential--Algebraic Systems Using Diagonally Implicit Runge--Kutta Methods

Diagonally Implicit Runge—Kutta (DIRK) methods are developedand applied to differential—algebraic systems arisingfrom dynamic process simulation. In particular, an embeddedfamily of DIRK methods is developed for implementation as avariable-step variable-order algorithm. The methods developedallow easy assessment of local solution error as well as theability to change the order of approximation. The stabilityproperties of the methods are chosen to make them suitable foruse on stiff systems. Some important aspects of implementation of DIRK methods arediscussed within the context of the solution of differential—algebraicsystems. The performance of this algorithm is compared withan alternative variable-order approach based on "triples" whichallows the patching together of several fixed-order formulae.The results indicate that the fully embedded DIRK algorithmis generally more efficient than the algorithm based on "triples".Areas of further investigation in the context of differential—algebraicsystems are outlined.     

Block-Arnoldi and Davidson methods for unsymmetric large eigenvalue problems

Summary We present two methods for computing the leading eigenpairs of large sparse unsymmetric matrices. Namely the block-Arnoldi method and an adaptation of the Davidson method to unsymmetric matrices. We give some theoretical results concerning the convergence and discuss implementation aspects of the two methods. Finally some results of numerical tests on a variety of matrices, in which we compare these two methods are reported.     

Minimization methods for approximating tensors and their comparison

Application of various minimization methods to trilinear approximation of tensors is considered. These methods are compared based on numerical calculations. For the Gauss-Newton method, an efficient implementation is proposed, and the local rate of convergence is estimated for the case of completely symmetric tensors.     

Exact and heuristic solutions of the global supply chain problem with transfer pricing

We examine the example of a multinational corporation that attempts to maximize its global after tax profits by determining the flow of goods, the transfer prices, and the transportation cost allocation between each of its subsidiaries. Vidal and Goetschalckx [Vidal, C.J., Goetschalckx, M., 2001. A global supply chain model with transfer pricing and transportation cost allocation. European Journal of Operational Research 129 (1), 134–158] proposed a bilinear model of this problem and solved it by an Alternate heuristic. We propose a reformulation of this model reducing the number of bilinear terms and accelerating considerably the exact solution. We also present three other solution methods: an implementation of Variable Neighborhood Search (VNS) designed for any bilinear model, an implementation of VNS specifically designed for the problem considered here and an exact method based on a branch and cut algorithm. The solution methods are tested on artificial instances. These results show that our implementation of VNS outperforms the two other heuristics. The exact method found the optimal solution of all small instances and of 26% of medium instances.     

Effective cell-centred time-domain Maxwell’s equations numerical solvers

This research work analyses techniques for implementing a cell-centred finite-volume time-domain (ccFV-TD) computational methodology for the purpose of studying microwave heating. Various state-of-the-art spatial and temporal discretisation methods employed to solve Maxwell’s equations on multi-dimensional structured grid networks are investigated, and the dispersive and dissipative errors inherent in those techniques examined. Both staggered and unstaggered grid approaches are considered. Upwind schemes using a Riemann solver and intensity vector splitting are studied and evaluated. Staggered and unstaggered Leapfrog and Runge–Kutta time integration methods are analysed in terms of phase and amplitude error to identify which method is the most accurate and efficient for simulating microwave heating processes. The implementation and migration of typical electromagnetic boundary conditions from staggered in space to cell-centred approaches also is deliberated. In particular, an existing perfectly matched layer absorbing boundary methodology is adapted to formulate a new cell-centred boundary implementation for the ccFV-TD solvers. Finally for microwave heating purposes, a comparison of analytical and numerical results for standard case studies in rectangular waveguides allows the accuracy of the developed methods to be assessed.     

Automatic step size and order control in implicit one-step extrapolation methods

A theory is presented for implicit one-step extrapolation methods for ordinary differential equations. The computational schemes used in such methods are based on the implicit Runge-Kutta methods. An efficient implementation of implicit extrapolation is based on the combined step size and order control. The emphasis is placed on calculating and controlling the global error of the numerical solution. The aim is to achieve the user-prescribed accuracy in an automatic mode (ignoring round-off errors). All the theoretical conclusions of this paper are supported by the numerical results obtained for test problems.     

Parallel iterated methods based on multistep Runge-Kutta methods of Radau type

This paper investigates iterated Multistep Runge-Kutta methods of Radau type as a class of explicit methods suitable for parallel implementation. Using the idea of van der Houwen and Sommeijer [18], the method is designed in such a way that the right-hand side evaluations can be computed in parallel. We use stepsize control and variable order based on iterated approximation of the solution. A code is developed and its performance is compared with codes based on iterated Runge-Kutta methods of Gauss type and various Dormand and Prince pairs [15]. The accuracy of some of our methods are comparable with the PIRK10 methods of van der Houwen and Sommeijer [18], but require fewer processors. In addition at very stringent tolerances these new methods are competitive with RK78 pairs in a sequential implementation.     

Introducing <Emphasis FontCategory="NonProportional" Type="Bold">libeemd</Emphasis>: a program package for performing the ensemble empirical mode decomposition

The ensemble empirical mode decomposition (EEMD) and its complete variant (CEEMDAN) are adaptive, noise-assisted data analysis methods that improve on the ordinary empirical mode decomposition (EMD). All these methods decompose possibly nonlinear and/or nonstationary time series data into a finite amount of components separated by instantaneous frequencies. This decomposition provides a powerful method to look into the different processes behind a given time series data, and provides a way to separate short time-scale events from a general trend. We present a free software implementation of EMD, EEMD and CEEMDAN and give an overview of the EMD methodology and the algorithms used in the decomposition. We release our implementation, libeemd, with the aim of providing a user-friendly, fast, stable, well-documented and easily extensible EEMD library for anyone interested in using (E)EMD in the analysis of time series data. While written in C for numerical efficiency, our implementation includes interfaces to the Python and R languages, and interfaces to other languages are straightforward.     

A Parallel Preconditioned Iterative Realization of the Panel Method in 3D

The parallel version of precondition iterative techniques is developed for matrices arising from the panel boundary element method for three-dimensional simple connected domains with Dirichlet boundary conditions. Results were obtained on an nCube-2 parallel computer showing that preconditioned iterative methods are very well suited also in three-dimensional cases for implementation on an MIMD computer and that they are much more efficient than usual direct solution techniques.     

Computational results of an interior point algorithm for large scale linear programming

This paper gives computational results for an efficient implementation of a variant of dual projective algorithm for linear programming. The implementation uses the preconditioned conjugate gradient method for computing projections. Our computational experience reported in this paper indicates that this algorithm has potential as an alternative for solving very large LPs in which the direct methods fail due to memory and CPU time requirements. The conjugate gradient algorithm was able to find very accurate directions even when the system was ill-conditioned. The paper also discusses a new mathematical technique called the reciprocal estimates for estimating the primal variables. We have conducted extensive computational experiments on problems representative of large classes of applications of current interest. We have also chosen instances of the problems of future potential interest, which could not be solved in the past due to the weakness of the prior solution methods, but which represent a large class of new applications. The hypergraph model is such an example. Comparison of our implementation with MINOS 5.1 shows that our implementation is orders of magnitude faster than MINOS 5.1 for these problems.     

Singly implicit diagonally extended Runge-Kutta methods of fourth order

Singly implicit diagonally extended Runge-Kutta methods make it possible to combine the merits of diagonally implicit methods (namely, the simplicity of implementation) and fully implicit ones (high stage order). Due to this combination, they can be very efficient at solving stiff and differential-algebraic problems. In this paper, fourth-order methods with an explicit first stage are examined. The methods have the third or fourth stage order. Consideration is given to an efficient implementation of these methods. The results of tests in which the proposed methods were compared with the fifth-order RADAU IIA method are presented.