At 04:34 PM 9/28/98 +0100, Nick Leaton wrote:
>Should the matrix class understand eigen systems?
>I don't think that it should. Have a separate class that
>is called eigensystem, that takes a matrix as an input.
>Matrix inversion is similar. Then you can have a sequential
>eigen system, parallel eigensytem all inheriting from the
>abstract base class.
This kind of thing has come up over and over again in our work on OO for
numerics, and I believe Nick Leaton has the right general approach. If
your matrix class has member functions for eigensystems, inversions,
whatever, then extending the system will always involve modifications to
existing classes. Algorithms like these should be separate classes that
use primitive operations in the matrix class(es). This is an example of
the Open/Closed principle (open to extension, closed to modification) first
described by Bertand Meyer (I think) in the context of Eiffel.
The key thing to doing this sort of thing well is to give your array
class(es) the right set of primitive operations, and write your algorthms
using them.
For example, fortran arrays support reading and writing individual elements
as its primitive operations. F90 adds reading and writing array sections
to the primitive operations. For linear algebra, BLAS 1, 2 and 3 were
defined and then a great many algorithms written in terms of those.
Those primitive operations would then be implemented in specialized ways
for different types of matrices. Different higher level algorithms (like
eigenstuff, inverses and so on) would be written in terms of the primitive
operations.
What does this buy you? Suppose you have N types of matrix and M types of
high level algorithms. If you put the algorithms in the matrices, the
amount of work you have to do scales like N*M. If you separate them the
amount of work scales like N+M. It lets you add new matrices without
touching your algorithms. It lets you add algorithms without touching your
matrices. Different people can be doing different kinds of projects.
These are all Good Things.
There are different ways of implementing this sort of thing. One of the
big questions (if you're using C++) is whether you want to use compile-time
or run-time polymorphism.
Run-time polymorphism means using the techniques in C++ Classic: abstract
base classes, virtual functions, inheritance heirarchies, etc. You'd have
a heirarchy of matrices and possibly many heirarchies of algorithms.
Compile-time polymorphism means using templates, and is less familiar to
people since robust template implementations are a relatively recent
development.
A technique my team uses extensively is called "policy template
parameters". With this you define what I call an "abstract template class"
like
template<class StoragePolicy> Matrix {};
Then you define particular kinds of Matrix by defining tag classes for the
particular kinds of Matrix (like dense, upper triangular, banded, etc) like
so:
class Dense {};
class UpperTriangular {};
class Banded {};
Then you specialize Matrix for each of these:
template<> Matrix<Dense> { ... };
template<> Matrix<UpperTriangular> { ... };
template<> Matrix<Banded> { ... };
The implementations for each of these are then fully specific to that
style, and can be implemented with inline functions if necessary for
efficiency. You can define various kinds of parallel matrices the same
way. As long as you support the same interface among all of these things,
you can use them interchangably.
The algorithms can then be written leaving that policy parameter arbitrary:
template<class Policy>
void foo(Matrix<Policy>& m)
{
...
}
[For the case of a complex eigenstuff library you would want classes
instead of the free floating function here, but you get the idea.]
Then foo will work with whatever type of Matrix you have, and unlike
run-time polymorphism, it can be optimized using knowledge of the leaf
class instead of the root class.
Unfortunately, C++ doesn't provide a way to ensure that the interfaces of
the leaf classes are the same in the way that it ensures that with abstract
base classes. The package we're working on, POOMA II, ensures that through
a level of indirection we call Engines. I can describe that if people are
interested.
We'll be releasing a developer beta of POOMA II in November at the
Supercomputing conference in Orlando. The structures and patterns in that
library grew out of the same kinds of questions about how to structure
class heirarchies for doing scienific computation that Eric Noulard was
asking about.
-Steve Karmesin
======================================================================
Pooma Team Leader
karmesin@lanl.gov
phone: (505)665-6019
fax : (505)667-4939
For PGP public key, send message with subject line: SEND PUBLIC KEY
======================================================================
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 03:20:06 EST