Re: [boost] matrix and vector iterators

From: Andrew Lumsdaine (lums_at_[hidden])
Date: 2001-06-10 16:35:50

> > One thing that we need to decide on is the interface for traversing
> > vectors and matrices. I see a couple options:
> >
> > 1. MTL-style iterators
> >
> >   *i          // return a matrix/vector element
> >   i.index()   // return the index of the element pointed to (for
vectors)
> >   i.row()     // return the row index (for matrices)
> >   i.colummn() // return the column indexn (for matrices)
> >
> > 2. iterator over "matrix elements" or "vector elements"
> >   *i         // returns a matrix or vector element object.
> >   value(*i)  // return the element value
> >   row(*i)    // returns the row index
> >   column(*i) // return the column index
> >   index(*i)  // return the index (for vectors)
> 3.
>   *i  ??    // returns a matrix or vector element object.
>   value(i)  // return the element value
>   row(i)    // returns the row index
>   column(i) // return the column index
>   index(i)  // return the index (for vectors)

One question to keep in mind while discussing this is what are the semantic
differences between these approaches and what are just syntactic
differences.  What are the ramifications of the semantic differences in
terms of what it is we want to do with MTL?

There is also the issue of the multi-level iterators which I think has
already been alluded to.

> I think Joerg and Mathias are concerned with higher-dimensional matrices
> (although I am not), which makes "row" and "column" hard to justify.
> Personally, I find it very helpful to have "row" and "column".

I am not sure what a higher-dimensional matrix is?  I know what a
higher-dimensional array is.  You can probably guess what I am getting at
here.  It is important to define exactly what we mean by a matrix.  I think
there are several different layers, which may benefit from different
interfaces.

At one level, a matrix is a linear transformation from one
finite-dimensional vector space to another.  It provides an operator*() to
apply that transformation.  One might call the thing at this layer a 'linear
operator' rather than a matrix.

At the next level down, a matrix has doubly-indexed values.  How should
those be traversed?  Note that not all matrices need to provide this
interface in order to provide operator*().  For instance, the matrix-free
operators I mentioned in an email long ago can apply a linear transformation
via operator*() but do not store any values.  I think at this level, for
matrices that do have doubly-indexed values, it might be appropriate to have
some kind of interface that provides access to the values and to the
indices.  I am not sure I like row and column though, those are to weighted.
One could take a mathematical point of view and talk about the range index
and the domain index (i.e., the indices of the vector spaces that the
operator is mapping between -- this is really where the indices in the
matrix come from).  Matrix is probably the best name for the thing at this
layer.

At the next level down, the values of the matrix are stored in some fashion
(dense or sparse, e.g., row oriented or column oriented).  Here it might
make sense to think about iterators specific to the type of representation
being used.  In the dense case, the thing we are talking about here could be
called an array perhaps.

At yet another level down are issues of the storage, allocators, and so
forth.

> I very much like the free function interface, but there is always the LWG
> 225/229 issue to be aware of. I am growing more convinced that under the
> current language rules, having an additional tag argument which ties the
> function to a namespace's semantics is the only good answer... but that's
> another topic I guess.

I also prefer the free function interface and there seem to be some good
semantic (as opposed to syntactic) reasons to go that way.

> > Also, there is the question of whether the above iterator interfaces
> > should apply to both sparse and dense, or whether the dense iterators
> > would just use the normal STL style iterators which would leave the
> > indices out of the picture (for dense iterators they can be obtained in
> > other ways).
>
> I guess it depends how often you feel you can write an algorithm that
works
> well for both sparse and dense by using the index interface. My feeling is
> that it is very useful during library development, because you can often
use
> common code to get something that works, although you may want to
customize
> it later for sparse/dense... and it's not unreasonable to think that a
very
> good compiler will not make you pay for such an interface to dense
> iterators. On the other hand, it does make the iterator twice as big, so
> it's not unreasonable to imagine that you will pay something for it.

This is actually quite a compelling reason to keep sparse and dense unified,
I think.

> I think the best efficiency is available with two iterator interfaces: one
> which allows indexing, the other which doesn't. This would have benefits
> even for some sparse codes: when you have the result structure
> pre-calculated, you don't need to look at the indices.

Per the above, I think there is a reasonable place in a layered approach to
matrices where having two kinds of iterators might make sense.

It would be quite helpful in a discussion like this to have a few example
algorithms to chew on.  The usual suspects would be something like vector
addition, matrix-vector product, matrix-matrix product, and LU
factorization.