[Blitz-support] Is this a bug, or am I missing something?

[Fernando P. <fp...@co...>]

Hi all,

the strange behavior I was seeing with the old Blitz code is still there, so 
I'd like to know if I am misusing the library in some way.

I am trying to compute the product of a rank-2 array M with a rank-d array T, 
by summing over the 2nd index of M and the last one of T:

U(i1,i2,...,id) = sum(M(i1,j)*T(i2,i3,...,id,j),j),

or in tensor notation

U             = M      T
  i1,i2,...,id     i1,j  i2,i3,...,j

I had initially written code to do this by hand-rolling nested loops with 
preprocessor macros (I need d=1..6), but this is super-ugly.  I later read the 
documentation, and found that Blitz supports essentially the above notation 
unchanged.  So I wrote a new version using this, but it doesn't work correctly 
for d>1 (d=1 is ok).

I'd like to know if I am misusing the library (case in which perhaps the docs 
need some clarification, b/c I read them fairly carefully many times), or if 
this is indeed a problem.

The relevant part of my code is (showing d=1,2,3):

// 1d version
void mat_ten_inner0(Array<double,2>& M,
                    Array<double,1>& T, Array<double,1>& U ) {
     firstIndex i0;
     secondIndex j;

     U = sum(M(i0,j)*T(j),j);
}

// 2d version
void mat_ten_inner0(Array<double,2>& M,
                    Array<double,2>& T, Array<double,2>& U ) {
     firstIndex i0;
     secondIndex j;

     secondIndex i1;
     // We want for(i0,i1) U(i0,i1) = sum(...,j): j is the sum index
     U = sum(M(i0,j)*T(i1,j),j);
}

// 3d version
void mat_ten_inner0(Array<double,2>& M,
                    Array<double,3>& T, Array<double,3>& U ) {
     firstIndex i0;
     secondIndex j;

     secondIndex i1;
     thirdIndex i2;

     // We want for(i0,i1,i2) U(i0,i1,i2) = sum(...)
     U = sum(M(i0,j)*T(i1,i2,j),j);
}

I can post a fully buildable test case if anyone wants to run it without 
having to write anything themselves.

I made a debug build, and here are the results (a non-debug build gives 
incorrect results for d=2 and segfaults for d=3):

planck[bug]> ./mtinner
Rank 1 size 3500 dif 4.40307e-12   // Dif is the difference with hand loops
[Blitz++] Precondition failure: Module 
/home/fperez/usr/local/blitz/include/blitz/array/eval.cc line 129
Assigned rank 1 expression to rank 2 array.
mtinner: /home/fperez/usr/local/blitz/include/blitz/array/eval.cc:129: 
blitz::Array<T, N>& blitz::Array<T, N>::evaluate(T_expr, T_update) [with 
T_expr = 
blitz::_bz_ArrayExpr<blitz::_bz_ArrayExprReduce<blitz::_bz_ArrayExpr<blitz::_bz_ArrayExprOp<blitz::_bz_ArrayExpr<blitz::ArrayIndexMapping<double, 
2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0> >, 
blitz::_bz_ArrayExpr<blitz::ArrayIndexMapping<double, 2, 1, 1, 0, 0, 0, 0, 0, 
0, 0, 0, 0> >, blitz::Multiply<double, double> > >, 1, 
blitz::ReduceSum<double, double> > >, T_update = blitz::_bz_update<double, 
double>, P_numtype = double, int N_rank = 2]: Assertion `0' failed.
Abort

Platform info: gcc 3.2.2, from a stock RedHat 9 install (x86).  Blitz from 
today's CVS.

I don't understand why it says that I'm assigning a rank 1 expression to a 
rank 2 array.  It seems pretty clear to me that in the 2d case, there's 2 free 
indices (i0,i1), so the resulting expression should be rank 2. No?

I'd very much appreciate any advice on this, as I'm rather new to C++ templates.

Regards,

Fernando.