Beautiful Code [207]
I have used this code since I was young, and I am still fascinated by it. The beauty of this code is that it contains many complicated mathematical algorithms that are very well tested and difficult to reproduce. You can reuse it and rely on it without having to worry about whether it works.
The library's high accuracy and reliability are why—even after all these years—it remains the mathematical library of choice for anyone who needs to accurately and reliably solve equations. There are other mathematical libraries available, but they can't compete with the CERN library's proven reliability and accuracy.
The first part of this chapter discusses this code's outer beauty: the things that make it so accurate and reliable that developers want to reuse it and the elements that make this reuse as simple as possible. The second part explores the inner beauty of its implementation details.
The Long-Term Benefits of Beautiful Design > Outer Beauty
15.3. Outer Beauty
If you've ever tried to solve a system of linear equations or perform an equally complicated mathematical operation, you know that many times the code you write to achieve this does not deliver the correct results. One of the greatest problems with mathematical libraries is that rounding errors and floating-point operations lead to solution instabilities and incorrect results.
If you design a mathematical library, you need to carefully define the range in which each algorithm will work. You need to write each algorithm in such a way that it will adhere to these conditions, and you also need to write it in such a way that the rounding errors will cancel out. This can be very complicated.
In the CERN library, the algorithms are specified in a very precise way. Basically, if you look at any routine, you will notice that it has a description of what it is going to do. It really doesn't matter in which language the routine is written. In fact, these routines were written in Fortran but have interfaces that allow them to be called from almost any other place. That's also a beautiful thing. In some sense, the routine is a black box: you don't care what goes on inside, only that it delivers the appropriate results for your inputs. It carefully defines what every routine is doing, under which conditions it is working, what input data it accepts, and what constraints must be put on the input data in order to get the correct answer.
For example, let's look at the LAPACK library's SGBSV routine, which solves a system of linear equations for a banded matrix. If you try to solve a system of linear equations numerically, you use different algorithms. Different algorithms operate better in different domains, and you need to know the structure of the matrix to choose the best one. For instance, you would want to use one algorithm to solve the problem if you had a banded matrix (a matrix where most of the elements are around the diagonal), and a different one if you had a sparse matrix (a matrix that has a lot of zeros and few numbers).
Because different routines are optimized for different situations, the best routine to use really depends on the matrix structure that you have. However, in order to understand the range of this, you need to understand how to input data to these routines. Sometimes you input data in the form of a matrix. Sometimes—for instance, with a banded matrix— you send it like a very narrow array. Each of these routines and their requirements are described very clearly in the library:
Code View: Scroll / Show All
SUBROUTINE SGBSV( N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO )