Beautiful Code [209]
END IF
*
* Compute the LU factorization of the band matrix A.
*
CALL SGBTRF ( N, N, KL, KU, AB, LDAB, IPIV, INFO )
IF( INFO.EQ.0 ) THEN
*
* Solve the system A*X = B, overwriting B with X.
*
CALL SGBTRS( 'No transpose', N, KL, KU, NRHS, AB, LDAB, IPIV,
$ B, LDB, INFO )
END IF
RETURN
*
* End of SGBSV
*
END
One of the first things to notice in the code for the SGBSV routine is that it starts with a long comment that describes the routine's purpose and use. In fact, the comment is exactly the same as the manual page for that routine. Having the full documentation of the routine's usage in the code is important because it connects the routine's internal structure with its usage. In many other cases, I have found that the manual description and code documentation have nothing in common. I think this practice of marrying the two is one thing that makes code beautiful.
Following the initial comments, the algorithm that the routine uses is detailed in the routine's description. This helps anyone using the code to understand what the code will do and how it should react. Next comes a detailed description of the arguments, with their ranges explicitly specified.
The AB argument is an interesting one to consider. This argument contains the elements of the matrix A. Because the matrix is banded, it contains a lot of zero values, which are not clustered close to the diagonal. In principle, the input to the routine could be presented as a two-dimensional array of the dimensions of the matrix. However, this would be a waste of memory space. Instead, the AB argument contains only nonzero elements of the matrix next to its diagonal.
The format of the AB argument not only conserves memory space, but also has another purpose. In this routine, the algorithm is using the properties of the system of equations to solve the problem in a more efficient way. This means that the algorithm relies on the user to provide the correct matrix type as an input. If the AB argument contains all the elements of the matrix, one or more of the elements outside the band could accidentally be set to nonzero. This could lead to errors in the solution. The format chosen for AB makes it impossible to make this mistake. This was done on purpose, and it contributes to the code's beauty.
The AB argument also plays another role: it serves as an output argument as well as an input argument. In this context, the design solves a different problem. By having the routine reuse the space that the original program allocated, the developers of this code ensured that the routine would work as long as the original program had sufficient memory. If it had been written so that the routine needed additional memory allocation, then it might not run if the system was unable to allocate more memory. This can be especially problematic when there is a really large system of equations, and the routine needs a significant amount of memory space in which it can operate. The sample code is immune to such problems because it was written so that the routine can return the solution as long as the original program has sufficient memory space to store the problem. This is very beautiful.
Before I move on to the other arguments, I want to discuss this issue further. I have seen a lot of code written in my lifetime. Very often, developers write code and unconsciously place intrinsic restrictions on it. Most commonly, they restrict the size of the problem that can be solved. This occurs as a result of the following thought process:
I have something to write.
I'll write it fast so I can see if it works.
Once it works, I will generalize it for the real problem.
This process prompts developers to build restrictions into code, which very often leads to difficult-to-find errors that may take years to clean. During this process, developers commonly place explicit or implicit restrictions on the size of the problem they may solve. For example, an explicit restriction may be the definition of a large data space, which should be large enough for