A parallel electrostatic solver for the VORPAL code

We are interested in the design of efficient algebraic multigrid (AMG) methods for the solution of large sparse systems of linear equations arising from finite element (FE) discretization of second-order elliptic partial differential equations (PDEs). In particular, we introduce the concept of so-called “edge matrices”, which–in the present context–are extracted from the individual element matrices. This allows for the construction of spectrally equivalent approximations of the original stiffness matrix that can be utilized in the framework of AMG. The edge matrices give rise to modify the definition of “strong” and “weak” connections (edges), which provides a basis for selecting the coarse-grid nodes in algebraic multigrid methods. Moreover, a reproduction of edge matrices on coarse levels offers the opportunity to combine classical coarsening algorithms with effective (energy minimizing) interpolation principles involving small-sized “computational molecules” (small collections of edge matrices). This yields a flexible and robust new variant of AMG, which we refer to as AMGm.

The polarizable continuum model (PCM) is a well established method for computing solvation effects. Its attractiveness comes from the fact that the solution is modelled as a continuum instead of modelling each atom individually. This allows for more complex simulations. Nevertheless, there are still a couple of challenges to face. One of them being the computational time required at the limit of large systems, or when high accuracy is needed. The wavelet boundary element method can be used to overcome these problems, provided a reliable cavity generator is used. The PCM calculates the molecular free energy in solution as the sum of electrostatic, dispersion-repulsion and the cavitation energy by solving the underlying partial differential equations. Those differential equations can be transformed into integral equations, solely defined on the boundary of the molecular cavity, by applying the integral equation formalism (IEFPCM). The integral equations can then be discretized using the wavelet boundary element method. The resulting sparse system of linear equations can be solved reliably by iterative solvers, leading to high accuracy solutions even for large systems. The challenge of this approach lies in the generation and refinement of the molecular cavities on which the interactions take place. Since wavelets are defined as a tensor product of one-dimensional functions, the wavelet boundary element method needs a discretization into quadrangular patches of the molecular cavity. This feature is missing in common commercial mesh generation tools, which focus on generating triangular meshes. Making use of characteristic functions for defining which points are inside the cavity, a model for generating the molecular surface independently of its exact description is achieved. We present here a way for generating triangular meshes which are subsequently merged into quadrangular ones. We apply geometrical quantities that measure the fitness of the shapes involved and improve the position of individual points iteratively. Thereby, we generate quadrangular meshes with an unprecedented quality, which ultimately can be refined in a hierarchical manner for the use of wavelets. Overall, this method can generate reliable parametrizations on a variety of smooth cavities. With the help of the wavelet boundary element method, the PCM equations can thus be solved with high accuracy also on large systems.

Many geo-scientific applications involve boundary value problems arising in simulating electrostatic and electromagnetic fields for geophysical prospecting and subsurface imaging of electrical resistivity. Modeling complex geological media with three-dimensional finite-difference grids gives rise to large sparse linear systems of equations. For such systems, we have implemented three common iterative Krylov solution methods on graphics processing units and compared their performance with parallel host-based versions. The benchmarks show that the device efficiency improves with increasing grid sizes. Limitations are currently given by the device memory resources.

Performance assessment and field-development planning for carbon geo-sequestration relies on reservoir simulation. The standard approach is to use regular (corner-point) grids in space and a fully implicit discretisation in time, solving for all primary unknowns at once using Newton’s method. While the former limits the physical realism of the simulation model built from the geomodel, the latter leads to a large ill-conditioned system of equations that is inefficient to solve. We overcome the regularisation issue via a hierarchy of fully unstructured, geobody-conforming, finite element meshes which can be refined until mesh convergence is achieved. To address the time discretisation issue we employ, for the first time, the linearly implicit extrapolation scheme (LIMEX) to solve the highly non-linear coupled two-phase flow and reactive transport equations. To solve the arising large sparse systems of linear equations, we apply the geometric multigrid (GMG) method that demonstrates optimal, linear complexity and allows an efficient parallelization on supercomputers. Another novel feature of our formulation is the consideration of the kinetics of mass transfer between the carbonic and aqueous phases. This approach removes the first-order dependence on mesh refinement of CO2 that is dissolved, a characteristic feature of standard equilibrium models. We demonstrate the parallel scalability of our simulation framework with an implementation based on the UG4 platform. Proof-of-concept results accurately capture key features of CO2 migration including filtration by capillary barriers and convective dissolution of CO2 at the base of the plume.

Computer subroutines for solving large sparse systems of linear equations usually involve fairly complicated data structures and storage management. In most cases, the user of such subroutines simply wants top solve his problem and does not have to understand how the storage management is done, or how the matrix components are stored. Thus software for solving sparse systems of linear equations should be designed in such a way that it can insulate the user from these considerations, while still allowing the software to be used in a variety of ways.

Many applications in science and engineering give rise to large sparse linear systems of equations that need to be solved as efficiently as possible. As the size of the problems of interest increases, it can become necessary to consider exploiting multiprocessors to solve these systems. We report on the design and development of parallel frontal solvers for the numerical solution of large sparse linear systems. Three codes have been developed for the mathematical software library HSL (www.cse.clrc.ac.uk/Activity/HSL). The first is for unsymmetric finite-element problems; the second is for symmetric positive definite finite-element problems; and the third is for highly unsymmetric linear systems such as those that arise in chemical process engineering. In each case, the problem is subdivided into a small number of loosely connected subproblems and a frontal method is then applied to each of the subproblems in parallel. We discuss how our software is designed to achieve the goals of portability, ease of use, efficiency, and flexibility, and illustrate the performance using problems arising from real applications.

Solving a large sparse system of linear equations (SSLE)over the binary field is an important step in the general number field seive (GNFS) algorithm that is often used to factorize large integers used in the RSA algorithm. Block Lanczos method proposed by Montgomery is an efficient and reliable method for solving such systems. A number of implementations of Block Lanczos method on distributed systems have been proposed in past. This paper discusses the parallel Montgomery’s Block Lanczos method for solving binary system of equations on shared multiprocessors(SMP). A simple experiment shows that the speed of convergence of this algorithm is dependent on a wise choice of initial guess for kick-starting the algorithm. A somewhat dense choice of the initial guess is suggested rather than a random choice for faster convergence. Devoid of any communication overheads of a loosely coupled system, the improved method gives good performance on SMP systems which can provide an able alternative to the more popular distributed systems. Details of implementation of this algorithm on SMP and experimental results are also provided in the paper.

This paper presents a novel hardware solution for solving sparse matrices using emulation technology. The proposed solution introduces a bio-inspired technique, which is Genetic Algorithm (GA) to solve the sparse system of linear equations (SLEs). The proposed solution can be used with any physical system that can be modeled as SLEs. Computations involved in Finite Element Method (FEM) consume too much time, which affects the final time-to-market value. Profiling shows that the most time-consuming part in the simulation process is the solver part which is responsible for solving the resultant SLEs generated from the FEM. The total number of equations may reach thousands or millions of linear equations. Hardware implementation of GA and making a good use of parallelism in implementing the architecture accelerate the time taken by the solver part. Solving the system by classical numerical methods, such as conjugate gradient (CG) implemented in software is performed sequentially, here the architecture allows the parallelism in performing GA operators, such as selection, crossover, and mutation. The parallelism in the proposed architecture is the superiority point which we rely on. The results show that our proposed method to solve the SLEs outperforms the classical CG software implementation on run time and number of equations.

Element-by-element iterative solution methods have been shown to be effective for solving large, sparse linear systems of equations with a positive definite stiffness matrix. In this paper the iterative method GMRES will be described. This solution method is also suited for solving nonsymmetric systems of equations. As acceleration technique an element-by-element preconditioner has been employed. Attention is focused mainly on the convergence behaviour of the method. The calculations described in this paper have been performed with the finite element package DIANA.

We present one of the first attempt at implementing a massively parallel frequency-domain full-waveform inversion algorithm for imaging 3D acoustic media. The inverse method is based on a classic steepest-descent algorithm. The algorithm was designed so that one or several frequencies are inverted at a time. Wave propagation modeling, a key component of the inversion algorithm, is performed with a finite-difference frequency-domain method which requires 4 grid points per wavelength. Frequency-domain methods for wave propagation modeling requires to solve a huge sparse system of linear equations whose solutions are the computed wavefields and right-hand side (RHS) terms are the sources. We solve this system with a massively parallel direct solver in order to compute efficiently multiple-source solutions. The massively parallel direct solver leaves the multiple solutions distributed over the processors. We take advantage of this distributed storage to compute in parallel the gradient of the cost function in the inversion. The algorithm was validated with preliminary synthetic examples of limited dimensions.

During the past few years, algorithmic improvements alone have reduced the time required for the direct solution of unsymmetric sparse systems of linear equations by almost an order of magnitude. This paper compares the performance of some well-known software packages for solving general sparse systems. In particular, it demonstrates the consistently high level of performance achieved by WSMP---the most recent of such solvers. It compares the various algorithmic components of these solvers and discusses their impact on solver performance. Our experiments show that the algorithmic choices made in WSMP enable it to run more than twice as fast as the best among similar solvers and that WSMP can factor some of the largest sparse matrices available from real applications in only a few seconds on a 4-CPU workstation. Thus, the combination of advances in hardware and algorithms makes it possible to solve those general sparse linear systems quickly and easily that might have been considered too large until recently.

BackgroundThe aim of electroencephalogram (EEG) source localization is to find the brain areas responsible for EEG waves of interest. It consists of solving forward and inverse problems. The forward problem is solved by starting from a given electrical source and calculating the potentials at the electrodes. These evaluations are necessary to solve the inverse problem which is defined as finding brain sources which are responsible for the measured potentials at the EEG electrodes.MethodsWhile other reviews give an extensive summary of the both forward and inverse problem, this review article focuses on different aspects of solving the forward problem and it is intended for newcomers in this research field.ResultsIt starts with focusing on the generators of the EEG: the post-synaptic potentials in the apical dendrites of pyramidal neurons. These cells generate an extracellular current which can be modeled by Poisson's differential equation, and Neumann and Dirichlet boundary conditions. The compartments in which these currents flow can be anisotropic (e.g. skull and white matter). In a three-shell spherical head model an analytical expression exists to solve the forward problem. During the last two decades researchers have tried to solve Poisson's equation in a realistically shaped head model obtained from 3D medical images, which requires numerical methods. The following methods are compared with each other: the boundary element method (BEM), the finite element method (FEM) and the finite difference method (FDM). In the last two methods anisotropic conducting compartments can conveniently be introduced. Then the focus will be set on the use of reciprocity in EEG source localization. It is introduced to speed up the forward calculations which are here performed for each electrode position rather than for each dipole position. Solving Poisson's equation utilizing FEM and FDM corresponds to solving a large sparse linear system. Iterative methods are required to solve these sparse linear systems. The following iterative methods are discussed: successive over-relaxation, conjugate gradients method and algebraic multigrid method.ConclusionSolving the forward problem has been well documented in the past decades. In the past simplified spherical head models are used, whereas nowadays a combination of imaging modalities are used to accurately describe the geometry of the head model. Efforts have been done on realistically describing the shape of the head model, as well as the heterogenity of the tissue types and realistically determining the conductivity. However, the determination and validation of the in vivo conductivity values is still an important topic in this field. In addition, more studies have to be done on the influence of all the parameters of the head model and of the numerical techniques on the solution of the forward problem.

A projection method for solving nonsymmetric linear systems on multiprocessors

P aS tiX: a high-performance parallel direct solver for sparse symmetric positive definite systems

Numerical study on incomplete orthogonal factorization preconditioners