Matrix libraries for C and C++

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

This article presents some of the currently available options for open source C/C++ matrix libraries employable within a Linux environment. Particular libraries discussed are Meschach, which provides routines for operating on matrices and vectors for projects coded in C, the Cooperware Matrix (CwMtx) for C++ coding, and Blitz, which provides an n-dimensional array class for C++ with integral, floating, complex, and well-behaved, user-defined types. Andrew Blais, who has contributed several articles to developerWorks, is a researcher and writer for Gnosis, Inc., and does work in neural nets.

Share:

Andrew Blais (onlymice@attbi.com), Researcher and writer

Andrew L. Blais divides his time between home schooling his son, A. Van Blais, being a husband to K.J. Krawczyk-Blais, researching and writing for Gnosis Software, Inc., and teaching philosophy and religion for Anna Maria College. On occasion, he breathes. You can contact Andrew at onlymice@attbi.com.



01 July 2002

Also available in Japanese

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 Resources 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 Resources 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
PlatformCompilerOut-of-cacheIn-cache
HPC-160KAI C++100.2%97.5%
Pentium IIegcs98.4%79.6%
Cray T3EKAI C++95.7%98.1%
Origin 2000KAI 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 Resources).


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
Library2x23x34x55x56x6
Blitz0.400.480.620.750.91
CwMtx2.643.574.585.606.60
Meschach0.170.270.450.600.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
Library2x23x34x55x56x6
Blitz0.871.712.834.346.13
CwMtx3.675.878.5911.9315.61
Meschach0.420.801.321.862.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
Library2x23x34x55x56x6
Blitz1.312.624.506.859.71
CwMtx3.675.878.5911.9315.61
Meschach1.172.454.286.639.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 Resources.

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

Click to see code listing

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.

Resources

Comments

developerWorks: Sign in

Required fields are indicated with an asterisk (*).


Need an IBM ID?
Forgot your IBM ID?


Forgot your password?
Change your password

By clicking Submit, you agree to the developerWorks terms of use.

 


The first time you sign into developerWorks, a profile is created for you. Information in your profile (your name, country/region, and company name) is displayed to the public and will accompany any content you post, unless you opt to hide your company name. You may update your IBM account at any time.

All information submitted is secure.

Choose your display name



The first time you sign in to developerWorks, a profile is created for you, so you need to choose a display name. Your display name accompanies the content you post on developerWorks.

Please choose a display name between 3-31 characters. Your display name must be unique in the developerWorks community and should not be your email address for privacy reasons.

Required fields are indicated with an asterisk (*).

(Must be between 3 – 31 characters.)

By clicking Submit, you agree to the developerWorks terms of use.

 


All information submitted is secure.

Dig deeper into Linux on developerWorks


static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=1
Zone=Linux, Open source
ArticleID=11226
ArticleTitle=Matrix libraries for C and C++
publish-date=07012002