Fwd: Re: [boost] ublas notation

From: Toon Knapen (toon.knapen_at_[hidden])
Date: 2002-07-15 07:58:12

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

> > >
> > > 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());`

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

> #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 ?
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())))
            ) );

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()))
        ) );

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),:) )