In this chapter we present discretizations for mixed initial value/boundary value problems (IVP/BVP) and for relaxation iterative solvers associated with such discretizations. The analogy between iterative procedures and equations of evolution, especially of parabolic type (diffusion), was realized about two centuries ago, but a rigorous connection was not established until the mid-1950s.

In the following, we first consider various mixed discretizations, and subsequently we derive some of the most popular iterative solvers. Our emphasis will be on parallel computing: A good algorithm is not simply the one that converges faster but also the one that is parallelizable. The Jacobi algorithm is such an example; forgotten for years in favor of the Gauss-Seidel algorithm, which converges twice as fast for about the same computational work, it was rediscovered during the past two decades as it is trivially parallelizable, and today it is used mostly as a preconditioner for multigrid methods. The Gauss-Seidel algorithm, although faster on a serial computer, is not parallelizable unless a special multicolor algorithm is employed, as we explain in Section 7.2.4. Based on these two basic algorithms, we present the multigrid method that exploits their good convergence properties but in a smart adaptive way.

On the parallel computing side, we introduce three new commands: MPI_Gather, MPI_Allgather, and MPI_Scatter. Both MPI_Gather and MPI_Allgather are used for gathering information from a collection of processes. MPI_Scatter is used to scatter data from one process to a collection of processes. In addition to providing syntax and usage information, we present the applicability of the gathering functions in the parallel implementation of the Jacobi method.

COMPILATION GUIDE

For the purposes of this compilation example, we will assume that we are using the GNU g++ compiler to compile a C++ program we have written contained within the file myprog.cpp. In the following examples, the argument following the “-o” flag designates the file name to be used for the output. If no “-o” option is specified, most compilers default to using the name “a.out.” We now present several different programming scenarios:

1. • No user-defined libraries or user-defined header files are needed, and no special system libraries (such as those associated with math.h) are needed):

2. g++ -o myprog myprog.cpp

3. • No user-defined libraries or user-defined header files are needed, but the special system library corresponding to math.h is needed:

4. g++ -o myprog myprog.cpp -lmath

5. • User-defined libraries, user-defined header files, and the special system library corresponding to math.h are needed:

6. g++ -o myprog myprog.cpp -I/users/kirby/includes -L/users/ kirby/libs -lSCmathlib -lmath

The string following the “-I” flag designates the location of the user-defined header files to be included. The string following the “-L” flag designates the location of the user-defined libraries to be included. The string “-lSCmathlib” links the program with the user-defined library we created, and the string “-lmath” links the program with the system math library corresponding to math.h.

In this chapter we consider implicit discretizations of space and time derivatives. Unlike the explicit discretizations presented in the previous chapter, here we express a derivative at one grid point in terms of function values as well as derivative values at adjacent grid points (spatial discretization) or in terms of previous and current time levels (temporal discretization). This, in turn, implies that there is implicit coupling, and thus matrix inversion is required to obtain the solution.

The material of this chapter serves to introduce solutions of tridiagonal systems and correspondingly parallel computing of sparse linear systems using MPI. We also introduce two new MPI functions: MPI_Barrier, used for synchronizing processes, and MPI_Wtime, used for obtaining the wall-clock timing information.

IMPLICIT SPACE DISCRETIZATIONS

The discretizations we present here are appropriate for any order spatial derivative involved in a partial differential equation, but they are particularly useful when high-order accuracy and locality of data is sought. The explicit finite differences could also lead to high accuracy but at the expense of long stencils, and this, in turn, implies coupling involving many grid points and consequently a substantial communications overhead. In contrast, the implicit finite differences employ very compact stencils and guarantee locality, which is the key to success of any parallel implementation. We only consider discretizations on uniform (i.e., equidistant) grids, as we assume that a mapping of the form presented in the previous chapter is always available to transform a nonuniform grid to a uniform one. We also present discretizations only for one-dimensional grids since multidimensional discretizations are accomplished using directional splitting, as before.

Two of the most common tasks in scientific computing are interpolation of discrete data and approximation by known functions of the numerical solution, the source terms, and the boundary or initial conditions. Therefore, we need to perform these tasks both accurately and efficiently. The data are not always nicely distributed on a uniform lattice or grid, and thus we must learn how to manage these situations as well. We often use polynomials to represent discrete data because they are easy to “manipulate,” that is, differentiate and integrate. However, sines and cosines as well as special functions called wavelets are very effective means to perform interpolation and approximation, and they have very interesting properties.

In this section, we will study various such representations and their corresponding C++ implementations. We consider cases where the data are just sufficient to determine exactly the representation (deterministic case) as well as cases where the data are more than the information needed (overdetermined case).

Finally, we will present a more detailed discussion of MPI_Send and MPI_Recv, the two fundamental building blocks of MPI.

In this chapter we introduce the main themes that we will cover in this book and provide an introduction for each of them. We begin with a brief overview of C++ and define the two basic concepts of functions and classes as well as other syntactic elements of the language. We then introduce basic mathematical concepts that include elements of linear algebra, vector orthogonalization, and corresponding codes and software. Finally, we introduce parallel programming and review some generic parallel architectures as well as standard parallel algorithms for basic operations (e.g., the fan-in algorithm for recursive doubling). We also provide a brief overview of the main MPI commands.

INTRODUCTION TO C++

An ancient proverb states that the beginning of a thousand-mile journey begins with a single step. For us, this single step will be a brief overview of the C++ programming language. This introduction is not designed to be all-inclusive, but rather it should provide the scaffolding from which we will build concepts throughout this book. Admittedly, what you will read now may seem daunting in its scope, but as you become more familiar with the concepts found herein, you will be able to use the C++ language as a tool for furthering your understanding of deeper mathematical and algorithmic concepts presented later in the book. With this in mind, let us begin our thousand-mile journey with this first step.

In this chapter we consider explicit discretizations of space and time derivatives. In such discretizations we can express directly a derivative at one grid point in terms of function values at adjacent grid points (spatial discretizations) or in terms of previous time levels (temporal discretizations). This, in turn, implies that there is no implicit coupling, and thus no matrix inversion is involved; instead, only simple daxpy type operations are required.

The material in this chapter is relatively easy to program both on serial as well as on parallel computers. It is appropriate for demonstrating fundamental concepts of discretization and primary constructs of the C++ language and of the MPI library. Specifically, we will demonstrate the use of loops, arrays, functions, and passing functions to functions. In addition to presenting MPI_Send and MPI_Recv implementations for finite differences, we also introduce MPI_Sendrecv and MPI_Sendrecv_replace as alternative advanced MPI function calls for parallelizing finite differences discretizations.

EXPLICIT SPACE DISCRETIZATIONS

Basics

The formulation of derivatives based on function values on a set of points, which we call the grid, dates back to Euler in the beginning of the eighteenth century. However, advances have of course been made since then. In this section, we will formulate ways to compute first- and higher order derivatives of a function using discrete data points. The key idea is to use Taylor expansions at a subset of adjacent points of the grid, as shown in Figure 5.1.

WHAT IS SIMULATION?

Science and engineering have undergone a major transformation at the research level as well as at the development and technology level. The modern scientist and engineer spend more and more time in front of a laptop, a workstation, or a parallel supercomputer and less and less time in the physical laboratory or in the workshop. The virtual wind tunnel and the virtual biology laboratory are not a thing of the future; they are here! The old approach of “cut and try” has been replaced by “simulate and analyze” in several key technological areas such as aerospace applications, synthesis of new materials, design of new drugs, and chip processing and microfabrication. The new discipline of nanotechnology will be based primarily on large-scale computations and numerical experiments. The methods of scientific analysis and engineering design are changing continuously, affecting both our approach to the phenomena that we study as well as the range of applications that we address. Whereas there is an abundance of software available to be used as almost a “black box,” working in new application areas requires good knowledge of fundamentals and mastering of effective new tools.

In the classical scientific approach, the physical system is first simplified and set in a form that suggests what type of phenomena and processes may be important and, correspondingly, what experiments are to be conducted. In the absence of any known type of governing equations, dimensional inter dependence between physical parameters can guide laboratory experiments in identifying key parametric studies.

We have already discussed how to solve tridiagonal linear systems of equations using direct solvers (the Thomas algorithm) in Chapter 6 and some iterative solvers (Jacobi, Gauss-Seidel, SOR, and multigrid) in Chapter 7. We have also discussed solutions of nonlinear and linear systems and have introduced the conjugate gradient method in Chapter 4. In the current chapter we revisit this subject and present general algorithms for the direct and iterative solution of large linear systems. We start with the classical Gaussian elimination (which is a fast solver) and then proceed with more sophisticated solvers and preconditioners for symmetric and nonsymmetric systems.

In parallel computing, we introduce the broadcasting command MPI_Bcast and demonstrate its usefulness in the context of Gaussian elimination. In addition, we reiterate the use of MPI_Send, MPI_Recv, MPI_Allgather, and MPI_Allreduce through example implementations of algorithms presented in this chapter.

GAUSSIAN ELIMINATION

Gaussian elimination is one of the most effective ways to solve the linear system

Ax = b.

The Thomas algorithm (see Section 6.1.4) is a special case of Gaussian elimination for tridiagonal systems.

Scientific Computing is by its very nature a practical subject – it requires tools and a lot of practice. To solve realistic problems we need not only fast algorithms but also a combination of good tools and fast computers. This is the subject of the current book, which emphasizes equally all three: algorithms, tools, and computers. Oftentimes such concepts and tools are taught serially across different courses and different textbooks, and hence the interconnections among them are not immediately apparent. We believe that such a close integration is important from the outset.

The book starts with a heavy dosage of C++ and basic mathematical and computational concepts, and it ends emphasizing advanced parallel algorithms that are used in modern simulations. We have tried to make this book fun to read, to somewhat demystify the subject, and thus the style is sometimes informal and personal. It may seem that this happens at the expense of rigor, and indeed we have tried to limit notation and the proving of theorems. Instead, we emphasize concepts and useful “tricks of the trade” with many code segments, remarks, reminders, and warnings throughout the book.

The material of this book has been taught at different times to students in engineering, physics, computer science, and applied mathematics at Princeton University, Brown University, and MIT over the past fifteen years. Different segments have been taught to undergraduates and graduates, and to novices as well as to experts. To this end, on all three subjects covered, we start with simple introductory concepts and proceed to more advanced topics; such bandwidth, we believe, is one strength of this book.

COMPILATION GUIDE

For the purposes of this compilation example, we will assume that we are using the GNU g++ compiler to compile a C++ program we have written contained within the file myprog.cpp. We will also assume that the machine on which you are trying to compile is a parallel machine with some version of MPI installed. You will need to contact your system administrator to find out the exact version of MPI that is available and the paths on your local architecture. In the following examples, the argument following the “-o” flag designates the file name to be used for the output. If no “-o” option is specified, most compilers default to using the name “a.out.” We now present several different programming scenarios:

1. • No user-defined libraries or user-defined header files are needed, and no special system libraries (such as those associated with math.h) are needed other than the MPI libraries:

2. g++ -o myprog myprog.cpp -lmpi

3. • No user-defined libraries or user-defined header files are needed, but the special system library corresponding to math.h is needed along with the MPI libraries:

4. g++ -o myprog myprog.cpp -lmath -lmpi

5. • User-defined libraries, user-defined header files, and the special system library corresponding to math.h are needed along with the MPI libraries:

6. g++ -o myprog myprog.cpp -I/users/kirby/includes -L/users/ kirby/libs -lSCmathlib -lmath -lmpi

The string following the flag designates the location of the user-defined header files to be included. The string following the “-L” flag designates the location of the user-defined libraries to be included. The string “-lSCmathlib” links the program with the user-defined library we created, and the string “-lmath” links the program with the system math library corresponding to math.h. The string “-lmpi” links the MPI libraries.

In this chapter we apply the approximation theory we presented in Chapter 3 to find solutions of linear and nonlinear equations and to perform integration of general functions. Both subjects are classical, but they serve as basic tools in scientific computing operations and in solving systems of ordinary and partial differential equations. With regard to root finding, we consider both scalar as well as systems of nonlinear equations. We present different versions of the Newton-Raphson method, the steepest descent method, and the conjugate gradient method (CGM); we will revisit the latter in Chapter 9. With regard to numerical integration we present some basic quadrature approaches, but we also consider advanced quadrature rules with singular integrands or in unbounded domains.

On the programming side, we first introduce the concept of passing a function to a function; in the previous chapter we were passing variables. This allows an easy implementation of recursion, which is so often encountered in scientific computing. We offer several C++ examples from root finding and numerical integration applications that make use of recursion, and we show an effective use of classes and overloaded operators. We also address parallel programming with emphasis on domain decomposition, specifically the concept of reduction operations. We introduce the MPI commands MPI_Reduce and MPI_Allreduce for accomplishing reduction operations among a collection of processes.

This overview of Fortran 90 (F90) features is presented as a series of tables that illustrate the syntax and abilities of F90. Frequently, comparisons are made with similar features in the C++ and F77 languages and the Matlab environment.

These tables show that F90 has significant improvements over F77 and matches or exceeds newer software capabilities found in C++ and Matlab for dynamic memory management, user-defined data structures, matrix operations, operator definition and overloading, intrinsics for vector and parallel processors and the basic requirements for object-oriented programming.

They are intended to serve as a condensed quick-reference guide for programming in F90 and for understanding programs developed by others.

Introduction

The programming process is similar in approach and creativity to writing a paper. In composition, you are writing to express ideas; in programming, you are expressing a computation. Both the programmer and the writer must adhere to the syntactic rules (grammar) of a particular language. In prose, the fundamental idea-expressing unit is the sentence; in programming, two units – statements and comments – are available.

Composition, from technical prose to fiction, should be organized broadly, usually through an outline. The outline should be expanded as the detail is elaborated and the whole reexamined and reorganized when structural or creative flaws arise. Once the outline settles, you begin the actual composition process using sentences to weave the fabric your outline expresses. Clarity in writing occurs when your sentences, both internally and globally, communicate the outline succinctly and clearly. We stress this approach here with the aim of developing a programming style that produces efficient programs humans can easily understand.

To a great degree, no matter which language you choose for your composition, the idea can be expressed with the same degree of clarity. Some subtleties can be better expressed in one language than another, but the fundamental reason for choosing your language is your audience: people do not know many languages, and if you want to address the American population, you had better choose English over Swahili.