[Perldl] threading example

[Chris Marshall]

One direction for improvement that I would like to see
for PDL is to simplify and clarify the use of PDL threading
to vectorize loops over multiple dimensions.  A recent
discussion on #pdl gave me an example to work from:

   http://irclog.perlgeek.de/pdl/2013-12-02

I've attached the input program and the script that
implements the CalcNormals computation with careful
annotations about the shape of the various arguments.
Getting the best layout of the dimension orders is key
to a compact implementation---however, due to the
density of pdl threading expressions, comprehension
can still be a problem.  In this case, the main result was
to convert this code:

   # Accumulate each triangle normal into each of the triangle vertices
   for ( $i = 0 ; $i < $IndexCount ; $i += 3 ) {

      $Index0 = $pIndices->at( $i );
      $Index1 = $pIndices->at( $i + 1 );
      $Index2 = $pIndices->at( $i + 2 );

      $v1     = $pVertices->slice( "0:2,($Index1)" ) -
$pVertices->slice( "0:2,($Index0)" );
      $v2     = $pVertices->slice( "0:2,($Index2)" ) -
$pVertices->slice( "0:2,($Index0)" );

      $Normal = $v1->crossp( $v2 )->norm;

      $pVertices->slice( "5:7,($Index0)" ) += $Normal;
      $pVertices->slice( "5:7,($Index1)" ) += $Normal;
      $pVertices->slice( "5:7,($Index2)" ) += $Normal;

   }

   # Normalize all the vertex normals
   for $i ( 0 .. $VertexCount - 1 ) {
      $pVertices->slice( "5:7,($i)" ) .= $pVertices->slice( "5:7,($i)" )->norm;
   }


into this:

   $tri = $Indices->splitdim(0,3)->copy;   # triangle/polys shape [
nC, nP ] = [ 3, 5 ]

   # calculate v10 and v20 for all triangles
   $v10 = $xyz->(:,$tri((1))) - $xyz->(:,$tri((0)));  #     shape [
nX, nP ] = [ 3, 5 ]
   $v20 = $xyz->(:,$tri((2))) - $xyz->(:,$tri((0)));  #     shape [
nX, nP ] = [ 3, 5 ]

   # calculate norms for all triangles,                     shape [
nX, nP ] = [ 5, 3 ]
   $N = $v10->crossp($v20)->norm;

   # use explicit loop to sum normals
   for my $i ( 0..2 ) { indadd $N(($i),*3)->flat, $tri->flat, $norm(($i)) }

   # renormalize to unit length
   $norm .= $norm->norm;


PDFs of the before and vectorized code are attached.  Comments and
suggestions for improved usability for threading with PDL are welcome.

Enjoy,
Chris

before.pdf
Description: Adobe PDF document

vectorized.pdf
Description: Adobe PDF document

_______________________________________________
Perldl mailing list
Perldl@jach.hawaii.edu
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl