Re: [boost] matrix and vector iterators

From: David Abrahams (abrahams_at_[hidden])
Date: 2001-06-10 18:27:28

----- Original Message -----
From: "Andrew Lumsdaine" <lums_at_[hidden]>
To: "David Abrahams" <abrahams_at_[hidden]>; <boost_at_[hidden]>

> 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.

Yes, but I confess that I never really understood why it would be important
to have such a thing (as opposed, say, to different views of the matrix and
the ability to transform one view's iterator into another).

> > 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.

Given the name "matrix-free" I'd rather call this "LinearTransformation"
than "Matrix".

Of course, you /can/ still provide double-indexing for matrix-free operators
by operating on unit vectors, though it might not be efficient. In my view
it will be a mistake to require lvalues for matrix element access, so that
part shouldn't be a problem.

> 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).

There's a lot of value in easily understood names, and algorithm writers
will commonly be thinking in terms of row and column. One nice thing about
free functions is that you can provide both interfaces without too much
clutter or intrusion.

> Matrix is probably the best name for the thing at this
> layer.

OK.

> > 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 agree, at least at some level.
>
> > 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.

Sure, but that sort of begs the questions, doesn't it? We're talking about
the terms in which these will be expressed.

-Dave