From: Stuart L. <sa...@il...> - 2014-08-20 17:11:30
|
Is there any concise Blitz way to do the equivalent of Array<float,3> A( Nz,Ny,Nx ); Array<float,1> yscale( Ny ); for i ... for j ... for k ... A(i,j,k) *= yscale(j) ? Do I just need to write it as an explicit loop over the 1-D axis? Or make a goofy (Nz,Ny,Nx) array with a hand-constructed stride vector that has the Z- and X-strides equal to zero? One reason I'd like this: trying to implement curl() over a grid with nonuniform spacing on each axis. If I can write this as an array expression, maybe it can fit in a stencil. Otherwise I'm making a mess of loops which looks a lot like the C style I was hoping to avoid... |
From: Patrik J. <co...@fa...> - 2014-08-20 17:21:40
|
Hi Stuart, I'm a bit rusty, but it seems like A = yscale(tensor::j); would work? Cheers, /Patrik On Wed, Aug 20, 2014 at 7:02 AM, Stuart Levy <sa...@il...> wrote: > Is there any concise Blitz way to do the equivalent of > > Array<float,3> A( Nz,Ny,Nx ); > Array<float,1> yscale( Ny ); > > for i ... > for j ... > for k ... > A(i,j,k) *= yscale(j) > > ? > > Do I just need to write it as an explicit loop over the 1-D axis? > > Or make a goofy (Nz,Ny,Nx) array with a hand-constructed stride vector > that has the Z- and X-strides equal to zero? > > One reason I'd like this: trying to implement curl() over a grid with > nonuniform spacing on each axis. > > If I can write this as an array expression, maybe it can fit in a > stencil. Otherwise I'm making a mess of loops which looks a lot like > the C style I was hoping to avoid... > > > ------------------------------------------------------------------------------ > Slashdot TV. > Video for Nerds. Stuff that matters. > http://tv.slashdot.org/ > _______________________________________________ > Blitz-support mailing list > Bli...@li... > https://lists.sourceforge.net/lists/listinfo/blitz-support |
From: Christopher S. <cs...@ma...> - 2014-08-20 18:23:17
|
On 08/20/14 13:21, Patrik Jonsson wrote: > Hi Stuart, > > I'm a bit rusty, but it seems like > > A = yscale(tensor::j); > > would work? It should mostly work, but with -DBZ_DEBUG the execution will fail with a rank error. The assertion checks that the rank of the left array precisely matches the rank of the right-hand array expression. This is overly restrictive, as the code actually works in practice without the debugging #define. Minimal example: $ cat bar.cpp #include <blitz/array.h> #include <iostream> int main() { blitz::Array<int,3> foo(2,2,2); blitz::Array<int,1> bar(2); bar = blitz::tensor::i; std::cout << "Running blitz++ version " << BZ_VERSION << std::endl; std::cout << "Rank 3 expression: " << std::endl; foo = bar(blitz::tensor::k); std::cout << foo << std::endl; std::cout << "Rank 2 expression: " << std::endl; foo = bar(blitz::tensor::j); std::cout << foo << std::endl; return 0; } In non-debug mode: $ g++ bar.cpp -lblitz; ./a.out Running blitz++ version 0.9 Rank 3 expression: (0,1) x (0,1) x (0,1) [ 0 1 0 1 0 1 0 1 ] Rank 2 expression: (0,1) x (0,1) x (0,1) [ 0 0 1 1 0 0 1 1 ] In debug mode: $ g++ bar.cpp -lblitz -DBZ_DEBUG; ./a.out Running blitz++ version 0.9 Rank 3 expression: (0,1) x (0,1) x (0,1) [ 0 1 0 1 0 1 0 1 ] Rank 2 expression: [Blitz++] Precondition failure: Module [...]/include/blitz/array/eval.cc line 129 Assigned rank 2 expression to rank 3 array. a.out: [...]/include/blitz/array/eval.cc:129: blitz::Array<P_numtype, N_rank>& blitz::Array<T, N>::evaluate(T_expr, T_update) [with T_expr = blitz::_bz_ArrayExpr<blitz::ArrayIndexMapping<blitz::FastArrayIterator<int, 1>, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0> >, T_update = blitz::_bz_update<int, int>, P_numtype = int, int N_rank = 3]: Assertion `0' failed. Aborted (core dumped) A debug-friendly workaround is to add +0*blitz::tensor::k (or the index corresponding to the highest dimension) in every tensor-using array expression. The 0*integer is a mathematical no-op, but it promotes the expression type to match the target's dimensionality. |
From: Stuart L. <sa...@il...> - 2014-08-20 18:14:28
|
Hmm, thanks, Patrik. That doesn't quite work - BZ_DEBUG reports that it's assigning a rank-2 expression to a rank-3 array. It does seem to work if I write all the tensor terms, like A *= ( 0*tensor::i + yscale( tensor::j ) + 0*tensor::k ); That does seem to create a 3-D expression with the required shape. But is it creating a whole 3D chunk of data to do it? I suspect it does and wish it wouldn't... It'd be very wasteful of cache, and I'm trying to do this efficiently on pretty large arrays. On 8/20/14 12:21 PM, Patrik Jonsson wrote: > Hi Stuart, > > I'm a bit rusty, but it seems like > > A = yscale(tensor::j); > > would work? > > Cheers, > > /Patrik > > > On Wed, Aug 20, 2014 at 7:02 AM, Stuart Levy <sa...@il...> wrote: >> Is there any concise Blitz way to do the equivalent of >> >> Array<float,3> A( Nz,Ny,Nx ); >> Array<float,1> yscale( Ny ); >> >> for i ... >> for j ... >> for k ... >> A(i,j,k) *= yscale(j) >> >> ? >> >> Do I just need to write it as an explicit loop over the 1-D axis? >> >> Or make a goofy (Nz,Ny,Nx) array with a hand-constructed stride vector >> that has the Z- and X-strides equal to zero? >> >> One reason I'd like this: trying to implement curl() over a grid with >> nonuniform spacing on each axis. >> >> If I can write this as an array expression, maybe it can fit in a >> stencil. Otherwise I'm making a mess of loops which looks a lot like >> the C style I was hoping to avoid... >> >> >> ------------------------------------------------------------------------------ >> Slashdot TV. >> Video for Nerds. Stuff that matters. >> http://tv.slashdot.org/ >> _______________________________________________ >> Blitz-support mailing list >> Bli...@li... >> https://lists.sourceforge.net/lists/listinfo/blitz-support > ------------------------------------------------------------------------------ > Slashdot TV. > Video for Nerds. Stuff that matters. > http://tv.slashdot.org/ > _______________________________________________ > Blitz-support mailing list > Bli...@li... > https://lists.sourceforge.net/lists/listinfo/blitz-support > |
From: Patrik J. <co...@fa...> - 2014-08-20 21:36:55
|
Yes, the tensor expressions must necessarily evaluate the rhs for every element of the lhs. I'm not aware of a way to take yscale and copy it into an arbitrary subset of A without a loop. I'm not even sure what the syntax for that would look like. On Wed, Aug 20, 2014 at 7:45 AM, Stuart Levy <sa...@il...> wrote: > Hmm, thanks, Patrik. > > That doesn't quite work - BZ_DEBUG reports that it's assigning a rank-2 > expression to a rank-3 array. > > It does seem to work if I write all the tensor terms, like > > A *= ( 0*tensor::i + yscale( tensor::j ) + 0*tensor::k ); > > That does seem to create a 3-D expression with the required shape. But > is it creating a whole 3D chunk of data to do it? I suspect it does > and wish it wouldn't... It'd be very wasteful of cache, and I'm trying > to do this efficiently on pretty large arrays. > > > On 8/20/14 12:21 PM, Patrik Jonsson wrote: >> Hi Stuart, >> >> I'm a bit rusty, but it seems like >> >> A = yscale(tensor::j); >> >> would work? >> >> Cheers, >> >> /Patrik >> >> >> On Wed, Aug 20, 2014 at 7:02 AM, Stuart Levy <sa...@il...> wrote: >>> Is there any concise Blitz way to do the equivalent of >>> >>> Array<float,3> A( Nz,Ny,Nx ); >>> Array<float,1> yscale( Ny ); >>> >>> for i ... >>> for j ... >>> for k ... >>> A(i,j,k) *= yscale(j) >>> >>> ? >>> >>> Do I just need to write it as an explicit loop over the 1-D axis? >>> >>> Or make a goofy (Nz,Ny,Nx) array with a hand-constructed stride vector >>> that has the Z- and X-strides equal to zero? >>> >>> One reason I'd like this: trying to implement curl() over a grid with >>> nonuniform spacing on each axis. >>> >>> If I can write this as an array expression, maybe it can fit in a >>> stencil. Otherwise I'm making a mess of loops which looks a lot like >>> the C style I was hoping to avoid... >>> >>> >>> ------------------------------------------------------------------------------ >>> Slashdot TV. >>> Video for Nerds. Stuff that matters. >>> http://tv.slashdot.org/ >>> _______________________________________________ >>> Blitz-support mailing list >>> Bli...@li... >>> https://lists.sourceforge.net/lists/listinfo/blitz-support >> ------------------------------------------------------------------------------ >> Slashdot TV. >> Video for Nerds. Stuff that matters. >> http://tv.slashdot.org/ >> _______________________________________________ >> Blitz-support mailing list >> Bli...@li... >> https://lists.sourceforge.net/lists/listinfo/blitz-support >> > > > ------------------------------------------------------------------------------ > Slashdot TV. > Video for Nerds. Stuff that matters. > http://tv.slashdot.org/ > _______________________________________________ > Blitz-support mailing list > Bli...@li... > https://lists.sourceforge.net/lists/listinfo/blitz-support |
From: Patrik J. <co...@fa...> - 2014-08-20 22:05:07
|
A clarification: blitz never creates temporary arrays (even if it would be advantageous to do so). But in this case it will *evaluate* the rhs for every value of i,j,k, even if not all of them appear in the expression. On Wed, Aug 20, 2014 at 11:36 AM, Patrik Jonsson <co...@fa...> wrote: > Yes, the tensor expressions must necessarily evaluate the rhs for > every element of the lhs. I'm not aware of a way to take yscale and > copy it into an arbitrary subset of A without a loop. I'm not even > sure what the syntax for that would look like. > > > On Wed, Aug 20, 2014 at 7:45 AM, Stuart Levy <sa...@il...> wrote: >> Hmm, thanks, Patrik. >> >> That doesn't quite work - BZ_DEBUG reports that it's assigning a rank-2 >> expression to a rank-3 array. >> >> It does seem to work if I write all the tensor terms, like >> >> A *= ( 0*tensor::i + yscale( tensor::j ) + 0*tensor::k ); >> >> That does seem to create a 3-D expression with the required shape. But >> is it creating a whole 3D chunk of data to do it? I suspect it does >> and wish it wouldn't... It'd be very wasteful of cache, and I'm trying >> to do this efficiently on pretty large arrays. >> >> >> On 8/20/14 12:21 PM, Patrik Jonsson wrote: >>> Hi Stuart, >>> >>> I'm a bit rusty, but it seems like >>> >>> A = yscale(tensor::j); >>> >>> would work? >>> >>> Cheers, >>> >>> /Patrik >>> >>> >>> On Wed, Aug 20, 2014 at 7:02 AM, Stuart Levy <sa...@il...> wrote: >>>> Is there any concise Blitz way to do the equivalent of >>>> >>>> Array<float,3> A( Nz,Ny,Nx ); >>>> Array<float,1> yscale( Ny ); >>>> >>>> for i ... >>>> for j ... >>>> for k ... >>>> A(i,j,k) *= yscale(j) >>>> >>>> ? >>>> >>>> Do I just need to write it as an explicit loop over the 1-D axis? >>>> >>>> Or make a goofy (Nz,Ny,Nx) array with a hand-constructed stride vector >>>> that has the Z- and X-strides equal to zero? >>>> >>>> One reason I'd like this: trying to implement curl() over a grid with >>>> nonuniform spacing on each axis. >>>> >>>> If I can write this as an array expression, maybe it can fit in a >>>> stencil. Otherwise I'm making a mess of loops which looks a lot like >>>> the C style I was hoping to avoid... >>>> >>>> >>>> ------------------------------------------------------------------------------ >>>> Slashdot TV. >>>> Video for Nerds. Stuff that matters. >>>> http://tv.slashdot.org/ >>>> _______________________________________________ >>>> Blitz-support mailing list >>>> Bli...@li... >>>> https://lists.sourceforge.net/lists/listinfo/blitz-support >>> ------------------------------------------------------------------------------ >>> Slashdot TV. >>> Video for Nerds. Stuff that matters. >>> http://tv.slashdot.org/ >>> _______________________________________________ >>> Blitz-support mailing list >>> Bli...@li... >>> https://lists.sourceforge.net/lists/listinfo/blitz-support >>> >> >> >> ------------------------------------------------------------------------------ >> Slashdot TV. >> Video for Nerds. Stuff that matters. >> http://tv.slashdot.org/ >> _______________________________________________ >> Blitz-support mailing list >> Bli...@li... >> https://lists.sourceforge.net/lists/listinfo/blitz-support |