Be aware of algorithms that can be useful to you. There are many textbooks on algorithms and data structure. One of the more encyclopedic is the book by Cormen, Leiserson, Rivest, and Stein [23].

Scientific and engineering software almost always needs to be efficient in both time and memory. You should first consider optimizing at a high level, choosing data structures and algorithms that are inherently efficient. At a lower level, you should understand how computers work, and how to efficiently use modern computer architectures.

Choosing good algorithms and data structures is the foundation of designing efficient software. Without this, no amount of low-level tweaking of the code will make it run fast. This is especially true for large problems. By tuning your code you could gain an improvement of anything from, say, a factor of two to a factor of ten, in the speed of an algorithm. But moving from an algorithm that is O(n3) in time and O(n2) in memory to one that is O(n2) in time and O(n) in memory can give you a much greater benefit, especially if n ≍ 10 000 or larger. Sometimes you can get approximate algorithms that are fast, but the speed will depend on how accurately you want the answer. Here you need to know the problem well in order to see how large an error can be tolerated.

Numerical algorithms

Since the development of electronic digital computers (and even well before then) there has been a great deal of development of numerical algorithms. Many are highly efficient, accurate, and robust.

Dynamic memory allocation

In the middle of a routine you realize that you need some scratch space n floating point numbers long, and n is a complicated function of the inputs to the routine. What do you do? Do you add another scratch space argument to your routine along with its length as an argument and check that it is big enough (stopping the program immediately if it is not)? Often a better idea is to allocate the memory needed. In Fortran 90 this is done using the allocate command; in C you use the malloc function; in C++ and Java the new operator will do this task; in Pascal the allocation is done by the new operator, but the syntax is different from C++ or Java. These can be used to dynamically allocate memory. All of these commands return, or set, a pointer to the allocated block of memory.

The allocated block of memory is taken from a global list of available memory which is ultimately controlled by the operating system. This block of memory remains allocated until it is no longer accessible (if garbage collection is used), explicitly de-allocated, or the program terminates. So dynamically allocated memory can be used to hold return values or returned data structures. Dynamically allocated memory can be passed to other routines, and treated like memory that has been statically allocated, or allocated on a stack in most respects.

The data structure that controls the allocation and de-allocation of memory is called a memory heap. A memory heap may contain a pair of linked lists of pointers to blocks of memory.

Picking algorithms vs. keyhole optimization

Optimizing the performance of software requires working with at least two different points of view. One is the global view where the questions are about how to structure the overall architecture of the system. Another view is how the individual routines are written, and even how each line of code is written. All of these views are important. Selecting good algorithms is as important as selecting good hardware and implementing algorithms efficiently.

When the term “optimization” is used in computing, it is often taken to mean something like picking compiler options or “writing tight code” that uses the least number of clock cycles or operations. This is clearly important in writing fast, effective software, but it is only a part of the process, which begins early in the design stage.

Usually the first part of the design stage is the design of the central data structures and databases that will be used. These should be chosen so that there are well-known efficient algorithms for handling these data structures, preferably with readily available implementations that perform correctly, reliably and efficiently. Then the algorithms to carry out the main tasks need to be selected. An important guide to selecting them is their asymptotic complexity or estimate of the time needed. However, this is not the only guide; see the last section of this chapter for more information about how to refine this information and to make sure that the asymptotic estimates are relevant in practice. The last part is the detailed design process where the algorithms are implemented, and this should be done with an eye on efficiency to ensure that the whole system works well.

The development of the Internet has contributed to the development of public libraries of scientific software. Some of this development has occurred through the efforts of individuals, through Internet collaborations (as with the development of Linux), and through government supported software development by academics and others (as with the development of LAPACK). There is a wide range of other software packages for scientific computing, many now written in C/C++, although Fortran and other languages are used for parts of many of these systems: PETSc (which supports the use of MPI for parallel computation), IML++ (an iterative methods library in C++), SparseLib++ (for handling sparse matrices in C++), and PLTMG (for solving partial differential equations).

In parallel with this, there has also been a tremendous development of commercial numerical software. Beginning in 1970 the Numerical Algorithms Group (NAG), based in the UK, developed libraries which have been sold commercially as the NAG libraries since 1976; the Harwell library was also developed in the UK; the IMSL libraries were developed commercially in the US. Another set of numerical libraries, called SLATEC, was developed by the Sandia and Los Alamos US National Laboratories and the US Air Force. These are available through netlib (see the next section). Perhaps the most spectacular example of commercial numerical software is the development of MATLAB. Initially a collection of Fortran 77 routines based on the early LINPACK and EISPACK libraries for dense matrix computations with a text interface, MATLAB has evolved into a full-featured interactive programming language with special support for numerical computation and scientific visualization.

The SGI Power C compiler (PCA) does not allow more threads than processors (cf. the document “Multiprocessing C Compiler Directives”). In this sense, programs execute like the fork() programming model.

The keyword critical corresponds most closely with mutex in that only one thread at a time can execute this code and all threads execute it. The keyword synchronize corresponds most closely with barrier in that all threads must arrive at this point before any thread can go on.

There are also additional directives. The directive one processor means that the first thread to reach this code executes it meanwhile other threads wait. After execution by the first thread, the code is skipped by subsequent threads. There is an enter gate and corresponding exit gate directive. Threads must wait at the exit gate until all threads have passed the matching enter gate.

Loops to run in parallel must be marked with the pfor directive. It takes the argument iterate (start index; number of times through the loop; increment/decrement amount).

A reduction variable is local to each thread and their contributions must be added in a critical section.

The need for speed. Since the beginning of the era of the modern digital computer in the early 1940s, computing power has increased at an exponential rate (see Fig. 1). Such an exponential growth is predicted by the well-known “Moore's Law,” first advanced in 1965 by Gordon Moore of Intel, asserting that the number of transistors per inch on integrated circuits will double every 18 months. Clearly there has been a great need for ever more computation. This need continues today unabated. The calculations performed by those original computers were in the fields of ballistics, nuclear fission, and cryptography. And, today these fields, in the form of computational fluid dynamics, advanced simulation for nuclear testing, and cryptography, are among computing's Grand Challenges.

In 1991, the U.S. Congress passed the High Performance Computing Act, which authorized The Federal High Performance Computing and Communications (HPCC) Program. A class of problems developed in conjunction with the HPCC Program was designated “Grand Challenge Problems” by Dr. Ken Wilson of Cornell University. These problems were characterized as “fundamental problems in science and engineering that have broad economic or scientific impact and whose solution can be advanced by applying high performance computing techniques and resources.” Since then various scientific and engineering committees and governmental agencies have added problems to the original list. As a result, today there are many Grand Challenge problems in engineering, mathematics, and all the fundamental sciences. The ambitious goals of recent Grand Challenge efforts strive to

• build more energy-efficient cars and airplanes,

• design better drugs,

• forecast weather and predict global climate change,

• improve environmental modeling,

• […]

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.