## Obituaries: Ralph A. Willoughby

**September 6, 2002**

Ralph A. Willoughby, 1923-2001

On receiving a PhD in mathematics from Berkeley, he accepted a teaching position at Georgia Tech, which he held until 1955, working during the summers at the Oak Ridge National Laboratory. In 1955, he joined Babcock and Wilcox in Lynchburg, Virginia, and became involved in atomic energy research. Two years later he joined the small Mathematics Group at IBM Research.

An initial focus of his work at IBM was circuit simulation. Presented with "stiff" equations, he and Werner Liniger realized that this "multiscale" phenomenon was not unique to circuit models. Together, they devised A-stable implicit differential integration formulae for general time-dependent nonlinear differential equations that included free parameters the user could exploit to mitigate the effects of the different time scales. This concept of "exponential fitting" is recorded in their paper "Efficient integration methods for stiff systems of ordinary differential equations" (*SIAM Journal on Numerical Analys*is, Vol. 7, 1970, 47-66), which was often quoted in the 1970s.

These new integration methods addressed the issue of stiffness and allowed practical time-step sizes. Use of these methods to solve nonlinear circuit equations, however, required the repeated solution of large systems of equations *Ax* = *b*. Large, at that point in time, meant several thousand variables, and such systems took hours to solve, occupying entire computers and requiring backup storage.

Ralph and Werner Liniger observed that although these equations were large, they were also sparse, with the sparsity structure fixed throughout the simulation. If a solver code that worked only with nonzero quantities could be produced, they reasoned, the amount of computation and memory required might be reduced to manageable levels. They proceeded to write an internal IBM memo detailing these thoughts.

Initially, special-purpose solver codes were built for each circuit problem. Subsequently seizing the challenge, Fred Gustavson, a colleague of Ralph and Liniger at IBM, devised a clever symbolic triangular factorization program that automatically generated a Fortran code specific for any given problem. Their frequently referenced joint work (Gustavson, Liniger, Willoughby) appeared as "Symbolic generation of an optimal Crout algorithm for sparse systems of linear equations" (*Journal of the Association for Computing Machinery*, Vol. 17, 1970, 87-109). Prior to this research, large sparse matrices were typically handled by iterative techniques; direct methods were considered appropriate only for small full matrices.

While Ralph and his IBM colleagues were engaged in solving circuit simulation problems, they became aware of a body of sparse matrix knowledge and techniques being developed in certain other application areas, e.g., power systems analysis, linear programming, and structural analysis. This knowledge was buried in application-specific codes or was being published only in applications journals that were not known to the numerical analysis community.

The IBM researchers were convinced from their own work of the emerging importance of sparse matrix problems. Ralph, following a path of the sort he pursued throughout his professional career, decided to organize a conference that would encourage people from these application areas to talk with numerical analysts and thereby stimulate research on this important new topic. He enlisted the financial and scientific support of IBM Research and organized the first Symposium on Sparse Matrices and Their Applications. The meeting was held at the IBM T.J. Watson Research Center in Yorktown Heights, New York, September 9-10, 1968.

The proceedings of this conference appeared only as an IBM Research Report, *Sparse Matrix Proceedings* (RA1 3-12-69). By 1970 several thousand requests for copies had come in from all over the world.

Ralph was a keynote speaker at a subsequent sparse matrix conference, organized by the Institute of Mathematics and Its Applications and held at Oxford, April 5-8, 1970. The proceedings of that conference were published as *Large Sparse Sets of Linear Equations* by Academic Press in 1970, with John Reid as editor.

A third meeting, organized jointly by Ralph and Donald Rose, was again held at the IBM T.J. Watson Research Center. At this meeting, held September 9-10, 1971, many of the talks were oriented toward algorithms for working with sparse matrices and toward direct methods for solving large sparse systems of linear equations. The meeting was part of the IBM Research Symposia Series, initiated in 1970, and the proceedings, *Sparse Matrices and Their Applications*, were published by Plenum Press in 1972, with Ralph and Donald Rose as co-editors.

During these years, Ralph was something of a pioneer in promoting the association of graphs and matrices. His work in this area focused on matrix reducibility and orderings, as in the IBM Research Technical Report *A Characterization of Matrix Irreducibility via Triangular Factorization* (RC 3931, 1972).

These sparse matrix conferences were remarkable in that many of the issues discussed are still considered important today. Examples include the effects of accessing various levels within a memory hierarchy on the efficiency of the computations, and the design of algorithms that are hierarchy-aware; the importance of matching data storage patterns to the pattern of use of that data within the algorithm, now known as "data locality," and blocking and reuse of data; the use of multiple functional units to speed computation; the need for computer architects to work with numerical analysts; the key role that preservation of structure can play in the repeated use of linear solvers with the nonlinear systems of equations that arise in most applications; reorderings to preserve sparsity; and the desirability of fast sparse vector operations, as in the sparse blas. I think it is fair to say that these meetings had a tremendous impact on the numerical linear algebra community and on its research agenda.

In 1973 Ralph organized another in the series of IBM Research Symposia, this time on stiff differential equations. The meeting was held in Wildbad, Germany, in October 1973, under the auspices of IBM Europe. The proceedings were published by Plenum Press.

In the latter part of the 1970s, Ralph's interests shifted to large-scale eigenvalue problems. We began working together and jointly developed the idea of using variants of Lanczos recursions as the basis for a series of algorithms for eigenvalue problems of various types. This work differed from other work in the literature by its refusal to conform to the accepted way of designing eigenvalue algorithms in terms of strictly orthogonal projections. The consequences of this approach, the benefits to be gained from its use, and discussions of practical implementations can be found in the two-volume book *Lanczos Algorithms for Large Symmetric Eigenvalue Computations* (Cullum and Willoughby, Birkhaüser Boston, 1985). These algorithms have been used in a wide variety of problems in computational physics and chemistry. Volume 1 will be republished this fall as part of the SIAM Classics Series. In 1985, Ralph and I organized an IBM European Institute Workshop titled Large Scale Eigenvalue Problems, which was held July 8-12, 1985, in Oberlech, Austria. Our paper, "A practical procedure for computing eigenvalues of large sparse nonsymmetric matrices," which demonstrated that a great deal of information about the spectrum of any nonsymmetric problem can be obtained by a simple generalization of the ideas used in our symmetric algorithms, was published in the proceedings of this workshop, *Large Scale Eigenvalue Problems*, which we co-edited (*Mathematical Studies*, Vol. 127, 1986, North-Holland, 193-240).

In 1989, working with Wolfgang Kerner, we extended this work to generalized eigenvalue problems for computing the Alfven spectrum of magneto-hydrodynamic models (see "A generalized nonsymmetric Lanczos procedure," *Computer Physics Communications*, Vol. 53, 1989, 19-48). In our paper "Computing eigenvalues of large matrices, some Lanczos algorithms and a shift and invert strategy" (in *Advances in Numerical Partial Differential Equations and Optimization*, S. Gómez, J.P. Hennart, R.A. Tapia, eds., SIAM, 1991, 198-246), we discussed the selection of suitable shifts for use in shift-and-invert eigenvalue algorithms.

In 1991, Ralph retired from IBM, and he and Nona moved to Walnut Creek to be near their two daughters, Carol and Ann.

Ralph will be remembered as a gentle person whose primary concern was never himself, but the welfare of his colleagues and family. He showed a particular interest in students, treating them as equal colleagues, asking questions to make them think about what they were doing, and showing interest in their progress by introducing them to leaders in the field and mentioning their work to others in the field.

A history buff, he seemed to know all there was to know about the history of numerical linear algebra, recent wars, and college football.

Ralph is survived by Nona, his wife of 54 years, his daughters Carol and Ann Willoughby, grandsons James and Paul Abers, two nephews (the sons of his late sister Eleanor McFadden), and by his brother, Peter Putnam, and his two children.

Contributions in his memory can be sent to the Scholarship Fund at the University of California, Berkeley.

*---Jane Cullum, Los Alamos National Laboratory, with the help of Iain Duff, Fred Gustavson, Werner Liniger, and Ann and Nona Willoughby.*