Re: Re: [boost] ublas notation

From: Joerg Walter (jhr.walter_at_[hidden])
Date: 2002-07-16 01:55:39

----- Original Message -----
From: "Toon Knapen" <toon.knapen_at_[hidden]>
To: "boost" <boost_at_[hidden]>
Sent: Monday, July 15, 2002 2:58 PM
Subject: Fwd: Re: [boost] ublas notation

> On Saturday 13 July 2002 20:10, Joerg Walter wrote:
> > > Would it be hard to support operator() for ranges? If not, I suggest
to
> > > do it: it is a short, clear and powerful notation.
> >
> > See below.
> >
> > > > A.project(range(0,3), range(4,6)).assign( prod(
> > >
> > > > B(range(4,6),range(0,10)),
> > >
> > > > C(range(9,19), range(0,2))));
> > > >
> > > > but we then realized, that we could use free functions to achieve
the
> > >
> > > > same expressiveness :
> > >
> > > > project (A, range(0,3), range(4,6)).assign( prod(
> > > > B(range(4,6),range(0,10)), C(range(9,19), range(0,2))));
>
> Joerg, in the rhs this example still uses an operator() with 2 range
> arguments (identical to our example). I suppose you mean here :
> project (A, range(0,3), range(4,6)).assign( prod(
> project( B, range(4,6),range(0,10)), project(C,range(9,19), range(0,2))));

Oops, yes.

> > > > This is the recommended way currently. The rationale to eliminate
the
> > > > earlier member functions was to decouple vector/matrix and vector
> > > > proxies/matrix proxies. See
> > > >
> > > > http://aspn.activestate.com/ASPN/Mail/Message/1156090
> > > >
> > > > > Another suggestion to make life simpler, is to define range() or
> > > > > norange() to indicate the full scope of the range, for example
> > > > > A(range(0,2), norange()) is then equivalent to
> > > > > A(range(0,2),range(0,A.size2()))
> > > >
> > > > Yes, I've seen similar in the past. We'll consider it.
>
> Blitz++ has Range::all() which is pretty clear IMO.
> BTW blitz++ also supports `A(Range(2,2), Range::all());`

Ok.

> > > > > A row can be selected as A(2, norange()) and a column as
A(norange(),
> > > > > A vector_range of a row is then A(2, range(2,4)), etc.
> > > >
> > > > Recommended:
> > > >
> > > > project (row (A, 2), range (2, 4)) etc.
> > >
> > > This looks complicated.
> > > Do not forget that mathematicians are no computer scientists.
> > > Mathematical expressions should be represented in a way as natural as
> > > possible, meaning: natural to mathematician, not to computer
> > > scientists. Mathematicians do not like long expressions.
> > > The designers of Matlab and Fortran90 have thought a long time about
> > > matrix
> > > operations, and the solutions they came up with have been accepted
> > > widely. I recommend to stay as closely as possible to that notation.
> >
> > Agreed, but I'm neither a language designer nor a compiler writer ;-)
I've
> > looked into this issue again and found, that we're able to support the
> > following two notations (with a slight modification of the current
version)
> > under sufficient standard compliant compilers:
>
> I admire your effort to keep ublas compile-able with MSVC but we OTOH we
> should'nt cripple the interface of such an important library because of 1
> compiler (and IIRC VC7 will be able to handle all necessary features)

I'm already working on the limit. So either I've to drop MSVC 6.0 support or
I've to find some workable compromises. VC7 doesn't support PTS.

> > #include <boost/numeric/ublas/config.h>
> > #include <boost/numeric/ublas/vector.h>
> > #include <boost/numeric/ublas/matrix.h>
> > #include <boost/numeric/ublas/io.h>
> >
> > int main () {
> > numerics::vector<double> v (1);
> > v.clear ();
> > std::cout << project (v, numerics::range (0, 1)) << std::endl;
> > std::cout << v.project (numerics::range (0, 1)) << std::endl;
> > std::cout << project (v, numerics::slice (0, 1, 1)) << std::endl;
> > std::cout << v.project (numerics::slice (0, 1, 1)) << std::endl;
> > numerics::matrix<double> m (1, 1);
> > m.clear ();
> > std::cout << row (m, 0) << std::endl;
> > std::cout << m.row (0) << std::endl;
> > std::cout << column (m, 0) << std::endl;
> > std::cout << m.column (0) << std::endl;
> > std::cout << project (m, numerics::range (0, 1), numerics::range (0,
> > 1)) << std::endl;
> > std::cout << m.project (numerics::range (0, 1), numerics::range (0, 1))
> > << std::endl;
> > std::cout << project (m, numerics::slice (0, 1, 1), numerics::slice (0,
> > 1, 1)) << std::endl;
> > std::cout << m.project (numerics::slice (0, 1, 1), numerics::slice (0,
> > 1, 1)) << std::endl;
> > return 0;
> > }
> >
> > *without* reintroducing the tight coupling of vector/matrix and
> > vector/matrix proxy classes. W.r.t. overloading operator() for ranges
I've
> > currently the problem that this would either reintroduce the earlier
tight
> > coupling of vector/matrix and vector/matrix proxy classes or break MSVC
6.0
> > compatibility, as far as I see. So this should be a feature for a future
> > version of ublas IMO.
>
> The project makes it easier indeed, it's not necessary anymore to declare
the
> types of all intermediate temporaries with typedefs.
> Still I wonder how
> 'm.project (numerics::range (0, 1), numerics::range (0, 1))' does not
> introduce tight coupling whereas 'm(numerics::range (0, 1),
numerics::range
> (0, 1)))' would ?
> And if I understand you correctly, your motivation is that users that
include
> matrix.h should not always have to include matrix_pr.h too because also
> having to include matrix_pr.h would be cause longer compiler times, right
?

Nope. The compromise is to add the project() functions to
vector_expression<>/matrix_expression<> for standard compliant compilers
instead of adding them to *every* statically derived class, which is
neccessary for MSVC 6.0 (that was tight coupling). If I change the project()
functions in vector_expression<>/matrix_expression<> to operator() then,
name resolution doesn't find them anymore due to element subscription via
operator(). So to achieve your goal I'd have to move element subscription
via operator() to vector_expression<>/matrix_expression<> too, but only for
non MSVC 6.0 environments :-(

> Well, 1)I find that ublas compiles surprisingly fast, 2) I think run-time
> performance and expressiveness is more important than compilation time, 3)
> most LA algorithms will also use ranges anyway and thus will include
> matrix_pr.h anyway.
>
> > > Anyway, I think that the `project' function would be an important step
> > > forward, although my concerns about simple expressions remain.
> > > I just wanted to express my concerns when mathematicians are going to
use
> > > ublas. It is your free choice to consider suggestions. The question
you
> > > have
> > > to pose is the following: are mathematians going to give up Fortran90
> > > when ublas is available?
> >
> > I simply do not know enough about the chances of Fortran90 (although I
> > believe, that it doesn't support generic programming like C++). My
current
> > goal is to get a good interface and to enable the reuse of platform
> > specific tuned (procedural) kernels at the same time.
>
> Agreed. Our concern is currently mostly the 'good interface'. I think, as
a
> cs, that the interface is clear and practical. But after years of
discussion
> with mathematicians about comparisons of C++, fortran and matlab, C++ is
more
> difficult when writing complicated LA expressions. And when writing an
> algorithm, the complexity of the syntax of the expression is proportional
to
> the number of bugs (due to typo's etc). Additionally, all mathematicians
are
> familiar with matlab thus supporting (almost) the same syntax makes ublas
> much easier to use.
>
> e.g. here's a small piece of code that is used in the inner product of the
> Gram-Schmidt orthogonalization in the Arnoldi method for the solution of
> the quadratic eigenvalue problem in ublas notation, short-ublas notation
and
> matlab notation :
>
> ublas notation:
> ----------------
>    typedef numerics::matrix_range<BlockType> RangeBlockType;
>    typedef const numerics::matrix_range<const BlockType>
ConstRangeBlockType;
>    typedef const numerics::matrix_range<const MatrixType>
ConstRangeMatType;
>    typedef numerics::matrix_range<MatrixType> RangeMatType;
>
>    RangeMatType(Temp, range(0,last-start), range(0,V.size2())).assign (
>       prod(herm(RangeMatType(WW, range(0,WW.size1()), range(start,
last))),
> ConstRangeBlockType(V, range(WW.size1(),VV.size1()), range(0,V.size2())))
>       + prod( herm(ConstRangeMatType(Hessenberg, range(start, next),
> range(start, last))),
>               prod(herm(RangeMatType(WW, range(0,WW.size1()), range(start,
> next))), ConstRangeBlockType(V, range(0, WW.size1()), range(0,V.size2())))
>             ) );

Ouch ;-)

> Shorter notation (as being discussed currently):
>
> Temp(range(0, last-start), range::all() ) .assign (
>   prod( herm(WW(range::all(), range(start, last))),
> V(range(WW.size1(),VV.size1()), range::all()))
>   + prod( herm(Hessenberg(range(start, next), range(start, last))),
>   prod( herm(WW(range::all(), range(start, next))), V(range(0,
WW.size1()),
> range::all()))
> ) );

Ok.

> Matlab notation:
>
> Temp(1:last-start, :) =
>    WW(:, start:last)' * V(size(WW,1)+1:size(VV,1), :)
>   + Hessenberg(start:next, start:last)' *
>       ( WW(:, start:next)' * V(1:size(WW,1),:) )

Inaccessible.

Regards

Joerg