Alan H. Karp

A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides

Explicit Runge-Kutta methods (RKMs) are among the most popular classes of formulas for the approximate numerical integration of nonstiff, initial value problems. However, high-order Runge-Kutta methods require more function evaluations per integration step than, for example, Adams methods used in PECE mode, and so, with RKMs, it is expecially important to avoid rejected steps. Steps are often rejected when certain derivatives of the solutions are very large for part of the region of integration. This corresponds, for example, to regions where the solution has a sharp front or, in the limit, some derivative of the solution is discontinuous. In these circumstances the assumption that the local truncation error is changing slowly is invalid, and so any step-choosing algorithm is likely to produce an unacceptable step. In this paper we derive a family of explicit Runge-Kutta formulas. Each formula is very efficient for problems with smooth solution as well as problems having rapidly varying solutions. Each member of this family consists of a fifty-order formula that contains imbedded formulas of all orders 1 through 4. By computing solutions at several different orders, it is possible to detect sharp fronts or discontinuities before all the function evaluations defining the full Runge-Kutta step have been computed. We can then either accpet a lower order solution or abort the step, depending on which course of action seems appropriate. The efficiency of the new algorithm is demonstrated on the DETEST test set as well as on some difficult test problems with sharp fronts or discontinuities.

Implementing linear algebra algorithms for dense matrices on a vector pipeline machine

This paper examines common implementations of linear algebra algorithms, such as matrix-vector multiplication, matrix-matrix multiplication and the solution of linear equations. The different versions are examined for efficiency on a computer architecture which uses vector processing and has pipelined instruction execution. By using the advanced architectural features of such machines, one can usually achieve maximum performance, and tremendous improvements in terms of execution speed can be seen over conventional computers.

A comparison of 12 parallel FORTRAN dialects

A simple program that approximates pi by numerical quadrature is rewritten to run on nine commercially available processors to illustrate the compilations that arise in parallel programming in FORTRAN. The machines used are the Alliant FX/8, BBN Butterfly, Cray X-MP/48, ELXSI 6400, Encore Multimax, Flex/32, IBM 3090/VF, Intel iPSC, and Sequent Balance. Some general impediments to using parallel processors to do production work are identified.<<ETX>>

Radiative transfer through an arbitrarily thick, scattering atmosphere

Abstract A method is presented for solving the equation of radiative transfer in a vertically inhomogeneous planetary atmosphere. The method, based on the spherical harmonics expansion, can be used to compute models with an arbitrarily large optical thickness and any scattering phase function. It is extremely efficient, requiring the equivalent of only two matrix multiplications per layer. This efficiency combined with its stability makes the method useful for computing realistic models of planetary atmospheres. To illustrate the range of validity of this method, we compute the plane albedo from model atmospheres containing clouds with optical thicknesses ranging from 0 to 106.

Efficient computation of spectral line shapes

Abstract A method for computing spectral line shapes based on the Fourier representation of the Voigt profile is presented. The linearity of Fourier transforms is used to eliminate the need to evaluate the Voigt function repeatedly. Instead, the method requires computing only sines, cosines, and exponentials for each line and two Fourier transforms for the entire frequency range. The method is shown to be computationally more efficient than previous approaches for a wide variety of problems, particularly when there are many overlapping lines. In contrast to other approximations which have fixed error bounds, the accuracy of the transform method can be controlled by varying the frequency interval and number of points. A part of the R -band of molecular oxygen at 13,160 cm -1 (7620A) is computed to demonstrate the advantages of the transform method.

Computing the angular dependence of the radiation of a planetary atmosphere

Abstract The numerical methods used for solving the transfer equation for a planetary atmosphere often give the specific intensity at a predetermined set of zenith angles, often the Gauss quadrature points. Since these ordinates are not always the ones at which the intensities are needed, some form of interpolation is needed. First, we show that the intensities computed from the moments derived from the spherical harmonics method are identical to those computed from the discrete ordinates method at the Gauss quadrature points. After a discussion of the various contributions to the errors, we compare three interpolation schemes. By far the most accurate way to interpolate the intensity curve is to use the integration of the source function method. We show how this method can be used with a minimum of computational effort. Although less accurate, interpolations are much faster. Not surprisingly, we find that Lagrange interpolation can produce very inaccurate values of the intensities and should not be used. Spline interpolation, on the other hand, always produces reasonable intensities. The accuracy of the spline interpolation can be improved greatly by supplementing the intensities at the Gauss points with values computed by the integration of the source function method at a few key ordinates.

Speeding up N-body Calculations on Machines without Hardware Square Root

The most time consuming part of an N-body simulation is computing the components of the accelerations of the particles. On most machines the slowest part of computing the acceleration is in evaluating

Data merging for shared-memory multiprocessors

r^{-3/2}

On the spherical harmonics and discrete ordinates methods for azimuth-dependent intensity calculations

, which is especially true on machines that do the square root in software. This note shows how to cut the time for this part of the calculation by a factor of 3 or more using standard Fortran.

1988 Gordon Bell Prize

An efficient software cache consistency mechanism for shared-memory multiprocessors that supports multiple writers and works for cache lines of any size is described. The mechanism relies on the fact that for a correct program only the global memory needs a consistent view of the shared data between synchronization points. The delayed consistency mechanism allows arbitrary use of data blocks between synchronizations. In contrast to other mechanisms, the mechanism needs no modification to the processor hardware or any assistance from the programmer or compiler. The processors can use normal cache management policies. Since no special action is needed to use the shared data, the processors are free to act almost as if they are all running out of a single cache. The global memory units are nearly identical to those on currently available machines. Only a small amount of hardware and/or software is needed to implement the mechanism. The mechanism can even be implemented using network-connected workstations.<<ETX>>