Parallel Iterative Solution Methods for Linear Systems arising from Discretized PDE's

In these notes we will present an overview of a number of related iterative methods for the solution of linear systems of equations. These methods are so-called Krylov projection type methods and the include popular methods as Conjugate Gradients, Bi-Conjugate Gradients, CGST Bi-CGSTAB, QMR, LSQR and GMRES. We will show how these methods can be derived from simple basic iteration formulas. We will not give convergence proofs, but we will refer for these, as far as available, to litterature. Iterative methods are often used in combination with so-called preconditioning operators (approximations for the inverses of the operator of the system to be solved). Since these preconditions are not essential in the derivation of the iterative methods, we will not give much attention to them in these notes. However, in most of the actual iteration schemes, we have included them in order to facilitate the use of these schemes in actual computations. For the application of the iterative schemes one usually thinks of linear sparse systems, e.g., like those arising in the finite element or finite difference approximatious of (systems of) partial differential equations. However, the structure of the operators plays no explicit role in any of these schemes, and these schemes might also successfully be used to solve certain large dense linear systems. Depending on the situation that might be attractive in terms of numbers of floating point operations. It will turn out that all of the iterative are parallelizable in a straight forward manner. However, especially for computers with a memory hierarchy (i.e. like cache or vector registers), and for distributed memory computers, the performance can often be improved significantly through rescheduling of the operations. We will discuss parallel implementations, and occasionally we will report on experimental findings

BiCGCR2: A new extension of conjugate residual method for solving non-Hermitian linear systems

In the present paper, we introduce a new extension of the conjugate residual (CR) for solving non-Hermitian linear systems with the aim of developing an alternative basic solver to the established biconjugate gradient (BiCG), biconjugate residual (BiCR) and biconjugate A-orthogonal residual (BiCOR) methods. The proposed Krylov subspace method, referred to as the BiCGCR2 method, is based on short-term vector recurrences and is mathematically equivalent to both BiCR and BiCOR. We demonstrate by extensive numerical experiments that the proposed iterative solver has often better convergence performance than BiCG, BiCR and BiCOR. Hence, it may be exploited for the development of new variants of non-optimal Krylov subspace methods

Dielectric and Spin Susceptibilities using Density Functional Theory

The response of a system to some external perturbation is almost ubiquitous in Physics. The application of perturbation theory through an electronic structure method such as Density Functional Theory has had significant contributions over the last few decades. Its implementation, aptly named Density Functional Perturbation Theory has seen use in a number of ab initio calculations on a variety of physical properties of materials which depend on their lattice-dynamical behaviour. Specific heats, thermal expansion, infrared, Raman and optical spectra are to name just a few. Understanding the complex phenomena has significantly corroborated the current understanding of the quantum picture of solids. The Sternheimer scheme falls under the umbrella of methods to compute response functions in Time-Dependent Density Functional Theory. Initially developed to study the electronic polarisability, it is now commonly utilised in the field of lattice dynamics to study phonons and related crystal properties. The Sternheimer equation has also been used to model spin wave excitations by computation of the magnetic susceptibility. The poles of the susceptibility are known to correspond to magnon excitations and these computations have been corroborated by experimental inelastic neutron scattering data. These excitations are of a transverse nature, in that they involve fluctuations of the magnetisation perpendicular to a chosen z axis. The lesser-known longitudinal excitations involve fluctuations of the magnetisation along z, an investigation of collective modes present in transition metals may be carried out from self-consistent computations of the Sternheimer equation. The dielectric response is an important linear response function in solid-state physics. Its computation from first principles provides an invaluable tool in the characterisation of optical properties and can be compared to the experimental method of spectroscopic ellipsometry.The work in this thesis concerns the implementation of the Sternheimer method in computing the dynamical response from either an external plane wave or spin-polarised perturbation. These response functions are the dielectric and spin (magnetisation) susceptibilities respectively.The scheme to compute the frequency-dependent dielectric response is implemented in a plane-wave pseudopotential DFT package. Calculations are performed on the semiconducting systems of Silicon, Gallium Arsenide, Zinc Oxide and perovskite Methylammonium Lead Triiodide. The overall shape of the dielectric spectra is in good agreement with spectroscopic ellipsometry data, however, there is a shift which is attributed to the limitations of DFT.The scheme developed to compute longitudinal spin dynamics is applied to the transitionmetal systems of body-centred cubic Iron and face-centred cubic Nickel. In a similar manner to another first principles approach, a single dominant peak is shown to be present in the magnetisation channel with the charge dynamics being effectively null in comparison. However, the exact position of these peaks is not in agreement with the other approach, a discussion is made regarding difficulties pertaining to self-consistent optimisation

Modular vs non-modular preconditioners for fluid-structure systems with large added-mass effect

In this article we address the numerical simulation of fluid–structure interaction (FSI) problems featuring large added-mass effect. We analyze different preconditioners for the coupled system matrix obtained after space–time discretization and linearization of the FSI problem. The classical Dirichlet–Neumann preconditioner has the advantage of “modularity” because it allows to reuse existing fluid and structure codes with minimum effort (simple interface communication). Unfortunately, its performance is very poor in case of large added-mass effects. Alternatively, we consider two non-modular approaches. The first one consists in preconditioning the coupled system with a suitable diagonal scaling combined with an ILUT preconditioner. The system is then solved by a Krylov method. The drawback of this procedure is that the combination of fluid and structure codes to solve the coupled system is not straightforward. The second non-modular approach we consider is a splitting technique based on an inexact block-LU factorization of the linear FSI system. The resulting algorithm computes the fluid velocity separately from the coupled pressure–structure system at each iteration, reducing the computational cost. Independently of the preconditioner, the efficiency of semi-implicit algorithms (i.e., those that treat geometric and fluid nonlinearities in an explicit way) is highlighted and their performance compared to the one of implicit algorithms. All the methods are tested on three-dimensional blood-vessel systems. The algorithm combining the non-modular ILUT preconditioner with Krylov methods proved to be the fastest

Deflation of the Finite Pointset Method

In this thesis a deflation method for the Finite Pointset Method (FPM) is presented. FPM is a particle method based on Lagrangian coordinates to solve problems in fluid dynamics. A strong formulation of the occuring differential equations is produced by FPM, and the linear system of equations obtained by an implicit approach is solved by an iterative method such as BiCGSTAB. To improve the convergence rate of BiCGSTAB, the computational domain is divided into a number of deflation cells and a projection between the deflated domain and the original domain is constructed with the help of different ansatz functions, either constant, linear or quadratic. Also, the Moore Penrose pseudoinverse of the projection is computed. Applying the projection and restriction to the linear FPM system, a deflated system is obtained which can easily be solved with a direct method. The deflated solution is then projected onto the full domain. The deflation is tested for a number of test cases in one and two dimensions. Constant ansatz functions provide acceptable results for Dirichlet problems, but give big errors when deflating problems with mixed boundary conditions. Linear ansatz functions provide good approximations which converge to the exact solution as the number of deflation cells increases. Quadratic ansatz functions provide deflated solutions as good as the exact solutions for all test cases but are computationally expensive. The BiCGSTAB convergence rate is improved when using a deflated solution as initial guess, compared to using the zero-vector. The size of the improvement varies between which ansatz functions are used. Overall, the proposed method provides an increased convergence rate in the BiCGSTAB algorithm for FPM. However, the computational effort in the deflation process should also be taken into account

Preconditioned fast solvers for large linear systems with specific sparse and/or Toeplitz-like structures and applications

In this thesis, the design of the preconditioners we propose starts from applications instead of treating the problem in a completely general way. The reason is that not all types of linear systems can be addressed with the same tools. In this sense, the techniques for designing efficient iterative solvers depends mostly on properties inherited from the continuous problem, that has originated the discretized sequence of matrices. Classical examples are locality, isotropy in the PDE context, whose discrete counterparts are sparsity and matrices constant along the diagonals, respectively. Therefore, it is often important to take into account the properties of the originating continuous model for obtaining better performances and for providing an accurate convergence analysis. We consider linear systems that arise in the solution of both linear and nonlinear partial differential equation of both integer and fractional type. For the latter case, an introduction to both the theory and the numerical treatment is given. All the algorithms and the strategies presented in this thesis are developed having in mind their parallel implementation. In particular, we consider the processor-co-processor framework, in which the main part of the computation is performed on a Graphics Processing Unit (GPU) accelerator. In Part I we introduce our proposal for sparse approximate inverse preconditioners for either the solution of time-dependent Partial Differential Equations (PDEs), Chapter 3, and Fractional Differential Equations (FDEs), containing both classical and fractional terms, Chapter 5. More precisely, we propose a new technique for updating preconditioners for dealing with sequences of linear systems for PDEs and FDEs, that can be used also to compute matrix functions of large matrices via quadrature formula in Chapter 4 and for optimal control of FDEs in Chapter 6. At last, in Part II, we consider structured preconditioners for quasi-Toeplitz systems. The focus is towards the numerical treatment of discretized convection-diffusion equations in Chapter 7 and on the solution of FDEs with linear multistep formula in boundary value form in Chapter 8

Lanczos-type solvers for nonsymmetric linear systems of equations

Among the iterative methods for solving large linear systems with a sparse (or, possibly, structured) nonsymmetric matrix, those that are based on the Lanczos process feature short recurrences for the generation of the Krylov space. This means low cost and low memory requirement. This review article introduces the reader not only to the basic forms of the Lanczos process and some of the related theory, but also describes in detail a number of solvers that are based on it, including those that are considered to be the most efficient ones. Possible breakdowns of the algorithms and ways to cure them by look-ahead are also discusse