Fleshing out the Matrix class - Embedded.com

Fleshing out the Matrix class

This month we continue with the process of developing a matrix class to accompany the vector class we developed earlier. As mathematical tools, vectors are a lot more useful when we can combine them with matrices. This is particularly true for the all-important matrix-vector product. In my last column, I gave you the skeleton of a matrix class. Now it's time to put flesh on the bones.

From class Vector , I've borrowed the notion of using class variables to define the default dimensions of a matrix. This important step lets us adjust the size of objects created by the default constructor. This, in turn, lets us create arrays of matrices. The first few additions to the class involve functions similar to those developed for class Vector , that let us examine and modify the class variables. These functions are:

    friend  size_t GetMatrixRows(void );    friend size_t GetMatrixCols(void );    friend void  SetMatrixRows(size_t m);    friend void  SetMatrixCols(size_t n);   

The usage should be pretty obvious.

The downside of using class variables is that they must be explicitly initialized–an initializer in the class definition won't do. What's more, the initialization must appear at the global level, in your main program. The form is rather specific:

    size_t Matrix::matClassRows = 3;    size_t Matrix::matClassCols = 3;   

Well, I suppose you don't have to initialize the variables. You could always leave them unset, as long as you remember to set them before the first matrix is constructed. However, as president and chairman of the board of the Society for the Elimination of Uninitialized Variables (SEUV), I don't think much of this idea.

Remember, these functions deal only with the dimensions of the default matrix–that is, a matrix created by the default constructor. Other matrices, created with explicit dimensions, or copied from somewhere else, are not affected by the class variables.

Once an instance of class Matrix is created, its dimensions are fixed. You can't change them, and changes to the class variables don't affect them either. You can, of course, find out what the dimensions are. I've included two forms of the query functions:

    size_t Rows(void );    size_t Cols(void );    friend  size_t Rows(const  Matrix & a);    friend  size_t Cols(const  Matrix & a);   

This approach lets us choose the usage we prefer:

    m = A.Rows( );   

or:

    m = Rows(A);   

Under construction
In my skeleton version of class Matrix , I only provided the default constructor. Any decent class needs at least a copy constructor, and we usually need others for constructing the object from other data types. In the case of matrices, the potential for different kinds of constructors is both rich and rather tricky, so I'll take the time to discuss the choices in depth.

The default constructor and the copy constructor are shown in Listing 1 .

View the full-size image

Note the use of the low-level functions vClear and vCopy , from the file SimpleVec.cpp . More about these functions later. The use of vClear is demanded by the charter of SEUV.

Array to matrix
Perhaps the most important and useful constructor is the third one, which builds a new vector, copying its data from an ordinary C array. Why is it so important? Because it's one of the few ways, short of reading an input file, that we can stuff data into the matrix elements. It's the closest we're likely to come to an initializer for objects of this class.

The constructor, shown in Listing 2 , borrows heavily from its cousin in class Vector .

View the full-size image

The difference is that, in the earlier constructor, the transfer of data is a simple element-by-element copy from one singly dimensioned array to another. In the new constructor, the transfer of data is a little more subtle. It also goes back to that issue of how array data are stored in memory, which I discussed in my last column. As before, it's worth taking the time to discuss the details. As we do so, we'll discover one last “gotcha” that the C syntax has in store for us.

Look carefully at the declaration of the constructor. The input array is declared to be double a[ ] , that is, an ordinary single-dimensioned array of doubles. To initialize that array, I might write:

    double a[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9};   

The declaration:

    Matrix A(a, 3, 3);   

would then construct a 3×3 matrix. In the process of building the matrix, we have transformed the shape of the data item from one-dimensional to two-dimensional.

Does this seem a little strange? Couldn't we, instead, define the array:

    double b[3][3] = {1, 2, 3, 4, 5, 6, 7, 8, 9};   

Now, you and I and the designers of C and the ISO/IEC C++ Standards Committee know that the two declarations have identical results. The C language doesn't store size or shape information anywhere in the executable software. Look at the two arrays with any debugger, and you'll see the exact same block of nine consecutive doubles.

What's more, we've all read the introductory C textbooks where we were told that any reference to an array is equivalent to a pointer to its first element. Thus:

    a     b    ptr_a = &a[0];   

and:

    ptr_b = &b[0][0];   

should all give the same result: a simple pointer to a double. Indeed, all through my code for class matrix, you'll see references like p[i] , where p is declared double* . It's a common practice and a fundamental part of the C idiom.

Despite all this obvious stuff, try passing a reference to the array b , and you'll find that some pointers are more equal than others. The declaration:

    Matrix B(b, 3, 3);   

doesn't work.

How can this be? If no size data are stored anywhere, the arrays a and b look identical in memory (and are perhaps in the same memory locations), and both references are passed as simple pointers, why then does the statement work differently for the two arrays?

The answer is simple enough: the compiler knows all. Although it doesn't find the need to store extra information anywhere, it knows that the two arrays are fundamentally different, syntactically. It won't let you treat them the same. To do so would be like equating two objects of different types.

What we'd really like to do is to define the constructor as:

    Matrix::Matrix(double a[ ][ ], size_t nr, size_t nc)   

but this doesn't work. The compiler has no way to tell the constructor the number of rows in a . Hence my earlier rant concerning C's lack of support for conformant arrays.

Of course, I can still use the converting constructor with a doubly dimensioned C array. I only have to make sure it's passed a pointer to its first element. The statement:

    Matrix B(&b[0][0], 3, 3);   

works just fine, albeit a little inscrutable.

Before we leave the topic, take one last look at the declaration:

    double b[3][3] = {1, 2, 3, 4, 5, 6, 7, 8, 9};   

I tend to use this kind of construct a lot, but strictly speaking, it's not good practice. I've declared the array to be doubly dimensioned, but used a singly dimensioned initializer. That it works at all is yet another indication of the on-again, off-again equivalence of singly and doubly dimensioned arrays. To be rigorously correct, I should have written:

    double b[3][3] = { {1, 2, 3}, {4, 5, 6}, {7, 8, 9}};   

To those who e-mailed me that the storage arrangement used by the compiler–row-major or column-major–shouldn't matter to the programmer, I can't resist pointing out that, in this declaration, it danged well does matter. The three inner brackets in that initializer refer to rows , not columns , and we had better not forget it.

What should we make of all this? Mostly, it's that trying to do useful matrix math in C/C++ is not fun. The C constructs, including esoterics like:

    *b++ = *a++;   

work great for singly dimensioned arrays, but not so great for multiply dimensioned ones. I've found it best to tell the compiler that all the arrays are singly dimensioned, and keep the truth to myself, hiding the details down inside functions or objects. That's why I wrote the class Matrix and its ancestors in the first place.

One word of caution concerning the converting constructor: there is not (nor can there be) any range checking on the size parameters. If you create a matrix from an array, you're expected to provide size parameters that don't violate the range of the input array. If you break this rule, you won't break the compiler or anything, but you will get garbage elements towards the end of the matrix.

Finally, I should mention that, because I'm so enamored with 3-vectors and 3×3 matrices, I sought to have the dimensions of this constructor default to 3×3, independently of the class variables. This would have let me write things like:

    Matrix A(a);   

without having to pass any size parameters. This, however, turned out to be a Terrible Idea. Because the call looks just like a call to the copy constructor, my compiler tried to typecast the array to an object of class Matrix .

I guess I could have overcome this obstacle with use of the verbatim keyword, but what if I had? Then I would have been able to write two declarations:

    Matrix A;    // default constructor    Matrix B(b); // convert from array   

that looked almost identical, but worked in fundamentally different ways. The default constructor gets its size parameters from the class variables, but the converting constructor gets them from the default values of the passed parameters. In all likelihood, I'd end up with matrices of different sizes. Not good.

So I gave up the notion of default values. If you want to create a matrix from an array, you have to give the dimensions explicitly, as in:

    Matrix B(b, 3, 3);   

In class Vector , I provided one more constructor, which was specific for vectors of length 3. It took the three scalar components (x , y , z ) of the vector as arguments. I'd like to have a similar constructor that would let me build a 3×3 matrix from its three columns, expressed as vectors. However, to build this constructor, I'd have to include the header files for both the vector and matrix classes, and cross-reference them. I've chosen to defer this step until next time. There's a change coming in the file layout that will make this step much easier.

Assignment statements
For this column, I've provided one assignment statement, which is shown in Listing 3 .

View the full-size image

As usual, the code contains a couple of subtleties that may not be immediately obvious. First, since the destination matrix already has its array allocated, we must check to make sure it's the right size. If not, we have to deallocate it and allocate a new one. In past versions, I've tended to retain the destination array if it's large enough, trading speed for a few wasted bytes. However, in this latest incarnation of the class, I've chosen to keep things clean by making sure the two arrays have the same size.

Note, however, that they don't necessarily have to have the same dimensions . If, say, one matrix is 3×4, and the other 2×6, that's okay. Both arrays have 12 elements, and the function will update the size parameters correctly.

One important point: when we implement assignment statements that involve objects that depend upon dynamic memory allocation, it's usually vitally important to make sure we're not assigning the object to itself. Otherwise, the function could end up deallocating its memory just before copying data from it. Not good. That way dragons be.

In this case, however, there's no problem. If we're foolish enough to assign an object to itself, the two sets of sizes will match, the if-test will fail, and the delete-new sequence will be skipped. We'll end up copying data to itself, but that's an acceptable performance hit.

In class Vector , I provided a second assignment statement that took, as its argument, a single scalar value. That value gets stuffed into every element of the array. It's not a very useful operation for vectors, because I can only think of a couple of cases where a vector with equal elements can be useful. The statement is handy, however, for the construct:

    v = 0.0;   

and that's why I supplied it.

A similar function makes even more sense for class Matrix . As with the vector class, a matrix with all its elements set to some scalar number has limited value. However, recall that there's an entity called the scalar matrix. It's a square matrix with equal values down the diagonal, and zeros elsewhere. The matrix:

(1)

has the property that the product Kx multiplies every element of x by the constant k . The scalar matrix is an often-used construct in matrix math. The statements:

    O = 0.0;   
or:    I = 1.0;   

have obvious value. Aside from that, a scalar matrix is often used as a starting point for more complicated computations.

Stream I/O
The finished class Matrix will definitely have this kind of assignment statement, but I'm deferring that until the next column. The last code I'll show you this month is for the two stream I/O functions. They are shown in Listing 4 .

View the full-size image

There's nothing fancy about these functions–they're about as minimal as one could provide. However, they do perform their intended duties. More elegant functions should probably be deferred to cases where we have specific uses, and therefore specific formats, in mind. For hints, you might check the I/O formats used in MATLAB.

I will call your attention to the indexing construct in operator <<( ) . It's the same as used in function Cell() . In truth, I might have been well-advised to use that function. I blushingly admit that I made an error in the indexing–an error that cost me an afternoon or so to find. I had used the column size parameter instead of the row parameter, incurring an error that only showed up when the matrix was non-square.

There could be no better way to make the point. My shame in admitting the error is well worth it if it drives home the point that the indexing math is best hidden well down inside software written for the purpose. To sprinkle such code through a main program is an invitation to disaster. It only takes one mis-written instance to make your life miserable.

Things to come
As usual, I'm uploading the code for the class Matrix and its test driver, to the web site. Check the test driver for examples of how to use the new functions.

I've included tests for the stream I/O. However, remember that my usual test approach is silent, in the sense that it gives no output unless something goes wrong. Therefore, in the default form of the test driver, the tests of the I/O functions are commented out.

You may have noticed that I've deferred a lot of items until the next column. There's a reason for that: it's because I'll be making three changes to the file structure that should make things easier in the future. I want to complete those changes before pressing on with more functions in the class.

First, note that, in several places, I've used calls to the functions vClear and vCopy . These functions are found in the file SimpleVec.cpp . That file was originally written to do low-level vector operations. However, they also work for matrices if the operations are to be done element-by-element. We only have to pass them an array size of m*n . For that matter, even the matrix add and subtract operations are done element-by-element, so I could use vAdd and vSub , as well.

I could, but I'm not going to. While I don't absolutely need to write new functions called mClear , mCopy , mAdd , mSub , and so forth, I also don't pay much of a penalty to do so. Most of the functions contain only one or two executable instructions, so it's not as though I'll be wasting lots of RAM space.

There's another thing. Although the low-level matrix functions I've mentioned will be identical to those for vectors, many others will be decidedly different. These include the generation of scalar matrices, the product algorithms, the inversion algorithm, and others.

Accordingly, I've decided to write another collection of functions parallelling those in SimpleVec.cpp . I'd call it SimpleMat.cpp , except for the next change.

With the introduction of SimpleMat.cpp and its header file, I'm up to eight files that must be included in every build. This is starting to get excessive. I've done things this way in the past and never had much of a problem, but I can't say as much for other users of my software. They tended to have problems collecting up all the needed files.

You might also remember that I didn't include, in this month's effort, any operations that involved both matrices and vectors. Since almost everything we do with vector/matrix math involves both types, we're going to have to make allowances for those operations, and we should do so very soon. In practice, this means adding references to vectors in class Matrix , and possibly the other way round as well.

To clean things up, I've decided to combine all the code–the classes Vector and Matrix , plus the low-level routines in SimpleVec/SimpleMat , into a single file, plus its corresponding header file. This will make the file a good bit larger than I'm comfortable with, but it won't use much code space because the individual functions are so tiny.

Look for those changes, coming soon to a magazine near you.

Jack Crenshaw is a systems engineer and the author of Math Toolkit for Real-Time Programming. He holds a PhD in physics from Auburn University. E-mail him at . For more information about Jack

Leave a Reply

This site uses Akismet to reduce spam. Learn how your comment data is processed.