A small C matrix library.

I am taking both Linear Algebra and Numeric Analysis this quarter and the latest assignment in Numerics is crossing over into Linear territory.  We have to calculate Spline interpolation, which involves solving a system of equations equal to n+1 in   size, which is a linear matrix of n+1 by n+1 in size.  So first I needed a matrix library that I could use to handle that part of the problem.

The code is here:
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.c
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix.h
https://github.com/BuckRogers1965/Examples/blob/master/Math/NumericAnalysis/matrix_testharness.c


Most of it is very strait forward and only took about an hour or two to write. The one section that took the most time was the part that built the inverse.

5 hours worth of coding.

This bit was difficult because the first way I tried to solve the problem I had 4 layers of for loops and it was just too confusing to figure out. So I deleted all that code and  tried again.

I broke the problem down into elementary row operations and had two layers of for loops in each layer, which I could hold in my mind to solve this function.

This section may be worth putting into two
different functions so I can reuse the parts to solve other problems.

I also want to add in a check here to test the inverse against the original matrix to make sure it is works.  If it doesn't work, then I should free the used memory and return a NULL.

Tested the matrix inverse function against another system to make sure it worked

After the inverse was working I wrote the multiply function.  I tested mult first by multiplying a matrix against the identity matrix and corrected a few issues with that.  When you multiply by the identity you should get the original matrix.   Once that was working I multiplied the matrix by its inverse and should have gotten the identity matrix, but I did not.  It took me a bit to figure out what I had reversed and then it was working.  After that matrix addition, and scalar multiplication was easy.

I converted the one file into three files so that the matrix.c  can be used as a library.  If you look at the matrix_testharness.c file, you can see that main() is located there, and that out of all the files there is only ever a single main().  You can see that at the top of matrix_testharness.c file there is a #include "matrix.h" file which brings in all the function headers.   The command to compile this 3 file example is here:

gcc matrix.c matrix_testharness.c

If you look at the top of matrix.h you can see that the data type for matrix is mapped into void *.  This keeps the contents of all the matrix objects private so they can only be accessed in the matrix.c file. This is important because you want a single interface to each of your objects.  If you let developers peek into the data structures they will start writing into them all nimbly bimbly and soon half the logic for the object is fuzzed across 50 files, with code copy pasted almost randomly between them all.  A single interface keeps things clean and is testable.  If there is a bug that needs fixed, it can be fixed in a single place.

I went through and cleaned up every warning.  I used -Wall in gcc to get most of them, then cppcheck and splint to find the rest.

The determinant

I added in a recursive function to calculate the determinate.  First of all, only square matrices can be calculated.  If the matrix is 2x2 then it returns ad-cb.  If the matrix is larger it picks the top row, and iterates over the columns, multiplying the cell in the top row against a new matrix made by removing the row and column you are currently in.  The sign alternates each time, starting positive.  The cell from the top line you are on is multiplied by the sign and the determinate of the new smaller matrix.

There is an alternative method that doesn't use recursion documented here:  http://www.purplemath.com/modules/determs2.htm  This looks fairly strait forward. and can be done in just a few loops without recursion.  This might be important for a 1000 x 1000 matrix, that would have to have a stack of 1000 function calls the original way.

I implemented this second  way of calculating the determinate. Because it is not calling N! subroutines.  However there is a problem.  Apparently this only works for 3x3 matrices. So switching back to the other method for now.

So, a few days later...

I ran this library against a 21 point spline interpolation and the determinate calculation took 2/10ths of a second.  Without the calculation the code ran in .007 seconds.  So I realized that the fractions I was pulling out of each line were the determinate, accumulated it and saved it to a new det field I created in the matrix object. This cached value is NAN until it is set.  The inverse matrix det is set to the original matrices 1/det.  If you multiply by a scalar the det increases by x^r, where x is the scalar and r is the number of rows, but only for a square matrix.   I rewrote Mat_Determinate() to use a more light weight variety of the same algorithm I used to calculate the inverse, stopping when I am in upper triangular form, with the pivots all one and the row multiplier options used to perform this accumulated into the det variable.

To do:

Write an automated test harness to test all the functions.  Worked toward that by creating a .h file and the start of a test harness file.

I noticed that the way I am trying to find the inverse requires that the pivot cannot contain a zero.  So I need to add in a recursive search to swap rows around to make this happen. So I also need the elementary row operation of swapping rows.  This needs to look down at the rows under the current row.  Because we already checked the determinate, we know there is a solution. 

Check for memory leaks.