K. B. Williams received a B.A. in Mathematics from California State University of Sacramento in 1955. He is a veteran of forty years in the scientific applications programming trenches and is the founder of King Baker Enterprises, of Stillwater, Oklahoma. He consults in applied numerical techniques. He can be reached at 802 South Ridge Drive, Stillwater, OK 74074, or via e-mail at Kbwms@aol.com.
Of the weapons in my software arsenal garnered over the past 40 years none has had more use and influence than the one described in this article: matmpy. matmpy began as an assembly-coded function in 1961 when I was solving problems in satellite tracking. Later I dressed it up in Fortran II clothes, followed by Fortran IV and ANSI Fortran, and finally C. As a C function matmpy plays a prominent role in Fading Memory and Kalman Filtering codes. This function has also seen wide use in solving structural analysis problems.
Function matmpy can perform the following operations:
1. A x B (Matrix Product)
2. A x BT (T means transpose)
3. ATx B
4. ATx BT
5. Negative Product (any of 1-4 above)
6. C += any of 1-4 above (replace-add)
7. C -= any of 1-4 above (replace-subtract)
I call this a Magic Matrix Multiply function because it performs transposition by indexing rather than creating transposes first. As a result, transposes get the same treatment as regular matrices. This technique saves a lot of memory space and computation time.
The code appears in Listing 1. matmpy calculates two-dimensional matrix indexes directly, in as efficient a manner as one-dimensional indexes. Consequently, the matrix multiply loops are fast.
The first two if/else statements set up the row and column controls. If the caller specifies ml as AB or ABT (#defined in the header file mpydefs.h, Listing 2) matmpy sets controls for accessing elements of a in the regular fashion. If not (ml is ATB or A-BT), matmpy sets up to access matrix a as a transpose. Matrix b gets similar treatment depending on whether the user specified mI as AB or ATB for a straight b matrix, or as ABT or ATBT to request transposition of matrix b.
The remainder of the program is straightforward, with the loops indexing over matrix a on a row-by-row basis and forming vector products with matrix b. Once the product is formed, matmpy first determines whether a negative product is desired (m2 is MP or CMP). Finally, matmpy determines whether a simple matrix product is desired (m2 is P or MP); if not, matmpy performs a replace-add operation to get the final result.
Sample Application
Listing 3 shows a code fragment taken from a function that adds data to a covariance matrix using the matrix lemma [1]. The equations are as follows:
H = SM(R + MSMT)-1 S = (I - HM)S(I - HM)T+ HRHTwhere:S is the j-by-j covariance matrix to be updated
M is the n-by-j row observation matrix
R is the n-by-n observation covariance matrix
H is the j-by-n filter gain matrix
Q is the temporary storage area dimensioned j-by-j
W is the temporary storage area dimensioned j-by-j
The code in this example calls functions avgvrt and madsub, which are stories in themselves. Function argvrt inverts a (possibly augmented) matrix that could be embedded in a larger matrix. Function madsub adds or subtracts matrices in combinations of transposes as specified by the caller (similar to matmpy).
Usage Notes
1. Function matmpy does not check for compatibility of dimensions of A and B.2. The software that accompanies this article includes a test driver, tstmpy.c, and a makefile to go with it, both included on this month's code disk.
Reference
1. Norman Morrison. Introduction to Sequential Smoothing and Prediction (McGraw-Hill, 1969), pp. 465-467, equations 12.2.17 and 12.3.7.