Matrix libraries for C and C++

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

Document options Print this page E-mail this page Congratulations! developerWorks wins Jolt Hall of Fame award Rate this page Help us improve this content

Level: Introductory

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

01 Jul 2002

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.

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.

                        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.

Back to top

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:

                        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.

                        #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.

Back to top

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:

                        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:

                        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 Resources).

Back to top

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:

                        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.

                        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:

                        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.

Back to top

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.

                        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
                        

                        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
                        

                        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
                        

Back to top

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.

Back to top

Resources

• 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.
  
  
  
  
• See The C++ Scalar, Vector, Matrix and Tensor Class Library Standard Page for proposals to define the matrix.
  
  
  
  
• Check out Harry Kuiper's page which contains more information on Cooperware Matrix.
  
  
  
  
• Read the Matrix and Quaternions FAQ, and find out what quaternions are.
  
  
  
  
• For more on quaternions and four-dimensional rotation, see Quaternion Powers. In adding random_fill to CwMtx, I employed drand48; for more information on this, see the man page, drand48.3.
  
  
  
  
• Dwnload Blitz, and check out the Blitz home page for more information.
  
  
  
  
• Take a look at an example of Blitz implementing matrices.
  
  
  
  
• To find out more about mailing lists, go to the Blitz mailing lists.
  
  
  
  
• 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).
  
  
  
  
• Improving the memory-system performance of sparse-matrix vector multiplication gives techniques for improving the speed of matrix multiplication by more than a factor of two on superscalar RISC processors.
  
  
  
  
• For more on improving sparse matrix algorithms, read Fast and effective algorithms for graph partitioning and sparse-matrix ordering.
  
  
  
  
• Download a trial version of the basic edition of WebSphere Studio Application Developer for Linux.
  
  
  
  
• Find more Linux articles and tutorials in the developerWorks Linux zone.
  

Back to top

About the author

		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.