For the past year I have been working on an Object-Oriented framework for
a special class of large-scale optimization algorithms implemented in C++
called rSQP++. I spent a lot of time looking from a good fondation for the
numerical computations but I was not satistied with anything I found. My
application domain must deal with both sparse and dense linear algebra
structures. As far as I can tell most of the work with OO numerical
libraries has been with dense matrices (2-D and larger) and vectors. When
it comes to dense vectors and matrices (2-D) it should not be too
difficult to standarize an interface to these data types. What will be
difficult is to define standards for sparse matrices and vectors.
In my application domain I have to call several Fortran numerical
subroutines (e.g. MA28) so compatability with them (which dominate the
computing time) is one of the
most important design constraints.
Because of all of this I have implemented my own dense linear algebra
class library called LinAlgPack and it is designed primarily for BLAS and
LAPACK compatability and uses STL like interfaces for VectorSlice objects
which can represent vectors, regions of vectors or rows, columns or
diagonals of 2-D matrices. These are compact abstract data types with
very little overhead and a lot of inlined funcitons.
I also have implemented several sparse matrix classes that primarily
support the Fortran compatable coordinate data structure
(val,ivect,jvect). These classes facilitate the permutation and
partitioning of these sparse matrices and interfaces to current Fortran
sparse linear solvers (MA28, MA48). These are also impemented as
concrete abstract data types.
Dealing with concrete data types at the application level however turned
out to be a nightmare so I was forced to come up with a set of abstract
interfaces for the matrices and their operations.
Now, here is what I think about standards for numerical computation in
C++:
1) Until most of the major compiler venders come up with stable
implementations of the C++ standard (e.g. templates), programmers like me
will not be able to experiment with all of the techniques that are
possible. And until the implementations of these techniques is stable and
portable (templates again are the main problem) no one will trust them
farther than you can
throw a 19'' monitor (including myself and I am pretty strong).
Therefore, not even dense linear
algebra standards are possible at this point.
2) In may real-world numerical applications like optimization and
simulation (both
important in Chemical Process Systems Engineering) the handling of sparse
linear algebra objects is a major concern and tend to dominate the run
time. In my limited experience, the
requirements for sparse matrix and vector manipulation are as diverse as
the application domains they are used in. Therefore, I don't think it
will ever be possible to totally standardize spare linear algebra. For
example, there are spare matrix classes in the TNT library
(http://math.nist.gov/tnt/) but a few minor details made them
inappropriate for my application. However, I do think that some common
types of data structures and operations should be standarized but done so
with some specific compatability constraints (e.g. Fortran). These would
be less than optimal but at least they would be avalible for quick use.
Once my application framework is more stable I will publish a paper on it
and make some of its building blocks publicly avalible.
--------------------------------------
| Roscoe Bartlett |
| Ph.D. student |
| Department of Chemical Engineering |
| Carnegie Mellon University |
| email: roscoe@andrew.cmu.edu |
--------------------------------------