From octave-maintainers-request at bevo dot che dot wisc dot edu  Fri Mar  6 02:42:02 1998
Subject: Sparse Matrices
From: "John W. Eaton" <jwe at bevo dot che dot wisc dot edu>
To: Andy Adler <adlera at njc dot org>
Cc: octave-maintainers at bevo dot che dot wisc dot edu
Date: Fri, 6 Mar 1998 02:42:15 -0600

On 12-Feb-1998, Andy Adler <adlera at njc dot org> wrote:

| I'm moving the discussion of Sparse Matrix support to here from
| help-octave. After several people suggested different sparse matrix
| libraries from which to build sparse matrix routines for octave,

OK.  Sorry for the delay in responding.  I am interested in having
Octave handle sparse matrices.

| I took a look at the SuperLU package
|  <ftp.cs.berkeley.edu/pub/src/lapack/SuperLU>
| 
| The package is fairly large
| -rw-rw-r--   1 andy     andy      1579843 Feb  2 16:57 superlu.tar.gz
| but contains demos and other extras.

The copy I got also had a 1.6MB executable in the EXAMPLE directory.
Just without that, it was less than 1MB tarred and compressed.

| The compressed source is
| -rw-rw-r--   1 andy     andy       347648 Feb 12 08:58 superlu_src.zip
| 
| It's written in C++,

Looks more like just C to me, but if it works, that's not a big
problem.

| and has code for sequential and multithreaded
| calculation (I suspect multithreading would not be easy to implement in
| octave). It contains code to
|   Factor:  [L,U,p]= superlu(A)  unit lower triangular L, upper triangular U,
|                                 and permutation vector p with A(p,:) = L*U.
|   Solve:   x = lusolve(A, b)    solution to the linear system A*x = b,
|   Estimate condition number:
|   Equilibrate:
|   Refine Solution:

OK.  I'm not an expert here, but this seems like a reasonable start.
What about other functions, like SVD?  I think Matlab offers that,
doesn't it?

| The license is:

Probably OK for use with Octave.

| The package includes *.mex files for linking to Matlab. I tried this on a
| Sparc10 under Solaris and ran the test routines which showed that SuperLU
| is 3.0 times faster than Matlab's sparse routines for large matrices
| (4720x4720, nnz(A)=32164).
| 
| QUESTION: Is this a good package to link to octave?
|           Does anyone have oppinions/better ideas?

I don't know what is the best code to use.  I hope that we can provide
a relatively generic internal interface for this sort of thing so that
it won't be difficult to switch later, if necessary.

| PROPOSAL: A first try link to octave.
| 
| Initially it's probably best to write *.m and *.oct files to do sparse
| math.  If these work well, jwe can probably be convinced to add a
| sparse matrix data format internally, as well as exending the operators
| to deal with sparse matrices.

Yes, though it might happen faster if someone else did the work.

| SPARSE MATRIX FORMAT: SuperLU uses the Harwell-Boeing or (compressed
| column) format. This is different to the Matlab internal sparse format.
| The H-B format should probably be used internally. 
| 
| Thus a sparse matrix S would be S.nrow  -> number of rows
|                                 S.ncol  -> number of cols
|                                 S.nnz   -> number of nonzero elements
|                                 S.nzval -> vector of nonzero elements
|                                 S.rowind-> vector index to rows
|                                 S.colptr-> vector pointer to rows

OK.

In Matlab, if I have an NxM sparse matrix S and I do something like

  S(i,j) = x;

with either i>N, or j>M, is S automatically resized (as happens
for dense matrices)?  Do range/vector indices work?  What about
assignment of an empty matrix to delete rows or columns?

| FUNCTIONS TO WRITE:
| 
| sparse -> *.m file based on Matlab's sparse

Please don't say `based on' if what you really mean is `compatible
with'.  To me, `based on' sounds like `derived from', which is no good
for Octave.

|   S= sparse(i,j,s,m,n,nzmax) uses the rows of [i,j,s] to generate an
|     m-by-n sparse matrix with space allocated for nzmax nonzeros.
|   S= sparse(A) generates sparse matrix S from full matrix A

OK.

| sp_add             -> corresponds to A+B
| sp_subtract        -> corresponds to A-B
| sp_dot_times       -> corresponds to A.*B
| sp_multiply        -> corresponds to A*B
| sp_dot_divide      -> corresponds to A./B
| sp_divide          -> corresponds to A/B

Do we really need the sp_* forms here?  Does Matlab provide those
functions?  Once you have the data type, it is fairly simple to make
the operator overloading work within Octave for either same-type or
mixed-type operations.

| full               -> like Matlab 
| 
| superlu            -> *.oct file from *.mex file with SuperLU
| lusolve            -> *.oct file from *.mex file with SuperLU

OK.

Any suggestions on what to do next?

Thanks,

jwe