Showing papers in "ACM Transactions on Mathematical Software in 2000"

A MATLAB differentiation matrix suite

Abstract: A software suite consisting of 17 MATLAB functions for solving differential equations by the spectral collocation (i.e., pseudospectral) method is presented. It includes functions for computing derivatives of arbitrary order corresponding to Chebyshev, Hermite, Laguerre, Fourier, and sinc interpolants. Auxiliary functions are included for incorporating boundary conditions, performing interpolation using barycentric formulas, and computing roots of orthogonal polynomials. It is demonstrated how to use the package for solving eigenvalue, boundary value, and initial value problems arising in the fields of special functions, quantum mechanics, nonlinear waves, and hydrodynamic stability.

Algorithm 799: revolve: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation

Abstract: In its basic form, the reverse mode of computational differentiation yields the gradient of a scalar-valued function at a cost that is a small multiple of the computational work needed to evaluate the function itself. However, the corresponding memory requirement is proportional to the run-time of the evaluation program. Therefore, the practical applicability of the reverse mode in its original formulation is limited despite the availability of ever larger memory systems. This observation leads to the development of checkpointing schedules to reduce the storage requirements. This article presents the function revolve, which generates checkpointing schedules that are provably optimal with regard to a primary and a secondary criterion. This routine is intended to be used as an explicit “controller” for running a time-dependent applications program.

Algorithm 806: SPRNG: a scalable library for pseudorandom number generation

Abstract: In this article we present background, rationale, and a description of the Scalable Parallel Random Number Generators (SPRNG) library. We begin by presenting some methods for parallel pseudorandom number generation. We will focus on methods based on parameterization, meaning that we will not consider splitting methods such as the leap-frog or blocking methods. We describe, in detail, parameterized versions of the following pseudorandom number generators: (i) linear congruential generators, (ii) shift-register generators, and (iii) lagged-Fibonacci generators. We briefly describe the methods, detail some advantages and disadvantages of each method, and recount results from number theory that impact our understanding of their quality in parallel applications. SPRNG was designed around the uniform implementation of different families of parameterized random number generators. We then present a short description of SPRNG. The description contained within this document is meant only to outline the rationale behind and the capabilities of SPRNG. Much more information, including examples and detailed documentation aimed at helping users with putting and using SPRNG on scalable systems is available at htt;//sprng.sc.fsu.edu. In this description of SPRNG we discuss the random-number generator library as well as the suite of tests of randomness that is an integral part of SPRNG. Random-number tools for parallel Monte Carlo applications must be subjected to classical as well as new types of empirical tests of randomness to eliminate generators that show defects when used in scalable envionments.

A simple method for generating gamma variables

Abstract: We offer a procedure for generating a gamma variate as the cube of a suitably scaled normal variate. It is fast and simple, assuming one has a fast way to generate normal variables. In brief: generate a normal variate x and a uniform variate U until In (U) 1. The gamma procedure is particularly fst for C implementation if the normal variate is generated in-line, via the #define feature. We include such an inline version, based on our ziggurat method. With it, and an inline uniform generator, gamma variates can be produced in 400MHz CPUs at better than 1.3 million per second, with the parameter α changing from call to call.

Algorithm 807: The SBR Toolbox—software for successive band reduction

Abstract: We present a software toolbox for symmetric band reduction via orthogonal transformations, together with a testing and timing program. The toolbox contains drivers and computational routines for the reduction of full symmetric matrices to banded form and the reduction of banded matrices to narrower banded or tridiagonal form, with optional accumulation of the orthogonal transformations, as well as repacking routines for storage rearrangement. The functionality and the calling sequences of the routines are described, with a detailed discussion of the “control” parameters that allow adaptation of the codes to particular machine and matrix characteristics. We also briefly describe the testing and timing program included in the toolbox.

A framework for symmetric band reduction

Abstract: We develop an algorithmic framework for reducing the bandwidth of symmetric matrices via orthogonal similarity transformations. This framework includes the reduction of full matrices to banded or tridiagonal form and the reduction of banded matrices to narrower banded or tridiagonal form, possibly in multiple steps. Our framework leads to algorithms that require fewer floating-point operations than do standard algorithms, if only the eigenvalues are required. In addition, it allows for space-time tradeoffs and enables or increases the use of blocked transformations.

PSBLAS: a library for parallel linear algebra computation on sparse matrices

Abstract: Many computationally intensive problems in engineering and science give rise to the solution of large, sparse, linear systems of equations. Fast and efficient methods for their soltion are very important because these systems usually occur in the innermost loop of the computational scheme. Parallelization is often necessary to achieve an acceptable level of performance. This paper presents the design, implementation, and interface of a library of Basic Linear Algebra Subroutines for sparse matrices (PSBLAS) which is specifically tailored to distributed-memory computers. PSBLAS enables easy, efficient, and portable implementations of parallel iterative solvers for linear systems. The interface keeps in view a Single Program Multiple Data programming model on distributed-memory machines. However, the architecture of the library does not exclude an implementation in different paradigms, such as those based on the shared-memory model.

NEOS and Condor: solving optimization problems over the Internet

Abstract: We discuss the use of Condor, a distributed resource management system, as a provider of computational resources for NEOS, an environment for solving optimization problems over the Internet. We also describe how problems are submitted and processed by NEOS, and then scheduled and solved by Condor on available (idle) workstations

PYTHIA-II: a knowledge/database system for managing performance data and recommending scientific software

Abstract: Often scientists need to locate appropriate software for their problems and then select from among many alternatives. We have previously proposed an approach for dealing with this task by processing performance data of the targeted software. This approach has been tested using a customized implementation referred to as PYTHIA. This experience made us realize the complexity of the algorithmic discovery of knowledge from performance data and of the management of these data together with the discovered knowledge. To address this issue, we created PYTHIA-II—a modular framework and system which combines a general knowledge discovery in databases (KDD) methodology and recommender system technologies to provide advice about scientific software/hardware artifacts. The functionality and effectiveness of the system is demonstrated for two existing performance studies using sets of software for solving partial differential equations. From the end-user perspective, PYTHIA-II allows users to specify the problem to be solved and their computational objectives. In turn, PYTHIA-II (i) selects the software available for the user's problem (ii) suggests parameter values, and (iii) assesses the recommendation provided. PYTHIA-II provides all the necessary facilities to set up database schemas for testing suites and associated performance data in order to test sets of software. Moreover, it allows easy interfacing of alternative data mining and recommendation facilities. PYTHIA-II is an open-ended system implemented on public domain software and has been used for performance evaluation in several different problem domains.

ADMIT-1: automatic differentiation and MATLAB interface toolbox

Abstract: ADMIT-1 enables the computation of sparse Jacobian and Hessian matrices, using automatic differentiation technology, from a MATLAB environment. Given a function to be differentiated, ADMIT-1 will exploit sparsity if present to yield sparse derivative matrices (in sparse MATLAB form). A generic automatic differentiation tool, subject to some functionality requirements, can be plugged into ADMIT-1; examples include ADOL-C (C/C++ target functions)and ADMAT (MATLAB target funcitons). ADMIT-1 also allows for the calculation of gradients and has several other related functions. This article provides an introduction to the design and usage of ADMIT-1.

Algorithms for the computation of all Mathieu functions of integer orders

Abstract: The article presents methods for the computation of all Mathieu functions of integer order, which cover a large range of n and h; previous algorithms were limited to small values of n. The algorithms are given in sufficient detail to enable straightforward implementation. The algorithms can handle a large range of the order n (0-200) and the parameter h (0-4n).

Algorithm 801: POLSYS_PLP: a partitioned linear product homotopy code for solving polynomial systems of equations

Abstract: Globally convergent, probability-one homotopy methods have proven to be very effective for finding all the isolated solutions to polynomial systems of equations. After many years of development, homotopy path trackers based on probability-one homotopy methods are reliable and fast. Now, theoretical advances reducing the number of homotopy paths that must be tracked, and in the handling of singular solutions, have made probability-one homotopy methods even more practical. POLSYS_PLP consists of Fortran 90 modules for finding all isolated solutions of a complex coefficient polynomial system of equations. The package is intended to be used in conjunction with HOMPACK90 (Algorithm 777), and makes extensive use of Fortran 90 derived data types to support a partitioned linear product (PLP) polynomial system structure. PLP structure is a generalization of m-homogeneous structure, whereby each component of the system can have a different m-homogeneous structure. The code requires a PLP structure as input, and although finding the optimal PLP structure is a difficult combinatorial problem, generally physical or engineering intuition about a problem yields a very good structure. POLSYS_PLP employs a sophisticated power series end game for handling singular solutions, and provides support for problem definition both at a high level and via hand-crafted code. Different PLP structures and their corresponding Bezout

Automatic sampling with the ratio-of-uniforms method

Abstract: Applying the ratio-of-uniforms method for generating random variates results in very efficient, fast, and easy-to-implement algorithms. However parameters for every particular type of density must be precalculated analytically. In this article we show, that the ratio-of-uniforms method is also useful for the design of a black-box algorithm suitable for a large class of distributions, including all with log-concave densities. Using polygonal envelopes and squeezes results in an algorithm that is extremely fast. In opposition to any other ratio-of-uniforms algorithm the expected number of uniform random numbers is less than two. Furthermore, we show that this method is in some sense equivalent to transformed density rejection.

Algorithm 805: computation and uses of the semidiscrete matrix decomposition

Abstract: We present algorithms for computing a semidiscrete approximation to a matrix in a weighted norm, with the Frobenius norm as a special case. The approximation is formed as a weighted sum of outer products of vectors whose elements are ±1 or 0, so the storage required by the approximation is quite small. We also present a related algorithm for approximation of a tensor. Applications of the algorithms are presented to data compression, filtering, and information retrieval; software is provided in C and in Matlab.

On stopping criteria in verified nonlinear systems or optimization algorithms

Abstract: Traditionally, iterative methods for nonlinear systems use heuristic domain and range stopping criteria to determine when accuracy tolerances have been met. However, such heuristics can cause stopping at points far from actual solutions, and can be unreliable due to the effects of round-off error or inaccuracies in data. In verified computations, rigorous determination of when a set of bounds has met a tolerance can be done analogously to the traditional approximate setting. Nonetheless, the range tolerance possibly cannot be met. If the criteria are used to determine when to stop subdivision of n-dimensional bounds into subregions, then failure of a range tolerance results in excessive, unnecessary subdivision, and could make the algorithm impractical. On the other hand, interval techniques can detect when inaccuracies or round-off will not permit residual bounds to be narrowed. These techniques can be incorporated into range thickness stopping criteria that complement the range stopping criteria. In this note, the issue is first introduced and illustrated with a simple example. The thickness stopping criterion is then formally introduced and analyzed. Third, inclusion of the criterion within a general verified global optimization algorithm is studied. An industrial example is presented. Finally, consequences and implications are discussed.

Mining and visualizing recommendation spaces for elliptic PDEs with continuous attributes

Abstract: In this paper we extend previous work in mining recommendation spaces based on symbolic problem features to PDE problems with continuous-valued attributes. We identify the research issues in mining such spaces, present a dynamic programming algorithm form the data-mining literature, and describe how a priori domain metaknowledge can be used to control the complexity of induction. A visualization aid for continuous-valued recommendation spaces is also outlined. Two case studies are presented to illustrate our approach and tools: (i) a comparison of an iterative and a direct linear system solver on nearly singular problems, and (ii) a comparison of two iterative solvers on problems posed on nonrectangular domains. Both case studies involve continuously varying problem and method parameters which strongly influence the choice of best algorithm in particular cases. By mining the results from thousands of PDE solves, we can gain valuable insight into the relative performance of these methods on similar problems.

Algorithm 800: Fortran 77 subroutines for computing the eigenvalues of Hamiltonian matrices. I: the square-reduced method

Abstract: This article describes LAPACK-based Fortran 77 subroutines for the reduction of a Hamiltonian matrix to square-reduced form and the approximation of all its eigenvalues using the implicit version of Van Loan's method. The transformation of the Hamiltonian matrix to a square-reduced form transforms a Hamiltonian eigenvalue problem of order 2n to a Hessenberg eigenvalue problem of order n. The eigenvalues of the Hamiltonian matrix are the square roots of those of the Hessenberg matrix. Symplectic scaling and norm scaling are provided, which, in some cases, improve the accuracy of the computed eigenvalues. We demonstrate the performance of the subroutines for several examples and show how they can be used to solve some control-theoretic problems.

Algorithm 804: subroutines for the computation of Mathieu functions of integer orders

Abstract: Computer subroutines in C++ for computing Mathieu functions of integer orders are described. The routines can handle a large range of the order n and the parameter h. Sample test results and graphs are given.

Band reduction algorithms revisited

Abstract: In this paper we explain some of the changes that have been incorporated in the latest version of the LAPACK subroutine for reducing a symmetric banded matrix to tridiagonal form. These modifications improve the performance for larger-bandwidth problems and reduce the number of operations when accumulating the transformations onto the identity matrix, by taking advantage of the structure of the initial matrix. We show that similar modifications can be made to the LAPACK subroutines for reducing a symmetric positive definite generalized eigenvalue problem to a standard symmetric banded eigenvalue problem and for reducing a general banded matrix to bidiagonal form to facilitate the computation of the singular values of the matrix.

Corrigendum: Algorithm 806: SPRNG: a scalable library for pseudorandom number generation

Abstract: In this article we present background, rationale, and a description of the Scalable Parallel Random Number Generators (SPRNG) library. We begin by presenting some methods for parallel pseudorandom number generation. We will focus on methods based on parameterization, meaning that we will not consider splitting methods such as the leap-frog or blocking methods. We describe, in detail, parameterized versions of the following pseudorandom number generators: (1) linear congruential generators, (ii) shift-register generators, and (iii) lagged-Fibonacci generators. We briefly describe the methods, detail some advantages and disadvantages of each method, and recount results from number theory that impact our understanding of their quality of parallel applications. SPRNG was designed around the uniform implementation of different families of parameterized random number generators. We then present a short description of SPRNG. The description contained within this document is meant only to outline the rationale behind and the capabilities of SPRNG. Much more information, including examples and detailed documentation aimed at helping users with putting and using SPRNG on scalable systems is available at http://sprng.cs.fsu.edu. In this description of SPRNG we discuss the random-number generator library as well as the suite of tests of randomness that is an integral part of SPRNG. Random-number tools for parallel Monte Carlo applications must be subjected to classical as well as new types of empirical tests of randomness to eliminate generators that show defects when used in scalable environments.

HPFBench: a high performance Fortran benchmark suite

Abstract: The high performance Fortran (HPF) benchmark suite HPFBench is designed for evaluating the HPF language and compilers on scalable architectures. The functionality of the benchmarks covers scientific software library functions and application kernels that reflect the computational structure and communication patterns in fluid dynamic simulations, fundamental physics, and molecular studies in chemistry and biology. The benchmarks are characterized in terms of FLOP count, memory usage, communication pattern, local memory accesses, array allocation mechanism, as well as operation and communication counts per iteration. The benchmarks output performance evaluation metrics in the form of elapsed times, FLOP rates, and communication time breakdowns. We also provide a benchmark guide to aid the choice of subsets of the benchmarks for evaluating particular aspects of an HPF compiler. Furthermore, we report an evaluation of an industry-leading HPF compiler from the Portland Group Inc. using the HPFBench benchmarks on the distributed-memory IBM SP2

Superconvergent interpolants for collocation methods applied to mixed-order BVODEs

Abstract: Continuous approximations to boundary value problems in ordinary differential equations (BVODEs), constructed using collocation at Gauss points, are more accurate at the mesh points than at off-mesh points. From these approximations, it is possible to construct improved continuous approximations by extending the high accuracy that is available at the mesh points to off-mesh points. One possibility is the bootstrap approach, which improves the accuracy of the approximate solution at the off-mesh points in a sequence of steps until the accuracy at the mesh points and off-mesh points is consistent. A bootstrap approach for systems of mixed-order BVODEs is developed to improve approximate solutions produced by COLNEW, a Gauss-collocation-based software package. An implementation of this approach is discussed and numerical results presented which confirm that the improved approximations satisfy the predicted error bounds and are relatively inexpensive to construct.

Note on generalization in experimental algorithmics

Abstract: A recurring theme in mathematical software evaluation is the generalization of rankings of algorithms on test problems to build knowledge-based recommender systems for algorithm selection. A key issue is to profile algorithms in terms of the qualitative characteristics of benchmark problems. In this methodological note, we adapt a novel all-pairs algorithm for the profiling task; given performance rankings for m algorithms on n problem instances, each described with p features, identify a (minimal) subset of p that is useful for assessing the selective superiority of an algorithm over another, for all pairs of m algorithms. We show how techniques presented in the mathematical software literature are inadequate for such profiling purposes. In conclusion, we also address various statistical issues underlying the effective application of this technique.

Accurate approximate solution of partial differential equations at off-mesh points

Abstract: Numerical methods for partial differential equations often determine approximations that are more accurate at the set of discrete meshpoints than they are at the “off-mesh” points in the domain of interest. These methods are generally most effective if they are allowed to adjust the location of the mesh points to match the local behavior of the solution. Different methods will typically generate their respective approximations on incompatible, unstructured meshes, and it can be difficult to evaluate the quality of a particular solution, or to visualize important properties of a solution. In this paper we will introduce a generic approach which can be used to generate approximate solution values at arbitrary points in the domain of interest for any method that determines approximations to the solution and low-order derivatives at meshpoints. This approach is based on associating a set of “collocation” points with each mesh element and requiring that the local approximation interpolate the meshpoint data and almost satisfy the partial differential equation at the collocation points. The accuracy associated with this interpolation/collocation approach is consistent with the “meshpoint accuracy” of the underlying method. The approach that we develop applies to a large class of methods and problems. It uses local information only and is therefore particularly suitable for implementation in a parallel or network computing environment. Numerical examples are given for some second-order problems in two and three dimensions.

Algorithm 802: an automatic generator for bivariate log-concave distributions

Abstract: Different automatic (also called universal or black-box) methods have been suggested to sample from univariate log-concave distributions. Our new automatic algorithm for bivariate log-concave distributions is based on the method of transformed density rejection. In order to construct a hat function for a rejection algorithm the bivariate density is transformed by the logarithm into a concave function. Then it is possible to construct a dominating function by taking the minimum of several tangent planes, which are by exponentiation transformed back into the original scale. The choice of the points of contact is automated using adaptive rejection sampling. This means that points that are rejected by the rejection algorithm can be used as additional points of contact. The article describes the details how this main idea can be used to construct Algorithm ALC2D that can generate random pairs from all bivariate log-concave distributions with known domain, computable density, and computable partial derivatives.

Remark on algorithm 746: new features of PCOMP, a Fortran code for automatic differentiation

Abstract: The software system PCOMP uses automatic differentiation to calculate derivatives of functions that are defined by the user in a modeling language similar to Fortran. This symbolical representation is converted into an intermediate code, which can be interpreted to calculate function and derivative values at run-time within machine accuracy. Furthermore, it is possible to generate Fortran code for function and gradient evaluation, which has to be compiled and linked separately. The first version of PCOMP was introduced in Dobmann et al. [1995]. In this article, we describe a series of extensions and additional features that have been implemented in the meantime.

Remark on Algorithm 761: scattered-data surface fitting that has the accuracy of a cubic polynomial

Abstract: Several improvements to the estimation of partial derivatives in Algorithm 761 are presented. The problems corrected are (1) in the calculation of the probability weight in subroutine SDPD3P which may result in overflow, (2) in the calculation of final weight in subroutine SDPD3P which may result in overflow, (3) in the computation of a determinant in subroutine SDLEQN which is not necessary, and (4) in the computation of the condition number of a matrix in subroutine SDLEQN which generates very different results for matrices that differ only in row order.

Algorithm 803: a simpler macro processor

Abstract: Macro processors have been in the computing tool chest since the late 1950's. Their use, though perhaps not what it was in the heyday of assembly language programming, is still widespread. In the past, producing a full-featured macro processor has required significant effort, similar to that required to implement the front-end to a compiler augmented by appropriate text substitution capabilities. The tool described here adopts a different approach. The text containing macro definitions and substitutions is, in a sense, “compiled” to produce a program, and this program must then be executed to produce the final output.

An observation on bisection software for the symmetric tridiagonal eigenvalue problem

Abstract: In this article we discuss a small modification of the bisection routines in EISPACK and LAPACK for finding a few of the eigenvalues of a symmetric tridiagonal matrix A. When the principal minors of the matrix A yield good approximations to the desired eigenvalues, these modifications can yield about 30% reduction in the computation times.

Editorial: special issue in honor of John Rice's 65th birthday

Abstract: This special issue of TOMS forms part of the proceedings of a symposium held May 22–26, 1999, at West Lafayette, Indiana, to honor John R. Rice for his many contributions to the field of computational science, as well as to computer sciences at Purdue University. The program reflected the breadth of Rice’s interests and accomplishments; there were tributes to Rice by his students, collaborators, and friends, perspectives on computational science, and research papers. Entitled the 1999 International Symposium on Computational Science, the symposium was organized by the Purdue University Department of Computer Sciences, and was sponsored by Working Group 2.5 (Numerical Software) of the International Federation for Information Processing (IFIP), and by the International Association for Mathematics and Computers in Simulation (IMACS). Over the years, John Rice enriched each of these organizations in many ways. Here we have collected some papers presented at the symposium in recognition of John Rice’s pivotal contributions to the field of mathematical software. His work in organizing the conferences Mathematical Software (1970) and Mathematical Software II (1974) at Purdue were key events which served to galvanize a community of researchers which continues to be productive to this day. In 1975, John Rice began an 18-year tenure as founding editor of the ACM Transactions on Mathematical Software. He remains as one of its Associate Editors. The symposium was held in conjunction with his 65th birthday. Those of us who have benefited from his insight, wisdom, mentoring, and presence have wanted for a long time to thank him and wish him the best for the future. We now take that opportunity with great pleasure.