Re: sparse matrices

[Andy Adler]

On Mon, 12 Apr 2004, John W. Eaton wrote:

> I'd like to get sparse matrices into the core Octave sources.

This will make a lot of people happy.

> I've looked at the sparse directory in octave-forge, and I think it
> would need a lot of work before I would include it as part of Octave.
> One big problem I see is that it makes extensive use of malloc and
> free, but for something in the core of Octave, it would really be best
> to use new and delete, and then only in constructors/destructors, so
> that we can avoid memory leaks.

The octave-forge sparse uses malloc/free because the underlying
SuperLU implementation does. Of course, this could be modified.

> For starters, I would just like to be able to store data in a sparse
> format and have some basic operations (assignment, indexing,
> arithmetic ops, etc.).  It seems that the SuperLU package provides
> some code for these things, but again, I think the implementation
> would need a lot of work before it could be included as a part of
> Octave.

SuperLU provides very minimal support for the basic sparse operations
you mentionned. Most of the rest was written by me. To be honest,
its surprisingly hard to write basic algebra operations on sparse
matrices. Keeping track of indexes while debugging the operations
was a huge challenge; not to mention trying to develop efficient
implementations.

> Does anyone know of some good general-purpose code for handling sparse
> matrix storage?

As of a few years ago, there were very few of these, and none under
a GPL compatible licence. I hope that that has changed.

SuperLU is fast, but I think somewhat less stable than that used
by matlab. Also, it is missing some column re-ordering functions.

> If something like this does not exist, then it seems to me that
> writing our own code for it should not be too hard.  The first
> question is, what is the best storage scheme to use?  It seems that
> something like compressed sparse column format is widely used, but it
> does not allow one to add or remove elements efficiently.  It seems
> that there must be better ways to manage the storage internally while
> we might be adding or removing elements, then we could convert to a
> particular storage format when needed to call a solver.

The compressed column sparse (Harwell-Boeing format) or compressed
row sparse seem to be by far the most common schemes. Also, this
is what Matlab uses internally. However, you are correct that
this format makes it hard to add or remove elements.

> Comments or suggestions?

By (biased) opinion is that we should:
 a) look for a good, free, sparse toolkit
and if that doesn't work
 b) modify the octave-forge stuff as appropriate.

As I mentionned, I think debugging new sparse matrix code is
harder than refactoring a working package.

Andy