The ACM Special Interest Group on Programming Languages, SIGPLAN, expects to hold the fourth in a series of conferences on the History of Programming Languages in 2020, see HOPL-IV. The first drafts of papers are to be submitted by August, 2018. That long lead time gives me the opportunity to write a detailed history of MATLAB. I plan to write the paper in sections, which I'll post in this blog as they are available. This is the second such installment.
LINPACK and EISPACK are Fortran subroutine packages for matrix computation developed in the 1970's. They are a follow up to the Algol 60 procedures that I described in my previous blog post about the Wilkinson and Reinsch Handbook. They led to the first version of MATLAB. This post is about EISPACK. A post about LINPACK will follow.
In 1970, even before the Handbook was published, a group at Argonne National Laboratory proposed to the U.S. National Science Foundation (NSF) to "explore the methodology, costs, and resources required to produce, test, and disseminate high-quality mathematical software and to test, certify, disseminate, and support packages of mathematical software in certain problem areas."
Argonne is a national research laboratory west of Chicago that has origins in the University of Chicago's work on the Manhattan project in the 1940s. In 1970 it was under the auspices of what was then called the Atomic Energy Commission (now the Department of Energy). As far as I know, this was the first time anybody at Argonne proposed a research project to the NSF.
At the time I was a summer consultant at Argonne and was spending a sabbatical year at Stanford. Yoshiko Ikebe, from the University of Texas, was also a participant. So initially we called the project NATS, for NSF, Argonne, Texas and Stanford, but that was soon changed to National Activity to Test Software.
NATS took on two projects. EISPACK, Matrix Eigensystem Package, translated the Algol procedures for eigenvalues from the second part of the Handbook to Fortran and worked extensively on testing and portability. FUNPACK, led by Jim Cody, developed rational approximations and software to compute special mathematical functions, including Bessel functions, gamma functions and the error function.
The Algol 60 procedures in the Wilkinson and Reinsch Handbook provide an excellent reference for the algorithms of matrix computation, but it was difficult to actually run the code because Algol compilers were not readily available.
In the 1970's the clear choice of a language for technical computation was Fortran. The EISPACK team rewrote the Algol in Fortran, being careful to retain the structure and the numerical properties. Wilkinson visited Argonne for a few weeks every summer. He didn't do any Fortran programming, but he provided advice about numerical behavior and testing procedures.
In 1971 the first release of the collection was ready to be sent to about two dozen universities and national laboratories where individuals had shown an interest in the project and agreed to test the software.
The Fortran subroutines in the first release were all derived from the Algol procedures in the second section of the Handbook, so the contents of the package was essentially the same as the contents I listed in my previous blog. Most of the subroutine names describe the underlying algorithm, not the problem being addressed. It is hard to tell from names like TRED1 and HQR2 what the subroutine actually does. There were 34 subroutines in the first release.
With seven authors of the users' guide for the first release, it took almost three years, until 1974, to complete the documentation. Since the project was considered to be language translation, the EISPACK Guide was published by Spring-Verlag, the publisher of the Handbook.
By 1976 a second release was ready for public distribution. The number of subroutines nearly doubled, to 70. The additional capabilities included
- The SVD algorithm of Golub, Kahan, and Reinsch for factoring rectangular matrices and solving overdetermined linear systems. The Algol procedures are in the first section of the Handbook.
- The QZ algorithm of Moler and Stewart for the generalized eigenvalue problem, $Ax = \lambda Bx$. We developed this in Fortran at the same time as the rest of EISPACK, so there was never any Algol version.
- A set of "driver routines" with names like RS, for real symmetric, and RSG, for real symmetric generalized, that reflect the problem being solved. These drivers collect the recommended subroutines from the underlying collection into one call.
A second Springer book, called the EISPACK Guide Extension, was published in 1976 to document the additional subroutines.
The hardware landscape that EISPACK faced initially was a diverse collection of mainframes from different manufacturers. The early 1970s was before the age of minicomputers. (The DEC VAX-11/780 was introduced in October, 1975.) It was also long before the IEEE 784 standard for floating point arithmetic. Each manufacturer had different word length and floating point format.
The most prevalent machines were from the IBM 360-370 family. These offered two floating point options. The 32-bit word with four bytes was designated REAL*4 and the 64-bit word was REAL*8. In contrast to today's floating point with the same word length, the 360 employed base 16. So all numbers between powers of 16, not powers of 2, were equally spaced. Jim Cody referred to this as "wobbling precision".
The Fortran distribution was available in several versions that differed in the setting of the machine accuracy parameter, MACHEP. Here are the different machines and corresponding values of MACHEP.
|IBM 360-370, REAL*4||16.**(-5)|
|IBM 360-370, REAL*8||16.D0**(-13)|
Algol does not have a complex data type, so the derived subroutines do not use Fortran complex arithmetic . Real and imaginary parts of complex quantities are stored in separate arrays.
The distributed software does not have a copyright notice or involve any kind of license. The term "open source" was not yet widely used. The closest we came to a license agreement was a statement about "certification". Both guides listed the machines, operating systems, and compilers of the 15 test sites. About half of the sites were universities and half were government labs. One was from Canada and one from Sweden. Various models from six different manufactures were involved.
This list of test sites was followed by the following statement about support.
Certification implies the full support of the NATS project in the sense that reports of poor or incorrect performance on at least the computer systems listed above will be examined and any necessary corrections made. This support holds only when the software is obtained through the channels indicated below and has not been modified; it will continue to hold throughout the useful life of EISPACK or until, in the estimation of the developers, the package is superseded or incorporated into other supported, widely-available program libraries. The following individual serves as a contact about EISPACK, accepting reports from users concerning performance:
This was followed by Burt Garbow's mailing address and phone number at Argonne. As far as I know, Burt never received any serious bug reports or complaints.
J. Dongarra and C. Moler, "EISPACK, A Package for Solving Matrix Eigenvalue Problems", chapter 4 in W. Cowell, Sources and Development of Mathematical Software, Prentice Hall, 1984, 10 pages, PDF available.
B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y. Ikebe, V. C. Klema, and C. B. Moler, Matrix Eigensystem Routines - EISPACK Guide, 2nd Edition, vol. 6 in Lecture Notes in Computer Science, Springer, 1974, 551 pages, EISPACK Guide.
B. S. Garbow, J. M. Boyle, J. J. Dongarra and C. B. Moler, Matrix Eigensystem Routines - EISPACK Guide Extension, vol. 51 in Lecture Notes in Computer Science, Springer, 1977, 343 pages, EISPACK Guide Extension.
Get the MATLAB code
Published with MATLAB® R2017a