As the attached example demonstrates, a chained reduction of an array expression can only be processed when the intermediate results are stored in an array. This is not feasible for large arrays.
#include <blitz/array.h>
using namespace blitz;
using namespace blitz::tensor;
#define EXPR a1(i, k) * a2(k, l) * a3(l, j)
int
main()
{
Array<float, 2> a1(2, 2), a2(2, 2), a3(2, 2);
Array<float, 4> a4(2, 2, 2, 2);
Array<float, 2> a5(2, 2), a6(2, 2);
a1 = 1, 0, 0, 1;
a2 = a1;
a3 = a1;
a4 = EXPR;
a5 = sum(sum(a4 , l), k); // this works
a6 = sum(sum(EXPR, l), k); // this fails, although it is semantically identical
return 0;
}
Further comments by Patrik Jonsson:
The problem seems to have something to do with (from the manual):
"Reductions have an important restriction: It is currently only possible to reduce over the last dimension of an array or array expression."
The reduction works if the expression is transposed so tensor::l is the last index of the last operand:
sum(a1(i, k) a2(l,k) a3(j, l), l)
However, transposing the operands shouldn't affect whether the expression can be reduced, so there's something strange going on here.
The problem is that compute_ordering() throws an assert if the ordering of the expression is indeterminate. In _bz_ArrayExprReduce::computeOrdering(), the ordering of the reduction is taken from the operands, except that missing values are just initialized in standard order, so just calling this for an expression with mixed order will throw the assert. It seems that mixed ordering of the operands here is not a fatal problem, since the expression will be evaluated using indices anyway. So perhaps if there's mixed order of the operands, we should just skip it and initialize that dim just as if it was missing?
Migrated from: https://sourceforge.net/p/blitz/bugs/17/
As reported by @mgrabner on 2008-08-18:
Further comments by Patrik Jonsson: