# Matrix libraries for C and C++

An evaluation and comparison of Meschach, Cooperware matrix, and Blitz

The following article assumes some familiarity with C/C++ and a preoccupation with the fact that C/C++ *per se* lacks matrix functionality. You may be analyzing econometric data or modeling rain forests. As for me, I work with neural nets, the implementation of which is greatly simplified by a matrix or two. Although C/C++ includes containers that can be considered matrices (such as arrays and, in the Standard Library, vectors, lists, and maps), a container that actually *is* a matrix would make the tasks at hand far easier. So we're going to look at three open source options that don't require you to build your matrices from scratch, but do allow you to tweak your matrix library. This is especially good if you anticipate needing to tweak it in as yet unforeseen ways.

## Meschach: a C option

For projects coded in C, Meschach (pronounced: me-shark) provides routines for operating on matrices and vectors. It has the virtue of compiling under Linux and most other operating systems, and is openly available under copyright, provided the customary acknowledgment is observed and errors are reported. Meschach was designed to solve systems of dense or sparse linear equations, compute eigenvalues and eigenvectors, and solve least squares problems, among other things. It provides nearly 400 functions for doubles and complex numbers. It also comes with a tutorial that introduces these functions in the context of illustrative mini case studies. David Stewart and Zbigniew Leyk introduce Meschach through a discussion of topics such as a generalized least square equation solver for over-determined equations, and problems involving sparse matrices. Their tutorial also includes the slightly more advanced topics of 3-dimensional matrices and error reporting.

Objects and class functions are often associated with code, and C structures may seem somewhat arcane, so C libraries are often dismissed as viable solutions. But in its defense, this library is very well organized and should not be dismissed by default. A quarter hour after downloading Meschach,
I was making, filling, and displaying matrices (the moral equivalent of creating a *Hello World!* program). An affordable paper manual, "Meschach: Matrix Computations in C," is available for reference (see Related topics later in this article). In particular, the test program "torture" contains a number of helpful clues.

Matrices can easily be sent to files or to the standard output. Meschach computes Fast Fourier Transforms, extracts columns and rows, and computes eigenvalues of symmetric matrices. You can fill a matrix with random ints and complexes. The library, believe it or not, even has facilities for adding matrices. One of the Meschach features that simplifies writing the obvious stepping stone programs is a function that returns a random double in [0,1). Although Meschach has a function that fills a matrix with 1.0s, it unfortunately lacks a function for filling a matrix with an arbitrary double or with random doubles. But these are easy to add.

##### Listing 1. Filling a matrix with an arbitrary double or with random doubles

MAT *m_fill( MAT *A, double x) /* MAT *m_random_fill( MAT *A ) */ { int i, j; for ( i = 0; i < A->m; i++ ) for ( j = 0; j < A->n; j++ ) { A->me[i][j] = x ; } /* { A->me[i][j] = m_random ; } */ return A; }

So, as we've seen, Meschach code is easily extensible, although in an ideal world more thoroughly commented code is preferable. But if you happen to be doing computational work in C that needs matrices, this is a highly useful library.

## Cooperware matrix: the basics

If you want to write object-oriented code in C++ and you prefer conceptual clarity to speed, Harry Kuiper's Cooperware Matrix (CwMtx) will work very well. Of the three matrix libraries discussed here, I found its conceptual architecture to be the most intuitive. In constructing a matrix, you use the obvious:

##### Listing 2. Constructing a matrix

CWMatrix A( rows, columns ) ;

Of the three libraries considered here, CwMtx had the worst performance in terms of three evaluative tasks, which we'll look at in detail in the section on Speed. But when clarity is more important than performance (when, for example you want to be sure that your data is being crunched properly), CwMtx is a great option. Make it right first, then make it fast.

The matrices in CwMtx include vectors and square matrices, where vectors include space vectors and quaternions. A matrix can be mapped into a matrix, filled with some element, transposed, and subjected to the usual mathematical operations. Kuiper originally used CwMtx to simulate systems built from discrete and interactive state machines. In addition to the obligatory matrix class, a quaternion class also exists. The answer to the obvious question is that q is a quaternion just when q = r + xi + yi + zi, where r is a real number, i is the square root of -1, and x, y, and z are complex numbers. Quaternions make it possible to extend the concept of rotation in three dimensions to four dimensions (see Related topics for links to quaternion references).

CwMtx has no built-in random number generator, and there is no class function that fills a matrix with random elements. But because it's free and released under the GNU LGPL license, you have the freedom to create these if you'd like. As far as filling a matrix with random elements, the following option is a good and easy way to go.

##### Listing 3. Filling a matrix with random elements

#include <stdlib.h> ... void random_fill( CWMatrix &M ) { int SIZE = M.GetRows() ; for ( int r=0; r<SIZE; r++ ) for ( int c=0; c<SIZE; c++ ) { M[r][c]= drand48(); } }

The documentation here is slight, but it is clear and well organized. It is easy to find such things as the class hierarchy and the constructor and member function options. Although there is no tutorial, you won't miss it; between the documentation and the test program, it won't be needed.

## Blitz: as fast as Fortran?

Blitz is another C++ library distributed under the GNU GPL, and with it you can freely create objects. It supports the KAI, Intel, gcc, Metroworks, and Cray 3.0.0.0 C++ compilers and provides an n-dimensional array class that can contain integral, floating, complex, and well-behaved user-defined types. Its constructor is more complex that the CwMtx constructor, as is evident in this example:

##### Listing 4. The Blitz constructor

Array<double,2> A(4,7) ;

This creates a 4x7 rank 2 array that contains doubles. But because it's a bit unclear, Blitz makes you think of an array as a matrix. Also, it does not implement many matrix functions. For example, no function returns the eigenvalues of a matrix. Also, no function fills an array or matrix with random doubles. But, Blitz does have two basic virtues.

One virtue lies in its breadth. By employing its native capabilities, it is easy to implement and construct a random double filling function, as you can see from this example:

##### Listing 5. A random double filling function in Blitz

template <class Array, class Uniform> void fill( Array &a, int size, Uniform &u ) { for ( int i=0; i<size; i++ ) for ( int j=0; j<size; j++ ) { a(i,j) = u.random() ; } }

Blitz's Uniform class provides a member function that returns a double in [0,1). It also provides three methods for accessing array elements: standard indexing, creating subarrays, and slicing, which produces a less dimensional view of an array segment. Blitz also possesses the standard calculator type functions, so arrays can be displayed on the standard output, and can be read from and sent to files. Laplacian, gradient, and Jacobian operators are just three examples of Blitz's stencil functions.

The other virtue lies in Blitz's speed. Depending on which compiler is employed, C++ can deliver performances at or about that of Fortran, which is renowned for its high performance in scientific and engineering computing. Take a look at the comparisons in the table below, but see the section on Speed below for an analysis of this data and the performance claims based upon it.

##### Table 1. Blitz performance on different platforms

Platform | Compiler | Out-of-cache | In-cache |
---|---|---|---|

HPC-160 | KAI C++ | 100.2% | 97.5% |

Pentium II | egcs | 98.4% | 79.6% |

Cray T3E | KAI C++ | 95.7% | 98.1% |

Origin 2000 | KAI C++ | 88.1% | 79.8% |

Blitz comes with a manual in HTML and Postscript, but unfortunately, there is no tutorial. There is, however, quite a bit of illustrative code from which the nuances of Blitz syntax can be gleaned. And the class reference provides the usual information. There are also several helpful mailing lists that are archived and searchable (see Related topics).

## Speed

Libraries can be evaluated in terms of their functional resources, documentation, the quality of their tutorials, how easily they can be extended, and so on. They can also be evaluated in terms of performance and/or speed. But comparing them is sometimes difficult, because (as is the case with our examples) they aren't all written in the same language, and they don't have much native functionality in common. In our case, there is sufficient overlap among the libraries to permit comparison in terms of speed based on three relatively simple yet revealing tasks, which are displayed and discussed in the following three examples:

##### Listing 6. Task 1

For ( d=2; d<7; d++) Construct 3 dxd matrices: A, B, C Start Clock: Do the following one million times: Fill A and B with 1.0s Let C = A + B Stop clock: Report elapsed time in seconds.

Implementing and executing this algorithm using our libraries produced the following results:

##### Table 2. Results of task 1

Library | 2x2 | 3x3 | 4x5 | 5x5 | 6x6 |
---|---|---|---|---|---|

Blitz | 0.40 | 0.48 | 0.62 | 0.75 | 0.91 |

CwMtx | 2.64 | 3.57 | 4.58 | 5.60 | 6.60 |

Meschach | 0.17 | 0.27 | 0.45 | 0.60 | 0.79 |

Given its intuitive architecture, it is unfortunate that CwMtx doesn't perform better here. And Although Blitz is outperformed by Meschach, it is remarkable that it so vastly out performs CwMtx, its object oriented competitor.

Recall that both Meschach and Blitz have functions (Random Number Generators) that deliver random doubles, and that CwMtx lacks native random generating capacity. Given the pivotal role of randomization in some matrix-based simulations, it is instructive to investigate how these libraries perform in situations that call on randomization.

##### Listing 7. Task 2

For ( d=2; d<7; d++) Construct 3 dxd matrices: A, B, C Start Clock: Do the following one million times: Fill A and B with random doubles (using library RNG, if any) Let C = A + B Stop clock: Report elapsed time in seconds.

Our libraries performed as follows:

##### Table 3. Results of task 2

Library | 2x2 | 3x3 | 4x5 | 5x5 | 6x6 |
---|---|---|---|---|---|

Blitz | 0.87 | 1.71 | 2.83 | 4.34 | 6.13 |

CwMtx | 3.67 | 5.87 | 8.59 | 11.93 | 15.61 |

Meschach | 0.42 | 0.80 | 1.32 | 1.86 | 2.52 |

Once again, Blitz is outperformed by Meschach, but outperforms its object oriented competitor CwMtx by an astounding margin. Lest you think that this is due to the performance of RNGs, let's take a look at this third evaluative task:

##### Listing 8. Task 3

For ( d=2; d<7; d++) Construct 3 dxd matrices: A, B, C Start Clock: Do the following one million times: Fill A and B with random doubles (using shared RNG) Let C = A + B Stop clock: Report elapsed time in seconds.

Our libraries retain the performance rankings of the previous two tasks:

##### Table 4. Results of task 3

Library | 2x2 | 3x3 | 4x5 | 5x5 | 6x6 |
---|---|---|---|---|---|

Blitz | 1.31 | 2.62 | 4.50 | 6.85 | 9.71 |

CwMtx | 3.67 | 5.87 | 8.59 | 11.93 | 15.61 |

Meschach | 1.17 | 2.45 | 4.28 | 6.63 | 9.40 |

As you would expect, CwMtx doesn't change its rank in performance capability. Moreover, Blitz and Meschach retain their relative rankings. If raw speed is the deciding factor, the order of rank among these libraries is now clear.

## Installing and compiling on Red Hat Linux 7.1

For convenience, I include the following notes on installing and compiling with these three libraries. Links to the downloads can be found in Related topics.

##### Listing 9. Blitz

tar zxf blitz-0.5beta3.tar.gz cd blitz-0.5beta3 ./configure --with-cxx=gcc make all cp -a blitz/ /usr/local/include/ ( or whatever you wish ) cp -a lib/ /usr/local/include/blitz/ Compile with: g++ -O2 pgm.cpp -o pgm

##### Listing 10. Meschach

unzip -q mesch12b.zip -d mesch12b cd mesch12b ./configure make basic mkdir /usr/local/include/meschach ( or whatever you wish ) cp *.h meschach.a /usr/local/include/meschach/ Compile with: gcc -O2 pgm.c /usr/local/include/meschach/meschach.a -o pgm

##### Listing 11. CwMtx

tar zxf cwmtx-0.3.0.tar.gz cd cwmtx-0.3.0 Open Makefile, and set INSTALL_DIR to /usr/local/include/cwmtx ( or whatever you wish )Execute: make install Compile with: g++ -O2 pgm.cpp -o pgm

## Summary

We've now looked at the features of three matrix libraries, and detailed their native capacities. We've also seen a number of functional shortcomings among them, and have looked at ways that these might be overcome. I have provided you with some simple tests that provide some crude quantitative data that might aid you in your coding choices, but in the end the decision is yours and should depend on the individual aspects of your project as well as the speed of the libraries in the given circumstances.

#### Downloadable resources

#### Related topics

- Download Meschach and get information on Meschach at David Stewart's site, or send e-mail to netlib@research.att.com with "send all from c/meschach" in the body (without the quotes). The site gives information on how to buy the Meschach manual, "Meschach: Matrix Computations in C". The random number generator is based on Knuth's lagged Fibonacci-based generator. See "Seminumerical Algorithms: The Art of Computer Programming" sections 3.2-3.3.
- Check out Harry Kuiper's page which contains more information on Cooperware Matrix.
- To learn more about the random number generator implemented in Blitz, read "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp. 3-30.
- Check out Matrix TCL Lite (free) or Matrix TCL Pro (peddled). The former worked with my compiler of choice, g++ 2.96, but since the library itself isn't optimized, it wouldn't fit in with my analysis. The latter is optimized, but since it wouldn't work with my compiler, it also wouldn't fit into my analysis.
- Get a list of object oriented libraries for scientific computing.
- Read Andrew's article "An introduction to neural networks" (
*developerWorks*, July 2001). - Download a trial version of the basic edition of WebSphere Studio Application Developer for Linux.