develooper Front page | | Postings from September 2000

RFC 116 (v2) Efficient numerics with perl

Perl6 RFC Librarian
September 11, 2000 11:49
RFC 116 (v2) Efficient numerics with perl
Message ID:
This and other RFCs are available on the web at

=head1 TITLE

Efficient numerics with perl

=head1 VERSION

  Maintainer: pdl-porters team <>
  Date: 16 August 2000
  Last Modified: 11 September 2000
  Mailing List:
  Number: 116
  Version: 2
  Status: Developing


This RFC describes basic requirements of a numerical object/variable
type for perl. A reference implementation (that could almost certainly
do with some syntactical sweetening) is provided. It is currently
available as the PDL distribution. A few key concepts of the reference
implementation will be described. Any efforts to get numerical
support into the perl core should provide means to either easily
hook a numerical core into perl to achieve this or supply the
features at the core design level. Current implementation
aspects are sketched.

It is also the pdl-porters' contribution to I<overstuff everybodies
brains> [1].


PDL provides an object class to compactly store and speedily
manipulate the large N-dimensional data arrays which are the bread and
butter of scientific computing.

Apart from the object representation itself PDL implements
three key concepts which its users have found critical for
efficient and enjoyable numerical computing with perl:

=over 3

=item *

slices or I<smart references>

=item *

PDL threading -- automated optimized looping over elementary operations

=item *

Code generation and/or algorithmic abstraction


=head2 Compact typed multi-dimensional arrays

The basic data object of PDL is a zero- to multidimensional array 
of numbers, all of the same data type (for efficiency reasons, 
it is desirable to store the numbers in the native binary format and 
consecutively in memory instead of having a Perl SV for each one).
Perl6 could come up with a representation for such types as part
of its core (or be sufficiently flexible in its SV design).

It is also critical to be able to specify a range of types,
for example PDL offers support for single and double precision
(float/double) as well as a variety of integer types ('byte',
'ushort', 'short', 'long'). A future implementation should
support a complex type at the core level. In the current
implemenetation this is difficult due to parser limitations
of the code generator (C<PDL::PP>) and C's lack of an intrinsic
complex type.

=head2 Smart references and transformations: slicing and dicing

During data analysis it is often desirable to be able to refer to
parts of the data, for example I<every other number in this vector>
without allocating new storage: if the original vector is 100Mb long,
we don't want to allocate and copy the data for another 50Mb
vector. For this purpose PDL supports "virtual piddles", which are
pointers to parts of other piddles (the implementation is transparent,
virtual piddles appearing just as any other piddle object to the

  $a = zeroes(1000,1000,100); # more than 100MB worth of data
  $b = $a->slice('0:-1:2');   # virtual piddle of every 2nd element
  $b *= 2;                    # changes the elements in $a

C<slice> is just one example of this type of operation. C<$b> contains
only information how to access C<$a>'s data rather than a copy. Similarly,
a transposed object with the first two dimensions exchanged is easily

  $c = $a->xchg(0,1); # exchange dimension 0 and 1

Again, C<$c> only contains information how to access C<$a>'s data.

Similarly we can create an index:

  $i = pdl(1,9,27); # a number of indices
  $c = $a->index($i);

Changing C<$c> changes C<$a> (and vice versa). However in this
case PDL makes no nice optimisation and manipulating C<$c> uses memory
for temps even though the changes are propagated back. The crucial point
is that C<$c> and C<$a> are tightly linked objects that propagate
changes in either to the other (in a way similar to dataflow). Find
out more about this in the implementation section.
An optimisation is in principle possible, and would be useful
for all array types not just PDL ones.

Similar optimizations for sparse and non-rectangular arrays would
be desirable, this would be especially useful in a perl for numerics.

=head2 Signatures: threading over elementary operations

In a large number of practical situations encountered in scientific
numerical programming fundamental (simple) operations are
iterated over many (small) chunks of (multi-dimensional) data.

A very simple examples is addition

  $c = $a + $b

of two multi-dimensional arrays (PDL objects). The fundamental
operation is the addition of two scalar elements. This fundamental
operation is just iterated over all pairs of elements of C<$a> and
C<$b> and results placed element by element into C<$c>.

Based on this observation PDL introduced the notion of I<PDL threading>
to implement such functionality in a flexible and efficient way:
it provides tools (C<PDL::PP>) to define those optimized innermost
loops (the fundamental operations) and establishes a framework to
efficiently execute loops over these fundamental operations
automatically as prescribed by the dimensionality of the input data. Most
importantly, those loops over the innermost fundamental operations
can generally be run in parallel (i.e. run in separate I<threads>).
PDL currently has experimental support to exploit this observation.

Two examples using the inner product as the fundamental
operation may serve to illustrate the idea. To convert an RGB image into a
greyscale image an inner product with a triple of weights with
each RGB triple is performed:

  $grey = inner $rgb, pdl(0.301,0.586,0.113);

The key is here that regardless of the number of trailing dimensions
of the RGB object C<$rgb> (which has a shape [3,.....]) the
fundamental operation (the inner product) is speedily iterated
over all triples of the data. So C<$rgb> could be a single RGB
pixel (shape [3]), a line of RGB data (shape [3,x]), an RGB image
(shape [3,x,y]), a stack of RGB images ([3,x,y,z]) and so on.

Each function that can be used in such a way in PDL has a I<signature>
that symbolically states dimensions consumed and generated in a fundamental
operation. Illustrating this by the two functions we have used
above the signatures of addition and inner product are

  a(); b(); [o] c()    # addition, we use ';'s as separators
  a(n); b(n); [o] c()  # inner product  

In this symbolic notation C<[o]> flags arguments which
are created (output arguments). The symbols in parentheses name
any dimensions of the variables involved. Empty parentheses refer to
scalars (0-dim). So addition takes two scalar input arguments
and generates one output scalar. The inner product consumes, as
expected, two one-dimensional input args (symbolized by C<a(n)>
and C<b(n)>) to produce a scalar output (C<[o] c()>).

To complete the picture PDL has rules how
these fundamental operations consume chunks of data of the
input args to produce a muti-dimensional result of chunks
of output data (the possibly multi-dimensional result piddle).
The user decides exactly how those operations are
iterated over his data by manipulating the dimensions
of the input data appropriately.
That is illustrated in the promised second example.

Those familar with linear algebra will know that
a matrix product consists of the result of inner products
between all pairs of row and column vectors of two
input matrices. Not surprisingly we can therefore use the
inner product to compute the matrix product of two PDL
objects C<$a> and C<$b>:

 $c =  inner $a->dummy(1), $b->xchg(0,1)->dummy(2);

C<dummy> and C<xchg> are dimension manipulation functions
(of the smart reference type) that are used to change the
dimensionality of the input piddles so that the inner
product is iterated over the desired pairs of sub-vectors
and placed into the resulting piddle C<$c>. We do not expect
the casual reader to follow the details of this example
but to appreciate the underlying idea.

The essence of I<PDL threading> is here that (1) it is a useful
feature in I<many> practical situations (2) can be very efficiently
implemented and (3) exposes parallelism in the underlying algorithm.

One could go as far as saying that the ideas underlying
I<threading> are PDL's main programming paradigm and
the user is expected to rethink his approach accordingly for
optimal performance.

=head2 Defining new PDL functions -- Glue code generation

One of the foremost design criteria was to implement a small set of
core structures, routines and code generators that would depend
heavily on each other. The rest of the package (implementing the real
numerical functionality) should only depend on the external interface
of these core functions.

Encapsulation can be achieved by consciously looking for
code that would be repeated in subroutines and either making a new
subroutine or a code generator do the job for the programmer. This 
way, when the repeated code suddenly needs to be different, there
is no painstaking need of changing all instances of it - simply
changing the subroutine or code generator will do it. We emphasize this 
because this kind of thinking seems to be forgotten all too often.

PDL introduced the code generator PDL::PP to create
subroutines from concise descriptions to allow the programmer
to focus on the algorithmic side of things. In many ways PDL::PP
can be thought of as the XS of PDL: it is complex to learn, powerful,
allows easy interfacing at the C-level and is, once mastered, straightforward
to use. Beyond that using PDL::PP ensures that new routines follow uniform
calling conventions and support PDL threading, slicing, etc.

In the current implementation it suffers from similar problems as XS
in perl5: the syntax is strange and hard to learn, it is also feature
incomplete (i.e. certain things can not be done right now). However,
it implements a number of ideas that we find are absolutely crucial
to pertain and develop further in a numerically friendly perl:
algorithmic abstraction, encapsulation and provision of a tool
for easy interfacing to the hundreds of existing numerical, graphics
and IO libraries that are to be tapped into.

As an illustration for the abstraction that is achieved, here is the
definition for an inner product function. This definition tells the
code generator (C<PDL::PP>) to generate the appropriate XS code that
will then be compiled as part of a new loadable PDL module:

        Pars => 'a(n); b(n); [o]c(); ', # the signature, see above
        Code => 'double tmp = 0;
                 loop(n) %{ tmp += $a() * $b(); %}
                 $c() = tmp;' );

The above code specifies firstly that the two first vectors must have
the same first dimension and the third vector will not have that
dimension. Inside the loop it is possible to use C<$a()> without
specifying the index because we have specified that we are looping
over the dimension called C<n>. Notice that we do not do any
checking that C<$a> and C<$b> have the same sized first
dimension - we just tell PDL::PP that we expect this to be so and it
checks whether any data passed to the C<inner> routine obeys this
assumption and C<croak>s if it does not.

After installation the user can then use

	$c = inner $a,$b;

to compute the inner product of two piddles. 

The code clearly illustrates the abstraction achieved, making no
reference to an actual data type or how the vectors are stored. This
kind of abstraction is similar to the ideas used in the C++ Standard
Template Library.  As a bonus, the use of PP means a
lot of the tedious housekeeping of programming can be automated.
Almost all index troubles vanish when writing and using these
subroutines: all the necessary index checks are done automatically by
the code generator.

The code generator also ensures that calling the function with
arguments that have more dimensions than the function expects works as
well - the function is simply repeatedly called in a loop over all
these dimensions (I<threaded> over the dimensions in PDL parlance),
providing a fast way of performing many similar operations in a row
without going back to the Perl level (with the latest versions of PDL,
it is possible to use multiple OS-level threads for these loops,
enabling an easy way of gaining from a multiprocessor computer). This
I<threading> allows many algorithms to be expressed concisely and PDL
scripts rarely have cause for loops over elements on the Perl level.

In a future perl this functionality might be integrated with a perl (JIT?)
compiler. Code could be specified in perl(-like?) syntax and compiled into an
optimized C level transformation as required.

If this sounds rather vague it is not by accident -- we still lack a clear
idea how to optimally achieve this goal. A nice perlish syntax to
replace the PP gibberish should certainly be the starting point.


The PDL core is currently suffering from similar problems as
perl 5: it works but is poorly documented and many hacks make it
increasingly difficult to modify and maintain. But it works C<;)>.

In any case, some of the perl6 guys have asked for a rundown of
the principal guts of the implementation. Here it goes...

=head2 Compact, typed multi-dimensional arrays

On the Perl side, an opaque scalar reference turned out to be the most
practical representation: Perl5 allows overloading of operators so
C<$c=$a+$b> can mean "add the objects (here representing matrices)
C<$a> and C<$b> and store the result in the scalar C<$c> with whatever
type that the result has".  As with any perl object, new objects can
be easily derived from piddles for specialised purposes.

In Perl6 it may be desirable to use C<@syntax> for all arrays - whether
compact, implementing smart references, or not. The use of C<$a> for
PDL arrays and C<@a> for normal perl arrays B<is> confusing, we would
like to see some unification. Ideally perl arrays would support some,
or all, of the PDL features whether they are multi-dimensional
or not.

=over 5

=item What is in a Piddle.

The PDL structure is defined in C<pdl.h>

	struct pdl {
	   unsigned long magicno; /* Always stores PDL_MAGICNO as a sanity check */
	     /* This is first so most pointer accesses to wrong type are caught */
	   int state;        /* What's in this pdl */

	   pdl_trans *trans; /* Opaque pointer to internals of transformation from
				parent */

	   pdl_vaffine *vafftrans;

	   void*    sv;      /* (optional) pointer back to original sv.
				  ALWAYS check for non-null before use.
				  We cannot inc refcnt on this one or we'd
				  never get destroyed */

	   void *datasv;        /* Pointer to SV containing data. Refcnt inced */
	   void *data;            /* Null: no data alloced for this one */
	   int nvals;           /* How many values allocated */
	   int datatype;
	   PDL_Long   *dims;      /* Array of data dimensions */
	   PDL_Long   *dimincs;   /* Array of data default increments */
	   short    ndims;     /* Number of data dimensions */

	   unsigned char *threadids;  /* Starting index of the thread index set n */
	   unsigned char nthreadids;

	   pdl *progenitor; /* I'm in a mutated family. make_physical_now must
			       copy me to the new generation. */
	   pdl *future_me;  /* I'm the "then" pdl and this is my "now" (or more modern
			       version, anyway */

	   pdl_children children;

	   short living_for; /* perl side not referenced; delete me when */

	   PDL_Long   def_dims[PDL_NDIMS];   /* Preallocated space for efficiency */
	   PDL_Long   def_dimincs[PDL_NDIMS];   /* Preallocated space for efficiency */
	   unsigned char def_threadids[PDL_NTHREADIDS];

	   struct pdl_magic *magic;

	   void *hdrsv; /* "header", settable from outside */

This is quite a structure for just storing some data in - what is going on?

=item Data storage

We are going to start with some of the simpler members: first of all,
there is the member

	void *datasv;

which is really a pointer to a Perl SV structure (C<SV *>). The SV is
expected to be representing a string, in which the data of the piddle
is stored in a tightly packed form. This pointer counts as a reference
to the SV so the reference count has been incremented when the C<SV *>
was placed here. This pointer is allowed to have the value C<NULL> which 
means that there is no Perl SV for this data - for instance, the data
might be allocated by a C<mmap> operation. Note the use of an SV*
was purely for convenience, it allows easy transformation of
packed data from files into piddles. Other implementations are not

The actual pointer to data is stored in the member

	void *data;

which contains a pointer to a memory area with space for

	int nvals;

data items of the data type of this piddle.

The data type of the data is stored in the variable

	int datatype;

the values for this member are given in the enum C<pdl_datatypes>, discussed

=item Dimensions

The number of dimensions in the piddle is given by the member

	int ndims;

which shows how many entries there are in the arrays

	PDL_Long   *dims;      
	PDL_Long   *dimincs;

These arrays are intimately related: C<dims> gives the sizes of the dimensions
and C<dimincs> is always calculated by the code

	int inc = 1;
        for(i=0; i<it->ndims; i++) {
		it->dimincs[i] = inc; inc *= it->dims[i];

in the routine C<pdl_resize_defaultincs> in C<pdlapi.c>.
What this means is that the dimincs can be used to calculate the offset
by code like

	int offs = 0;
	for(i=0; i<it->ndims; i++) {
		offs += it->dimincs[i] * index[i];

but this is not always the right thing to do,
at least without checking for certain things first.

=item Default storage

Since the vast majority of piddles don't have more than 6 dimensions,
it is more efficient to have default storage for the dimensions and dimincs
inside the PDL struct.

   	PDL_Long   def_dims[PDL_NDIMS];   
   	PDL_Long   def_dimincs[PDL_NDIMS]; 

The C<dims> and C<dimincs> may be set to point to the beginning of these
arrays if C<ndims> is smaller than or equal to the compile-time constant
C<PDL_NDIMS>. This is important to note when freeing a piddle struct.
The same applies for the threadids:

   	unsigned char def_threadids[PDL_NTHREADIDS];

=item Magic

It is possible to attach magic to piddles, much like Perl's own magic
mechanism. If the member pointer

	   struct pdl_magic *magic;

is nonzero, the PDL has some magic attached to it. The implementation
of magic is discussed below.

=item State

One of the first members of the structure is 

	int state;

The possible flags and their meanings are given in C<pdl.h>.

=item Transformations and virtual affine transformations

As you should already know, piddles often carry information about
where they come from. For example, the code

	$b = $a->slice("2:5");
	$b .= 1;

will alter $a. This information is stored in the members

   	pdl_trans *trans; 
   	pdl_vaffine *vafftrans;

C<pdl_trans> and C<pdl_vaffine> are structures that we will look at in
more detail below.

=item The Perl SVs

When piddles are referred to through Perl SVs, we store an additional
reference to it in the member

	void*    sv;

in order to be able to return a reference to the user when he wants to 
inspect the transformation structure on the Perl side.

Also, we store an opaque

	void *hdrsv; 

which is just for use by the user to hook up arbitrary data with this sv.


=head2 Smart references and transformations: slicing and dicing

Smart references and most other fundamental functions
operating on piddles are implemented via I<transformations>
which are represented by the type C<pdl_trans> in PDL.

A transformation links input and output piddles and contains
all the infrastructure that defines how

=over 4

=item *

output piddles are obtained from input piddles

=item *

changes in smartly linked output piddles (e.g. the I<child>
of a sliced I<parent> piddle) are flown back to the input
piddle in transformations where this is supported (the most
often used example being C<slice> here).

=item *

datatype and size of output piddles that need to be created
are obtained


In general, executing a PDL function on a group of piddles
results in creation of a transformation of the requested
type that links all input and output arguments (at least
those that are piddles). In PDL functions that support
data flow between input and output args (e.g. C<slice>,
C<index>) this transformation links I<parent> (input) and
I<child> (output) piddles permanently until either the link is
explicitly broken by user request (C<sever> at the perl level)
or all parents and childen have been destroyed. In those
cases the transformation is lazy-evaluated, e.g. only executed
when piddle values are actually accessed.

In I<non-flowing> functions, for example addition (C<+>) and inner
products (C<inner>), the transformation is installed just as
in flowing functions but then the transformation is immediately
executed and destroyed (breaking the link between input and output args)
before the function returns.

It should be noted that the close link between input and output args
of a flowing function requires that piddle objects that are linked in
such a way be kept alive beyond the point where they have gone
out of scope from the point of view of perl:

  $a = zeroes(20);
  $b = $a->slice('2:4');
  undef $a;    # last reference to $a is now destroyed

Although $a should now be destroyed according to perl's rules
the underlying C<*pdl> must actually only be freed when C<$b>
also goes out of scope (since it still references
internally some of C<$a>'s data). This example demonstrates that such
a dataflow paradigm between PDL objects necessitates a special
destruction algorithm that takes the links between piddles
into account and couples the lifespan of those objects. The
non-trivial algorithm is implemented in the function
C<pdl_destroy> in F<pdlapi.c>.

=head2 What's in a C<pdl_trans>

All transformations are implemented as structures

  struct XXX_trans {
	int magicno; /* to detect memory overwrites */
	short flags; /* state of the trans */
	pdl_transvtable *vtable;   /* the all important vtable */
	void (*freeproc)(struct pdl_trans *);  /* Call to free this trans
		(in case we had to malloc some stuff dor this trans) */
        pdl *pdls[NP]; /* The pdls involved in the transformation */
	int __datatype; /* the type of the transformation */
        /* in general more members
        /* depending on the actual transformation (slice, add, etc)

The transformation identifies all C<pdl>s involved in the trans

  pdl *pdls[NP];

with C<NP> depending on the number of piddle args of the particular
trans. It records a state

  short flags;

and the datatype

  int __datatype;

of the trans (to which all piddles must be converted unless
they are explicitly typed). Most important is the pointer to
the vtable that contains the actual functionality

 pdl_transvtable *vtable;

The vtable structure in turn looks something like (slightly
simplified from F<pdl.h> for clarity)

  typedef struct pdl_transvtable {
	pdl_transtype transtype;
	int flags;
	int nparents;   /* number of parent pdls (input) */
	int npdls;      /* number of child pdls (output) */
	char *per_pdl_flags;  /* optimization flags */
	void (*redodims)(pdl_trans *tr);  /* figure out dims of children */
	void (*readdata)(pdl_trans *tr);  /* flow parents to children  */
	void (*writebackdata)(pdl_trans *tr); /* flow backwards */
	void (*freetrans)(pdl_trans *tr); /* Free both the contents and it of
					the trans member */
	pdl_trans *(*copy)(pdl_trans *tr); /* Full copy */
  	int structsize;
	char *name; /* For debuggers, mostly */
  } pdl_transvtable;

We focus on the callback functions:

  	void (*redodims)(pdl_trans *tr);

C<redodims> will work out the dimensions of piddles that need
to be created and is called from within the API function that
should be called to ensure that the dimensions of a piddle are
accessible (F<pdlapi.c>):

   void pdl_make_physdims(pdl *it)

C<readdata> and C<writebackdata> are responsible for the actual
computations of the child data from the parents or parent data
from those of the children, respectively (the dataflow aspect).
The PDL core makes sure that these are called as needed when
piddle data is accessed (lazy-evaluation). The general API
function to ensure that a piddle is up-to-date is

  void pdl_make_physvaffine(pdl *it)

which should be called before accessing piddle data from
XS/C (see F<Core.xs> for some examples).

C<freetrans> frees dynamically allocated memory associated
with the trans as needed and C<copy> can copy the transformation.

The transformation and vtable code is hardly ever written by
hand but rather generated by C<PDL::PP> from concise descriptions
(see the example of the definition of C<inner> above).

Certain types of transformations can be optimized very
efficiently obviating the need for explicit C<readdata>
and C<writebackdata> methods. Those transformations are
called I<pdl_vaffine>. Most dimension manipulating
functions (e.g., C<slice>, C<xchg>) belong to this class.

The basic trick is that parent and child of such a transformation work
on the same (shared) block of data which they just choose
to interpret differently (by dusing different C<dims>, C<dimincs> and
C<offs> on the same data, compare the C<pdl> structure above).
Each operation on a piddle sharing
data with another one in this way is therefore automatically flown
from child to parent and back -- they are reading and writing
the same block of memory. This is currently not perl thread safe --
no big loss since the whole PDL core is not reentrant
(perl threading C<!=> PDL threading).

A more detailed description of I<vaffines> will be added
in a later version of this document.

=head2 Signatures: threading over elementary operations

Most of that functionality is implemented in the file F<pdlthread.c>.
The C<PDL::PP> generated functions (in particular the
C<readdata> and C<writebackdata> callbacks) use this infrastructure to 
make sure that the fundamental operation implemented by the
trans is performed in agreement with PDL's threading semantics.

=head2 Defining new PDL functions -- Glue code generation

Please, see L<PDL::PP> and examples in the PDL distribution. Implementation
and syntax are currently far from perfect.

As the above example on the definition of C<inner> suggests,
the principal functionality
is implemented through the function C<pp_def> which takes a hash
of directives as a concise specification of the desired functionality.
The C<Code> and C<RedoDimsCode> fields contain the algorithmic specification
in a home grown pseudo code syntax that has some faint similarity to
a mixture of perl, XS and C. Here we could do with some good ideas
for improvement. From this description C<PDL::PP> produces C and XS code
(which will go away in a future perl?) by processing the hash through a
translation table.

The ability to write the fundamental algorithm in a perlish way
(or even in perl) that is optimized by (ideally JIT) compilation
is what we are really after.

=head1 SUMMARY

We have explained the key ideas which make PDL work, and work well. It remains
to be determined which of these might belong in the core of perl6, in order
to make PDL work more seamlessly, and which belong in the module.

Our take on this: perl6 should support a compact, numerical array type which is
overloadable and can be manipulated with nice syntax. A key point is to abolish the
visual distinction between perl arrays (C<@x>) and PDL arrays (C<$x>) which is confusing
but forced upon us by the current object paradigm. There is also considerable
scope for improving syntax to make numerical Perl more user-friendly - this is explained
in more detail by the RFCs on Ranges and Default Methods.

Compact arrays should support virtual slices when being manipulated.

It seems to us that the more detailed/generalised PDL threading paradigm is
beyond the core and belongs in a module. Some syntactical hooks would be
appreciated though so one could at least say code like this at the Perl level:

        sub mysub : PDL {
                my($a, $b) = @_;
                signature('a(n), b(n)');
                loop(n) {   # n a bareword, loop has prototype
                # Threaded code

=head1 SEE ALSO

L<PDL> (

perl6 RFCs on Ranges and Default Methods.


IDL (the one from rsinc)




Hords of other array-oriented languages, try

[1] LW:

    Let me just add that I don't mind the brainstorming at all.  To be a
    good language designer, you have to stuff your brain with what you
    *could* do before you can reasonably choose what you *will* do.  At the
    moment, I'm not only trying to follow along here; I'm also reading all
    the books on computer languaes I can get my hands on--not just to look
    for ideas to steal, but also to remind myself of the mindset Perl was
    designed to escape.

    In fact, I'd go as far as to say that it's imperative that you
    overstuff your brain to the point where you no longer feel tempted to
    overstuff the language.  Hypotheticality is your friend. Perl Programming lists via nntp and http.
Comments to Ask Bjørn Hansen at | Group listing | About