Numerical computations are a fundamental tool for engineers and scientists. The current practice of science and engineering demands that nontrivial computations be performed with both great speed and great accuracy. More and more, one finds that scientific insight and technologial breakthroughs are preceded by intense computational efforts such as modeling and simulation. It is clear that computing is, and will continue to be, central to the further development of science and technology.

As market forces and technological breakthroughs lowered the cost of computational power by several orders of magintude, there was a natural migration from large-scale mainframes to powerful desktop workstations. Vector processing and parallelism became possible, and this parallelism gave rise to a new collection of algorithms. Parallel architectures matured, in part driven by the demand created by the algorithms. Large computational codes were modified to take advantage of these parallel supercomputers. Of course, the term supercomputer has referred, at various times, to radically different parallel architectures. This includes vector processors, various shared memory architectures, distributed memory clusters, and even computational grids. Although the landscape of scientific computing changes frequently, there is one constant; namely, that there will always be a demand in the research community for high-performance computing.

When computations are first introduced in beginning courses, they are often straightforward “vanilla” computations, which are well understood and easily done using standard techniques and/or commercial software packages on desktop computers. However, sooner or later, a working scientist or engineer will be faced with a problem that requires advanced techniques, more specialized software (perhaps coded from scratch), and/or more powerful hardware.

High-performance computation was first realized in the form of SIMD parallelism with the introduction of the Cray and Cyber computers. At first these were single processor machines, but starting with the Cray XMP series, multiprocessor vector processors gained the further advantages of MIMD parallelism. Today, vector processing can be incorporated into the architecture of the CPU chip itself as is the case with the old AltiVec processor used in the MacIntosh.

The UNIX operating system introduced a design for shared memory MIMD parallel programming. The components of the system included multitasking, time slicing, semaphores, and the fork function. If the computer itself had only one CPU, then parallel execution was only apparent, called concurrent execution, nevertheless the C programming language allowed the creation of parallel code. Later multiprocessor machines came on line, and these parallel codes executed in true parallel.

Although these tools continue to be supported by operating systems today, the fork model to parallel programming proved too “expensive” in terms of startup time, memory usage, context switching, and overhead. Threads arose in the search for a better soluton, and resulted in a software revolution. The threads model neatly solves most of the low-level hardware and software implementation issues, leaving the programmer free to concentrate on the the essential logical or synchronization issues of a parallel program design. Today, all popular operating systems support thread style concurrent/parallel processing.

In this chapter we will explore vector and parallel programming in the context of scientific and engineering numerical applications. The threads model and indeed parallel programming in general is most easily implemented on the shared memory multiprocessor architecture.

Monte Carlo Methods for Optimization

Consider the problem of searching for the extremal values of an objective function f defined on a domain Ω, and equally important, for the points x ∈ Ω, where these values occur. An extremal value is called an optimum (maximum or minimum) while a point where an optimum occurs is called an optimizer (maximizer or minimizer).

If the domain is a subset of Euclidean space, we will assume f is differentiable. In this case gradient descent (or ascent) methods are used to locate local minima (or maxima). Whether or not a global extremum has been found depends upon the starting point of the search. Each local minimum (maximum) has its own basin of attraction and so it becomes a matter of starting in the right basin. Thus there is an element of chance involved if globally extreme values are desired.

On the other hand, we allow the possibility that Ω is a discrete, and possibly large, finite set. In this case downhill/uphill directional information is nonexistent and the search is forced to make due with objective values only. As the search proceeds from one point to the next, selecting the next point to try is often best left to chance.

A search process in which the next point or next starting point to try is randomly determined and may depend on the current location is, mathematically, a finite Markov Chain. Although the full resources of that theory may be brought to bear on the problem, only general assertions will be possible without knowing the nature of the specific objective function.

Many problems in scientific computation can be solved by reducing them to a problem in linear algebra. This turns out to be an extremely successful approach. Linear algebra problems often have a rich mathematical structure, which gives rise to a variety of highly efficient and well-optimized algorithms. Consequently, scientists frequently consider linear models or linear approximations to nonlinear models simply because the machinery for solving linear problems is so well developed.

Basic linear algebraic operations are so fundamental that many current computer architectures are designed to maximize performance of linear algebraic computations. Even the list of the top 500 fastest computers in the world (maintained at www.top500.org) uses the HPL benchmark for solving dense systems of linear equations as the main performance measure for ranking computer systems.

In 1973, Hanson, Krogh, and Lawson described the advantages of adopting a set of basic routines for problems in linear algebra. These basic linear algebra subprograms are commonly referred to as the Basic Linear Algebra Subprograms (BLAS), and they are typically divided into three heirarchical levels: level 1 BLAS consists of vector–vector operations, level 2 BLAS are matrix–vector operations, and level 3 BLAS are matrix–matrix operations. The BLAS have been standardized with an application programming interface (API). This allows hardware vendors, compiler writers, and other specialists to provide programmers with access to highly optimized kernel routines adapted to specific architectures. Profiling tools indicate that many scientific computations spend most of their time in those sections of the code that call the BLAS. Thus, even small improvements in the BLAS can yield substantial speedups.

The problem of solving linear systems of equations is central in scientific computation. Systems of linear equations arise in many areas of science, engineering, finance, commerce, and other disciplines. They emerge directly through mathematical models in these areas or indirectly in the numerical solution of mathematical models as for example in the solution of partial differential equations. Because of the importance of linear systems, a great deal of work has gone into methods for their solution.

A system of m equations in n unknowns may be written in matrix form as Ax = b in which the coefficient matrix A is m × n, while the unknown vector x and the right-hand side b are n-dimensional. The most important case is when the coefficient matix is square, corresponding to the same number of unknowns as equations. The more general m × n case can be reduced to this.

There are two main methods for solving these systems, direct and iterative. If arithmetic were exact, a direct algorithm would solve the system exactly in a predetermined finite number of steps. In the reality of inexact computation, a direct method still stops in the same number of steps but accepts some level of numerical error. A major consideration with respect to direct methods is mitigating this error. Direct methods are typically used on moderately sized systems in which the coefficient matrix is dense, meaning most of its elements are nonzero. By contrast, iterative methods are typically used on very large, sparse systems. Iterative methods asymptotically converge to the solution and so run until the approximation is deemed acceptable.

Directed Acyclic Graph Representation

Just as a graph is an effective way to understand a function, a directed acyclic graph is an effective way to understand a parallel computation. Such a graph shows when each calculation is done, which others can be done at the same time, what prior calculations are needed for it and into what subsequent calculations it feeds.

Starting from a directed acyclic graph and given a set of processors, then a schedule can be worked out. A schedule assigns each calculation to a specific processor to be done at a specified time.

From a schdule, the total time for a computation follows and, from this, we get the difficulty or complexity of the computation.

A Directed Acyclic Graph Defines a Computation

A computation can be accurately depicted by means of a directed acyclic graph (DAG), G = (N, A), consisting of a set of vertices N and a set of directed arcs A. In such a portrayal the vertices, or nodes, of the graph represent subtasks to be performed on the data, while the directed arcs indicate the flow of data from one subtask to the next. In particular, a directed arc (i, j) ∈ A from node i to j indicates that calculation j requires the result of calculation i.

The input data is shown at the top of the graph. We take the flow of data from the top to the bottom of the graph (or, less often, from left to right). This also becomes the flow of time; it follows that the graph can have no cycles. With this convention, we may omit the direction indicators on the arcs.