|
|
Blitz Support : |
From: Fernando Perez (fperez_at_[hidden])
Date: 2003-08-29 01:49:58
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.