 # Boost :

From: lums_at_[hidden]
Date: 2001-03-24 07:50:35

--- In boost_at_y..., hankel_o_fung_at_y... wrote:
> --- In boost_at_y..., "David Abrahams" <abrahams_at_m...> wrote:
> > > IMHO, it is natural to use both index styles in a single
problem,
> >
> > (!) How does /that/ situation come up?
> A typical scenario is that you want to index n+1 increasing time
> points as well as the n time intervals (or their lengths) bracketed
> by these points. In some problem domains (e.g., financial option
> pricing), it is natural to use 1-based indices for time periods
> and 0-based indices for the time points.

Our experience with developing (and using) MTL, interestingly enough,
is that once you start writing your algorithms in a generic
programming style, the fact that there is any base index at all
becomes surprisingly irrelevant. There are cases where one might
have to expose the internal representations of a matrix (say, in a
compressed sparse matrix) when dealing with external code (importing
or exporting the matrix data), but to a very large extent, you just
use iterators and the notion of an index having a value vanishes.

Here is an illustrative example.

template <class Matrix, class VecX, class VecZ>
void mult(const Matrix& A, const VecX& x, VecZ z)
{
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j;

for (i = A.begin(); i != A.end(); ++i)
for (j = (*i).begin(); j != (*i).end(); ++j)
z[row(j)] += *j * x[column(j)];
}

This a matrix-vector product written in a generic style. Note what
is and what is not going on. There is no assumption about the
storage format of the matrix -- this will work for dense or sparse,
row major or column major, triangular, banded, etc. There is also no
assumption of indexing. All that is required is that the vector
classes have a consistent domain/range relationship via the linear
transformation represented by the matrix. The vectors (at our
mathematical layer of abstraction) are essentially tuples and we use
some information from the linear operator to get to the element of
interest within the tuple. That we used the terms 'row' and 'column'
and somewhat a red herring.

Disclaimer: We are also somewhat mixing layers of abstraction in this
example. The iterator i for instance is iterating over a container
of containers. This really belongs more at the storage
representation layer. A more pedagogically correct implementation
might have been to first get a handle on the storage and then jump
into the loops:

for (i = A.storage.begin(); i != A.storage.end(); ++i)
for (j = (*i).begin(); j != (*i).end(); ++j)
z[row(j)] += *j * x[column(j)];

(Just a random thought)

In a more mature version of such a routine there are also concept
checks that one might want to add to assert, e.g., that the Matrix is
a transformation from the one Vector Space to the other (and, one
would also want checks that the vector quantities being passed in
were elements of a Vector Space).

> > > and a good array library should support both styles.
> > > A matrix library, however, can only and should at least support
> > > the 1-based index, because matrix indices are 1-based in math.
> >
> > Are you using some subtle semantic distinction between "array"
> > and "matrix" here?
> In terms of arrangement of quantities or symbols, a matrix is
> always an array.

Conceptually, yes. Since all finite dimensional vector spaces are
isomorphic to linear arrangements of numbers (tuples), all linear
operators mapping from one finite dimensional vector space to another
are isomorphic to a collection of numbers arranged in tabular
fashion -- "array" form. Representing a matrix in a computer program
as an actual array -- meaning there is space allocated for all
possible entries in such a tabular form -- is only one way of
representing a matrix (usually called "dense") and there are even
variations within that (row-major, column-major, e.g.).

cf also a previous exchange I had with Dave Abrahams on this list
about "matrix-free" linear operators (probably better-named "array-
free" operators or "storage-free" operators -- but "matrix-free" is
the name used in the literature).

> > > Yet, even if the index convention problem is fixed, we still
need
> > > to solve the row-major/column-major problem ...
> >
> > What problem is that? I thought row-major/column-major was just a
> > distinction describing underlying storage format, and thus an
> > implementation detail.
> Sorry for my poor presentation. I simply meant that one has to
> make a design decision on whether the i in A(i,j), for instance,
> means the i-th row, the i-th column, or both (depending on the
> storage format). This is not just an implementation detail because
> it concerns whether expressions like
> (A*u)(1) = A(1,1)*u(1) + ... + A(1,n)*u(n)
> are promised to be identities.

This example is not quite the same thing and is mixing two different
levels of concepts. To the extent that a matrix A is modeled as a
tabular arrangement of numbers A(i,j) is the element at row i and
column j. This is the mathematical convention. And, mathematically,
the expression you wrote would hold (modulo the normal disclaimers
about floating point arithmetic). That is, mathematically, the
relationship v(i) = sum(A(i,j)*u(j)) would hold when multiplying A*u.

However, as shown above, you would never want to write the matrix-
vector product as two nested loops over integers i and j. In between
our mathematical abstraction of a matrix -- with indexing convention
as above -- is a layer mapping from underlying storage (possibly an
array, possibly a sparse format of some kind, possibly something
else), to the mathematical.

The issue of row major versus column major comes up in this layer
when you have to decide how you actually want to stuff your elements
into storage. In the dense case this is an issue of mapping between
the A(i,j) abstraction of the matrix to some A_rep[i][j] or A_rep[j]
[i] or A_rep[i*lda+j] or A_rep[j*lda+i] array representation.

>
> > > I hope the Boost folks can draw a line between vectors, matrices
> > > and arrays in their matrix & array library.
> >
> > What sort of line did you have in mind?
> Uhh, I don't have one in mind, because I don't know how other people
> use arrays. At any rate, it is essential to identify different
> orthogonal concepts, and not to let the base array class (if any)
> assume too much. If I just want to interface my old legacy array
> architecture with the library functions (exchange two rows, slice,
> ... etc.), I don't want to define array element multiplication
> just because the library assumes that my array is a matrix.

Based on previous exchanges across this list, I think it is becoming
fairly certain that there will be a very careful distinction drawn
between matrix and array.

Regards,
Andrew Lumsdaine