When I first began the class Matrix, I hoped to be able to leverage this similarity by using calls to the SimpleVec functions. For example, to add two vectors, you simply add their elements, two by two. The same is true of matrices. It seemed a shame to repeat the same operation twice, when a call to vAdd( ) works for both cases.

Other operations, however, are completely different for the two cases. There are two multiplication operators for vectors, only one for matrices, and the algorithms for multiplying them are different. Vector math has no function of division; matrix math does (sort of). In the end, I found it much cleaner to separate the primitive operations into two files, SimpleVec and SimpleMat. So what if a couple of functions look alike? The code we're talking about is tiny either way.

This month, I've made a few additions to SimpleMat. First, I noticed that I'd left out comparison functions--functions that were in SimpleVec. These are the functions mIsEqual( ) and mIsUnequal( ). They work in the way you might expect: two matrices are equal if and only if all their elements are equal.

The usual caveat concerning comparison of floating point numbers apply. For those who don't know (yes, there are plenty), many numbers can't be represented exactly in floating point, so when you do arithmetic with them, you may be surprised that the result isn't quite what you expect. Thus:

(1)

because 1/3 cannot be represented exactly in floating point format. Indeed, no fraction can, except negative powers of two. Comparing real numbers for equality is only reliable if there is no intervening arithmetic to cause roundoff error.

In the case of the function mIsEqual, I still lean back on the vector function. The code is simply:

// Test for equality
bool mIsEqual(const double a[ ], const double
    b[ ], const size_t m, const size_t n)
{
  return vIsEqual(a, b, m*n);
}

The vector class also had tests for >, < , etc. They were based on the absolute value (norm) of the vector. We don't attempt similar tests for matrices, because the concept of one matrix being numerically (not dimensionally) larger than another is so fuzzy. I suppose we could use the determinant as a measure of size, but I can't think of a good reason to do so. Here, I'm going to invoke the Plauger doctrine.

I also noticed that I'd left out the transpose operator. This is important, because we use the operation a lot. To transpose a matrix, you simply exchange rows with columns, and vice versa. Mathematically, b_ji=a_ij, where i and j are the row and column indices, respectively. The code reflects this relationship, while taking into account the conversion from matrix indices to simple array index. The key line in the code is:

b[j*m+i] = a[i*n+j];

where m and n are the row and column dimensions of A.

Be warned that you can't transpose a matrix back on itself. That is, you can't write:

mTrans(a, a, m, n);

Doing so scrambles the matrix that you're working on.

Consider, however, what happens if the matrix is square. The diagonal elements don't move during the tranposition, and the off-diagonals only move two-by two, exchanging them symmetrically across the diagonal. It's definitely possible to write a function that will transpose a square matrix in place, though the indexing is a bit tricky.

Everything goes all to heck, though, if the matrix isn't square. The process goes from exchanging elements, two by two, to something more akin to rotating a Rubik's cube. It can be done, but the logic is just horrible, and not at all recommended. In the end, I chose not to provide a function that only works for square matrices. I may revisit this decision in the future; we'll see.