[comp.lang.c] Matrices

dl@rocky.oswego.edu (Doug Lea) (11/02/88)

It is true that a matrix package is planned for libg++. It is
also true that the guiding principle will be to allow the
use of stuff from the `Numerical Recipes' text for most of
the `interesting' algorithms, although the goal will be something
that allows easy translation of the algorithms, not the code!

While putting together a Matrix class looks straightforward, anyone
who has contemplated how to put together a set of Matrix classes that
support (1) C-like efficiency where it counts (2) a wide range of
Matrix types and representations (like tridiagonal, and delay
Matrices), (3) use of sensible Matrix constructs, and (4) a
natural, yet object-oriented programming style will appreciate the
fact that I and several other people with whom I have talked,
corresponded, and worked have already thrown out dozens of designs. A
few other designs are being pursued. January is optimistic but
possible. If you are interested in joining the fun, email to
me directly.

-Doug



Doug Lea, Computer Science Dept., SUNY Oswego, Oswego, NY, 13126 (315)341-2367
email: dl@rocky.oswego.edu        or dl%rocky.oswego.edu@nisc.nyser.net
UUCP :...cornell!devvax!oswego!dl or ...rutgers!sunybcs!oswego!dl

casey@admin.cognet.ucla.edu (Casey Leedom) (11/07/88)

[Doug Lea talks about matrix library he and several others are working on.
In particular he says that coming up with a matrix representations that
are flexible and efficient is very difficult.]

  Perhaps Doug or someone else could help us with a problem of a somewhat
similar nature.  We're developing a C++ development environment for
programming the Connection Machine (CM) from Thinking Machines
Incorporated in Massachusetts.

  First off, any CM object always has two parts: the actual storage used
by the object on the CM, and various description information stored on
the front end (FE) machine (the CM is driven by a FE computer which feeds
the CM a stream of CM instructions).  The CM object storage is obvious.
The FE object description is minimally the address of the objects storage
on the CM.

  The CM supports integers of *any* bit length.  We would like to define
a C++ class to support this which doesn't suffer from stupid performance
hits.  Specifically, we'd like to be able to do something like the
following:

	CM_int a(8, 0), b(16, 0), c(32, 0);

where the first numbers are the sizes of the declared objects.  It's easy
to see that we can simply add a length to the address description set,
but this incurs a performance hit.  Most objects are fully defined at
compile time.  In particular, we know that ``a'' (above) is 8 bits long.
The storage to remember this at run time is totally useless.  But this
isn't much of a concern for us.  We're much more concerned with the run
time checking used to determine which particular CM instruction to use
to, say, add two CM_ints.

  Essentially what we're looking for are compile time attributes of
objects.  One obvious way of doing this is to simply invent a bunch of
subclasses.  E.g.:

	CM_int_8 a(0);
	CM_int_16 b(0);
	CM_int_32 c(0);

but this requires that we define all possible CM_int_* subclasses (of
which there are literally thousands possible).  So then we use a CPP
macro to create subclasses, but it starts getting ugly here.

  I'll end this with a concrete example to try to forestall
miscommunication.  The CM implements two signed integer addition
instructions.  One is a three address, single length (the two sources and
the destination are all the same length), the other is three addresses
and three lengths (both sources and the destination may all be different
lengths).

  We could always use the second, but we take a bad performance hit if
we do.  We could always check the lengths of the arguments to "+" at run
time to decide which to use, but this uses up time on the front end (and,
as has already been pointed out above, is almost always already known at
compile time).  We could depend on automatic coersion to larger integer
subclasses and then always use the first CM instruction, but that
eliminates the speed advantage of letting the CM do it's work properly.

  We'd like to be able to implement "+" something like the following:

	inline CM_int operator+(CM_int &a, CM_int &b)
	{
		if (sizeof(a) == sizeof(b))
			return(CM_int(sizeof(a), CM_s_add_1L(...)));
		else
			return(CM_int(max(sizeof(a), sizeof(b),
				CM_s_add_3L(...)));
	}

Where the ``sizeof''s are only used to give you an impression of what we
want (this is actually true of the entire outline above).  Obviously it
would be acceptable to generate stupid code and then let the optimizer
toss null branches, etc.

  The reason I say that Doug's matrix problem has some of the same
problems can be seen in the simple case of matrix multiplication.  E.g.:

	matrix a(M,N), b(N,P);

obviously it's going to be impossible to declare all possible matrix
sizes/shapes.  Therefore, we have to declare the size/shape of each
matrix at object declaration time.  Also, most of the time we know
everything about the size/shape of a matrix at compile time.

  It should be perfectly possible to detect whether ``a*b'' is legal at
compile time for most matrix objects.  While this could easily be done at
run time via descriptors, this suffers both a performance penalty *and*
loses the benefits of strict compile time type checking.  A particular
line of code may well compile correctly, but never be executed till some
critical point of a programs life time causing the program to fail.

  Oh well, I've dribbled on enough here.  We're probably just being
stupid.  Hopefully someone will point out an obvious solution.  Flames
for stupidity welcome.

Casey