Ron Burk has a B.S.E.E. from the University of Kansas and has been a programmer for the past 10 years. Helen Custer holds degrees in computer science, English and psychology from the University of Kansas and is currently researching a book on operating systems for Microsoft Corporation. You may contact them both at Burk Labs, P.O. Box 3082, Redmond, WA 98073-3082.
Math.h++, from Rogue Wave Associates, is a library of C++ classes designed for efficient numerical computing. Most of the code performs various numerical operations on vectors and matrices. The library also supports Fast Fourier Transforms (FFTs), linear algebra, statistics and random number generation. This report examines the Math.h++ library, v3.3.
Why C++?
Several well-established program libraries already exist for mathematics and statistics, such as NAG and IMSL. They are typically written in FORTRAN but can usually be called from a variety of programming languages. This numerical C++ library simplifies the task of writing large programs that are readable and maintainable.
A principle of good computer program design is "problem-solution similarity". If this principle is followed, the code should contain representations of some of the same components that make up the problem being solved. C++ doesn't have built-in support for lists as Lisp does, or for arrays as APL does, or for complex numbers as FORTRAN does. However, with C++ the programmer can extend the language to support whatever constructs are needed.
Consider the Mandelbrot set, which is characterized by the following iterative formula:
Z = Z * Z + CYou might expect to see this formula in the code of a program that plots the Mandelbrot set. However, since C and C++ do not directly support complex numbers, the Mandelbrot formula is usually obscured by the details of complex arithmetic.
It is possible, however, to add a complex number type to C++. In fact, complex numbers are one of the simpler classes supplied by the Math.h++ class library. Listing 1 shows a Zortech C++ program that uses the Math.h++ complex number class and the Zortech graphics package to draw the Mandelbrot set on a PC.
Adding complex numbers to C++ is only a minor feature of Math.h++; the library's real strengths are vectors and matrices. For example, Figure 1 shows an electrical circuit and the matrix equation that characterizes it. Finding the circuit's currents requires inserting actual component values into the formula and multiplying the voltage vector by the inverse of the impedance matrix.
If you write a program to perform this simple matrix inversion and multiplication, you can end up with a lot of code just for ancillary details (i.e., reading the data into the matrices). Listing 2 shows a simple solution using C++ and the Math.h++ library.
The Math.h++ classes are designed to work with the C++ "stream" I/O package (which uses the << operator for input and the >> operator for output). Listing 3 shows a particular input and output for the CIRCUIT. C program in Listing 2. The Math.h++ class library expects you to enclose matrices and vectors in square brackets ([]) and to specify matrix dimensions (the 3x3 in this case) before the matrix data.
This program took only a few minutes to write and required no debugging. The program is not concerned with how large the vectors and matrix will be, or with whether matrices are stored in column or row order, or with any other details not directly related to the task. Larger programs take more effort, of course, but programming at a higher level of abstraction can significantly reduce the time and expense of development, debugging and maintenance.
Overview Of Classes
The library contains vectors and matrix classes for the following data types: signed char, unsigned char, int, float, double and double complex. With some effort, you can also compile new classes for vectors or matrices of other types. (Genericity is discussed later in this article.)
Math.h++ overloads arithmetic operators for addition, subtraction, multiplication, and division to work with vectors and matrices. For example, if A and B are matrices, then the following expressions are legal:
A = A + B; A += B;Both expressions add each element of B to the corresponding element of A. It is an error if both matrices do not have the same number of rows and columns; (A discussion of error handling appears later in this article.)
A variety of library functions operate on the class objects in the library. The following standard C++ functions are overloaded so that they also operate on vectors and matrices:
abs() acos() asin() atan() atan2() ceil() cos() cosh() exp() floor() log() log10() pow() sin() sinh() sqrt() tan() tanh()For example, given a matrix A, the following code produces a matrix in which each element in B is the absolute value of the corresponding element in A.
DGEMatrix B = abs(A);The library supplies specialized functions that operate on vectors and matrices, including functions for forming vector sums and cumulative sums, complex conjugates of complex vectors or matrices, determinants of matrices, dot products of vectors, the spectral variance of a complex vector and the transpose of a matrix.
Math.h++ implements the LINPACK algorithms for solving simultaneous linear equations. Math.h++ provides separate classes that contain the LU decomposition (a matrix transformation which speeds repetitive calculations) of a matrix. In fact, the code in Listing 2 implicitly performs an LU decomposition on the impedance matrix before solving the equation. If the program were going to solve the same equation multiple times with different inputs (voltages), it would be much faster to generate the LU decomposition once and use that result explicitly in the call to solve().
Math.h++ also provides support for other areas of numerical analysis. One set of classes support FFTs. This package also supplies random number generators that produce the following distributions: binomial, exponential, gamma, Gaussian, Poisson and uniform. Also included are a class that handles linear regressions and a histogram class that helps compile and display the distributions.
Slices And Picks
Another concept integral to the Math.h++ library is a vector slice. A vector slice is a subset of vector elements formed by starting at a certain position and taking every nth element until the desired number of elements are collected. Numerical algorithms can often be coded in terms of vector slices. Such solutions are more easily ported to a vector processing machine.
In the Math.h++ library, vectors and slices are the same — a vector is just a slice that starts at the 0th element and has a step size of 1. Since a slice is not a special case, it can be used wherever a vector can, and vice versa.
You can refer to arbitrary subsets of vectors using vector picks. Picks are implemented with a separate class and are not as efficient as slices, because the elements selected must be remembered explicitly, rather than calculated.
To form matrix picks, you specify a set of rows (a vector of integers) and a set of columns (another vector of integers). The picked elements are the intersections of these rows and columns, which enables you to perform operations on selected portions of a matrix. For example, you could use a matrix pick to set the lower half of a matrix to 0.
Slices do not apply to matrices, but there are matrix operators that return the rows, columns and diagonals of a matrix. These operators return lvalues, a feature that supports coding techniques such as:
// 10x10 matrix DGEMatrix a(10,10); // set column 1 to zero a.col (1) = 0; // set main diag to -1 a.diagonal(0) = -1;Treating rows, columns, diagonals, slices and picks as atomic elements eliminates a lot of for loops from your code.
Class Design Issues
Every class library design begins by answering at least two important questions: Does the inheritance structure form a tree (all descended from one root) or a forest (many independent trees)? Do class objects assume reference or value semantics?
The Math.h++ library is a forest rather than a tree of classes. Although there are disadvantages to both approaches, merging a forest of classes into an existing class tree is probably easier than merging two class trees. For example, suppose you want all of your classes to be descendants of a class called Object. You could create a new class for each Math.h++ class that was simply the inheritance of both the Object class and the individual Math.h++ class.
One of the more confusing aspects of C++ class libraries is reference versus value semantics. The following statement, which assigns one matrix to another, illustrates the problem:
A = BAfter this statement is executed, does A contain a copy of the matrix B, or doesA simply refer to the same matrix as B?
The answer is a little complicated. In C++, assignment and initialization are treated differently. The Math.h++ class library generally uses reference semantics for initialization and value semantics for assignment. In other words, assignment creates a distinct copy of the object (vector, matrix, and so on) that is being assigned. With initialization, however, you can create two references to the same matrix. The library also supplies functions for forcing a distinct copy of an object to be made. Using Math.h++, you can almost, but not quite, ignore the issue for much of your coding.
Portability
In addition to the portability problems of different environments, a C++ program must also cope with variations among currently available C++ compilers. The language itself continues to change, so it is difficult to choose between all the latest language features or a program that is portable to the largest number of compiler implementations.
With normal portability problems, Math.h++ does a fair job. The worst flaw is some gratuitous name-space pollution. For example, the header files define the symbols TRUE and FALSE as macros, which will cause a syntax error when used in programs that define these symbols as enums.
There are also some unnecessary dependencies on non-ANSI I/O functions for doing binary file I/O. On the other hand, conditional compilation is used in many places to adapt to host environment differences.
The authors seem to have worked hard to accommodate differences between C++ compilers. The code has been ported to compilers from AT&T, Borland, Zortech, Oregon Software, HCR and GNU. Conditional compilation allows the code to adapt to various compilers. For example, for compilers that supply a complex number class, the library generally uses that class rather than its own, to make it easier to use the library with existing code.
The fact that no standard class library for C++ is defined poses a significant problem for writing portable C++ class libraries. For example, AT&T ships with a revised and enhanced stream class for doing I/O, but some compilers do not. Math.h++ handles this by using whichever stream class the compiler provides.
Efficiency
In general, the design of the library makes efficiency a high priority. Since reference semantics are used for initialization, matrix copying is avoided in places where the compiler might generate a temporary matrix variable. Likewise, since slices are an integral part of the vector design, there is no significant penalty for using a slice instead of a vector.
The space efficiency of the library depends on how smart your linker is. The member functions are grouped into relatively few object modules, so many linkers will load the entire object module even if your program uses only one of the functions. The CIRCUIT.C program in Listing 2 results in a 55K executable under MS-DOS.
If space is critical, you could recompile the library in a more granular form, but this is tedious without an automated tool. Perhaps linkers will become smarter as object-oriented languages proliferate.
On a related subject, matrices are stored in column-order, as in FORTRAN. Since most numerical libraries are written in FORTRAN, this increases the odds that you will be able to mix the Math.h++ library and some other numerical library in the same program.
Genericity
One important feature that C++ lacks is genericity — the ability to define a single class for multiple data types. For example, the only difference between a float matrix class and a double matrix class is probably that each occurrence of the key word float is replaced by the key word double. Currently, the only portable way to accomplish this is with large, multi-line preprocessor macros that are difficult to read, use and debug.
The Math.h++ solves this problem by supplying its own program (a shell script, for UNIX) for expanding class "templates" into class definitions. These templates take a type name as a parameter and expand into the correct definitions for the given type. The supplied program is not set up to be used on new data types, but the basic task is simple text substitution, which can be done in an editor or with a simple sed script in UNIX.
Handling Errors
C++ detects many kinds of type errors, such as applying an operator to the wrong variable type. The Math.h++ class library can also detect several kinds of semantic errors, such as an attempt to access an element outside the bounds of a vector or matrix.
The library supplies bounds checking in a couple of forms. First, the subscript operator [] always performs bounds checking. Second, the function call operator () is overloaded to work just like the subscript operator, except it performs bounds checking only if the macro BOUNDS_CHECK is defined. The function call operator enables you to perform bounds checking during development and easily remove it to speed up the final product. Other kinds of errors, such as attempts to multiply matrices of incompatible sizes, are also detected. When the library encounters a runtime error, it prints a message on the screen and (depending upon the severity of the error) terminates the program.
Unfortunately, you can't easily replace this error function with your own. For example, if you use the library in a windowing environment, you may have to alter the error function in the library so that it displays errors in a window and terminates the program by posting the correct system event.
Summary
Math.h++ is not aimed at the same market as a down-and-dirty, wired-for-speed numerical library. Its primary strength is the level of abstraction it provides and the ease it provides in constructing correct, readable, maintainable programs. Its secondary strength is portability, although that is a moving target as C++ becomes more standardized. Also, since the library is supplied in source form, you can change or improve it, but you must maintain your changes through future releases of the library.
This is one of the earliest commercial C++ class libraries. Class library vendors face many problems, such as compiler incompatibilities and lack of support for genericity. Despite these problems, Math.h++ demonstrates that the promise of reusable class libraries is not all hype — the future looks bright.
Math.h++ Library
Rogue Wave Associates
P.O. Box 2328
Corvallis, OR 97339
(503) 745-5908
