Re: [ublas] [bindings] naming for sparse matrices

Subject: Re: [ublas] [bindings] naming for sparse matrices
From: Rutger ter Borg (rutger_at_[hidden])
Date: 2010-02-09 03:25:06

Thomas Klimpel wrote:

> The most common name for this array is "row_ptr" / "col_ptr". But a c/c++
> programmer will have an uneasy feeling about the "ptr". UBlas calls it
> "index1_data_", mtl4 calls it "starts", glas calls it "compressed_index_"
> and eigen2 calls it "m_outerIndex".

This is actually a nice showcase of the diversity of tackling this problem.

>>   begin_minor_index: I like this one, as you describe (or
>> begin_index_minor)?
> 
> How about "begin_index_minor" and "begin_compressed_index_major"?

Excellent. A bit longish, but it describes very well what it does and 
returns. Good.

> I'm not sure we should support the mapped stuff at all. I don't think that
> elementary matrix operations can be implemented efficiently with this
> storage format. So it would only be useful during matrix assembly. And I
> don't think that mtl4, glas or eigen2 offer the mapped stuff at all (but I
> could be wrong).

Ok, let's drop support for mapped stuff until the point where more than a 
handful users start yelling at us about where it is :-).

> I guess the tridiagonal format is a special case of a special banded
> format. This may be best realized by a "banded_view" over a normal matrix,
> like mtl4 offers. But I would have to study both ublas and mtl4 more in
> detail to see whether their offering for banded matrices supports this.
> Without support for tridiagonal matrices in the existing matrix libraries,
> the current bindings for tridiagonal matrices might be not so bad after
> all.

Not the most important case, yet, too?

> I tried to implement the sparse matrix stuff this weekend, but got buried
> in studying umfpack, ublas, glas, mtl4, eigen2 and the bindings themselves
> to learn how this would work. I didn't even had time to look into SuperLU
> and MUMPS in more detail. And I studied the tag-dispatching structure of
> the new traits system, and realized that the sparse stuff will also have
> to add new tags for its commands. The yet to be invented names for the new
> tags finally stopped me from doing any real work at all...

Let's get rid of the index<> tag first, it's too confusing with the sparse 
index stuff. What about renaming the index<> tag to dim<>? There are also 
some meta-funcs in index.hpp which should be renamed. I would add two static 
functions to compressed containers, 

begin_index1 and begin_index2,

and make the begin_index_major et al. dispatch according to the data_order 
(like done with size_major / size_minor). 

> And the coordinate_matrix format seems actually quite tricky to handle (in
> the case of umfpack), because the "sort" method must be executed, and the
> "compressed_index_major" array must be build and stored (and its lifetime
> must be "managed").

We can't use a trick like std::sort( begin_..., end... )? :-)  Sounds like a 
workspace thing, perhaps we can generalize workspaces somehow? Many FFT 
routines also use stuff like "plans", perhaps this falls into this category, 
too. 

> (Can we remove "has_rank.hpp" and "min_rank.hpp"? How do
> "begin_row"/"begin_column" actually work? I haven't understood how the row
> index is passed to the command.)

I was thinking of the usefulness of the rank function itself, I am seeing a 
couple of 

typename boost::enable_if< mpl::less< rank<T>, mpl::int_<2> >,
typename boost::enable_if< mpl::equal_to< rank<T>, mpl::int_<2> >,

which in the second case could be reduced to 

// implement the matrix-specific implementation
typename boost::enable_if< has_rank<T,2> > 

but indeed, perhaps this may be done better (independent of rank (which the 
bindings are trying to be), and more appropriate for the context of some 
algorithm). I have no problem with deleting them.

begin_row and begin_column aren't finished yet, it is something else I am 
thinking about: it can be an iterator over rows/columns (dereferencing the 
iterator will give a vector), or a shortcut for, e.g., begin( column( m, i ) 
), I guess the latter does what you mean. But I'm leaning towards

begin_column( m )             // vector-iterator
begin( column( m, i ) )       // same as begin_column( m, i )

why not begin_column( m, i )? Well, size_column( m, i ) also doesn't exist 
(if it's possible to keep begin/end/size consistent, would be nice, wouldn't 
it?). OTOH, begin_column( m, i ) adds to the user-friendliness, which it 
should be all about.

Erm... the dispatching may be a bit overwhelming at first. E.g., size1 goes 
like this:

size1(m) -> size( m, index<1> ) -> detail::size_impl< M, index<1> ->
 if index outside of rank of object ->
    return either 1 or 0 (as either ptrdiff_t, or mpl, this depends)
 else
    get< size_type<1> >( m ) ->
      if size_i is an mpl::int
         return compile time stuff
      else
         use PP trickery to translate index<1> to size1, call
         adaptor_access<T>::size1(m)

and, size_row is nothing else than 
size_row( m ) -> size( m, index<1> ) ....

etc.. With my head being fully into this, shall I try to glue in some of the 
sparse stuff?

Cheers,

Rutger