=head1 NAME
PDL::Internals - description of some aspects of the current internals
=head1 DESCRIPTION
=head2 Intro
This document explains various aspects of the current implementation
of PDL. If you just want to
use
PDL
for
something, you definitely
do
not need to
read
this. Even
if
you want to interface your C routines
to PDL or create new L<PDL::PP> functions, you
do
not need to
read
this
man page (though it may be informative). This document is primarily
intended
for
people interested in debugging or changing the internals
of PDL. To
read
this, a good understanding of the C language
and programming and data structures in general is required, as well
as some Perl understanding. If you
read
through this document
and understand all of it and are able to point what any part of
this document refers to in the PDL core sources and additionally
struggle to understand L<PDL::PP>, you will be awarded the title
"PDL Guru"
(of course, the current version of this document
is so incomplete that this is
next
to impossible from just these notes).
B<Warning:> If it seems that this document
has
gotten out of date,
please inform the PDL porters email list (pdl-devel
@lists
.sourceforge.net).
This may well happen.
=head2 ndarrays
The pdl data object is generally an opaque
scalar
reference into a
pdl structure in memory. Alternatively, it may be a hash reference
with
the C<PDL> field containing the
scalar
reference (this makes overloading
ndarrays easy, see L<PDL::Objects>). You can easily find out
at the Perl level which type of ndarray you are dealing
with
. The example
code below demonstrates how to
do
it:
die
"not an ndarray"
unless
UNIVERSAL::isa(
$pdl
,
'PDL'
);
if
(UNIVERSAL::isa(
$pdl
,
"HASH"
)) {
die
"not a valid PDL"
unless
exists
$pdl
->{PDL} &&
UNIVERSAL::isa(
$pdl
->{PDL},
'PDL'
);
print
"This is a hash reference,"
,
" the PDL field contains the scalar ref\n"
;
}
else
{
print
"This is a scalar ref that points to address $$pdl in memory\n"
;
}
The
scalar
reference points to the numeric address of a C structure of
type C<pdl> which is
defined
in F<pdl.h>. The mapping between the
object at the Perl level and the C structure containing the actual
data and structural that makes up an ndarray is done by the PDL typemap.
The functions used in the PDL typemap are
defined
pretty much at the top
of the file F<pdlcore.h>. So what does the structure look like:
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 */
PDL_Indx nvals; /* How many
values
allocated */
int
datatype;
PDL_Indx
*dims
; /* Array of data dimensions */
PDL_Indx
*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_children children;
PDL_Indx def_dims[PDL_NDIMS]; /* Preallocated space
for
efficiency */
PDL_Indx 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?
=over 5
=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 ndarray
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 reference count business
has
to
do
with
Perl's
garbage collection mechanism -- don
't worry if this doesn'
t mean much
to you). This pointer is allowed to have the value C<NULL> which
means that there is
no
actual 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 ndarrays. Other implementations are not
excluded.
The actual pointer to data is stored in the member
void
*data
;
which contains a pointer to a memory area
with
space
for
PDL_Indx nvals;
data items of the data type of this ndarray. PDL_Indx is either
'long'
or
'long long'
depending on whether your perl is 64bit or not.
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> (see
F<pdl.h>). Currently we have byte, short, unsigned short, long,
index
(either long or long long), long long, float and double types, see
also L<PDL::Types>.
=item Dimensions
The number of dimensions in the ndarray is
given
by the member
int
ndims;
which shows how many entries there are in the arrays
PDL_Indx
*dims
;
PDL_Indx
*dimincs
;
These arrays are intimately related: C<dims> gives the sizes of the dimensions
and C<dimincs> is always calculated by the code
PDL_Indx 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
PDL_Indx 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 ndarrays 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_Indx def_dims[PDL_NDIMS];
PDL_Indx 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 an ndarray struct.
The same applies
for
the threadids:
unsigned char def_threadids[PDL_NTHREADIDS];
=item Magic
It is possible to attach magic to ndarrays, 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 can be gleaned from the file F<pdlmagic.c> in the distribution.
=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>.
These are mainly used to implement the lazy evaluation mechanism
and keep track of ndarrays in these operations.
=item Transformations and virtual affine transformations
As you should already know, ndarrays often carry information about
where they come from. For example, the code
$y
=
$x
->slice(
"2:5"
);
$y
.= 1;
will alter
$x
. So C<
$y
> and C<
$x
> I<know> that they are connected
via a C<slice>-transformation. This information is stored in the members
pdl_trans
*trans
;
pdl_vaffine
*vafftrans
;
Both C<
$x
> (the I<parent>) and C<
$y
> (the child) store this information
about the transformation in appropriate slots of the C<pdl> structure.
C<pdl_trans> and C<pdl_vaffine> are structures that we will look at in
more detail below.
=item The Perl SVs
When ndarrays 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.
This one is generally manipulated through L<sethdr|PDL::Core/sethdr> and
L<gethdr|PDL::Core/gethdr> calls.
=back
=head2 Smart references and transformations: slicing and dicing
Smart references and most other fundamental functions
operating on ndarrays are implemented via I<transformations>
(as mentioned above) which are represented by the type C<pdl_trans> in PDL.
A transformation links input and output ndarrays and contains
all the infrastructure that defines how:
=over 4
=item *
output ndarrays are obtained from input ndarrays;
=item *
changes in smartly linked output ndarrays (e.g. the I<child>
of a sliced I<parent> ndarray) are flown back to the input
ndarray in transformations where this is supported (the most
often used example being C<slice> here);
=item *
datatype and size of output ndarrays that need to be created
are obtained.
=back
In general, executing a PDL function on a group of ndarrays
results in creation of a transformation of the requested
type that links all input and output arguments (at least
those that are ndarrays). 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) ndarrays permanently
until
either the
link
is
explicitly broken by user request (C<sever> at the Perl level)
or all parents and children have been destroyed. In those
cases the transformation is lazy-evaluated, e.g. only executed
when
ndarray
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 (like L<slice|PDL::Slices/slice>) requires
that ndarray 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:
$x
= zeroes(20);
$y
=
$x
->slice(
'2:4'
);
undef
$x
;
Although
$x
should now be destroyed according to Perl's rules
the underlying C<pdl> structure must actually only be freed
when
C<
$y
>
also goes out of scope (since it still references
internally some of C<
$x
>'s data). This example demonstrates that such
a dataflow paradigm between PDL objects necessitates a special
destruction algorithm that takes the links between ndarrays
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>. In fact, most of the code
in F<pdlapi.c> and F<pdlfamily.c> is concerned
with
making sure that ndarrays (C<pdl *>s) are created, updated
and freed at the right
times
depending on interactions
with
other ndarrays via PDL transformations (remember, C<pdl_trans>).
=head2 Accessing children and parents of an ndarray
When ndarrays are dynamically linked via transformations as
suggested above input and output ndarrays are referred to as parents
and children, respectively.
An example of processing the children of an ndarray is provided
by the C<baddata> method of L<PDL::Bad>.
Consider the following situation:
pdl>
$x
= rvals(7,7,{
Centre
=>[3,4]});
pdl>
$y
=
$x
->slice(
'2:4,3:5'
);
pdl> ? vars
PDL variables in
package
main::
Name Type Dimension Flow State Mem
----------------------------------------------------------------
$x
Double D [7,7] P 0.38Kb
$y
Double D [3,3] -C 0.00Kb
Now,
if
I suddenly decide that C<
$x
> should be flagged as possibly
containing bad
values
, using
pdl>
$x
->badflag(1)
then I want the state of C<
$y
> - it's I<child> - to be changed as
well (since it will either share or inherit some of C<
$x
>'s data and
so be also I<bad>), so that I get a
'B'
in the I<State> field:
pdl> ? vars
PDL variables in
package
main::
Name Type Dimension Flow State Mem
----------------------------------------------------------------
$x
Double D [7,7] PB 0.38Kb
$y
Double D [3,3] -CB 0.00Kb
This bit of magic is performed by the C<propagate_badflag> function,
which is listed below:
/* newval = 1 means set flag, 0 means clear it */
/* thanks to Christian Soeller
for
this */
void propagate_badflag( pdl
*it
,
int
newval ) {
PDL_DECL_CHILDLOOP(it)
PDL_START_CHILDLOOP(it)
{
pdl_trans
*trans
= PDL_CHILDLOOP_THISCHILD(it);
int
i;
for
( i = trans->vtable->nparents;
i < trans->vtable->npdls;
i++ ) {
pdl
*child
= trans->pdls[i];
if
( newval ) child->state |= PDL_BADVAL;
else
child->state &= ~PDL_BADVAL;
/* make sure we propagate to grandchildren, etc */
propagate_badflag( child, newval );
} /*
for
: i */
}
PDL_END_CHILDLOOP(it)
} /* propagate_badflag */
Given an ndarray (C<pdl
*it
>), the routine loops through
each
C<pdl_trans> structure, where access to this structure is provided by the
C<PDL_CHILDLOOP_THISCHILD> macro.
The I<children> of the ndarray are stored in the C<pdls> array,
after
the
I<parents>, hence the loop from C<i = ...nparents> to
C<i = ...npdls - 1>.
Once we have the pointer to the child ndarray, we can
do
what we want to
it; here we change the value of the C<state> variable, but the details
are unimportant).
What B<is> important is that we call C<propagate_badflag> on this
ndarray, to ensure we loop through its children. This recursion
ensures we get to all the I<offspring> of a particular ndarray.
Access to I<parents> is similar,
with
the C<
for
> loop replaced by:
for
( i = 0;
i < trans->vtable->nparents;
i++ ) {
/*
do
stuff
with
parent
}
=head2 What's in a transformation (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
for
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 ndarray args of the particular
trans. It records a state
short flags;
and the datatype
int
__datatype;
of the trans (to which all ndarrays must be converted
unless
they are explicitly typed, PDL functions created
with
L<PDL::PP>
make sure that these conversions are done as necessary). Most important is
the pointer to the vtable (virtual table) 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 ndarrays that need
to be created and is called from within the API function that
should be called to ensure that the dimensions of an ndarray 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
ndarray data is accessed (lazy-evaluation). The general API
function to ensure that an ndarray is up-to-date is
void pdl_make_physvaffine(pdl
*it
)
which should be called
before
accessing ndarray 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.
Again, functions built
with
L<PDL::PP> make sure that copying
and freeing via these callbacks happens at the right
times
. (If they
fail to
do
that we have got a memory leak -- this
has
happened in
the past ;).
The transformation and vtable code is hardly ever written by
hand but rather generated by L<PDL::PP> from concise descriptions.
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 using different C<dims>, C<dimincs> and
C<offs> on the same data, compare the C<pdl> structure above).
Each operation on an ndarray sharing
data
with
another one in this way is therefore automatically flown
from child to parent and back --
after
all 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!).
=head2 Signatures: threading over elementary operations
Most of that functionality of PDL threading (automatic iteration
of elementary operations over multi-dim ndarrays) is implemented in the
file F<pdlthread.c>.
The L<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 but it does a good job!
=head2 The Core struct
As discussed in L<PDL::API>, PDL uses a pointer to a structure
to allow PDL modules access to its core routines. The definition of this
structure (the C<Core> struct) is in F<pdlcore.h> (created by
F<pdlcore.h.PL> in F<Basic/Core>) and looks something like
/* Structure to hold pointers core PDL routines so as to be used by
* many modules
*/
struct Core {
I32 Version;
pdl* (
*SvPDLV
) ( SV* );
void (
*SetSV_PDL
) ( SV
*sv
, pdl
*it
);
pdl* (
*pdlnew
) ( );
pdl* (
*tmp
) ( );
pdl* (
*create
) (
int
type);
void (
*destroy
) (pdl
*it
);
...
}
typedef struct Core Core;
The first field of the structure (C<Version>) is used to ensure
consistency between modules at run
time
; the following code
is placed in the BOOT section of the generated xs code:
if
(PDL->Version != PDL_CORE_VERSION)
Perl_croak(aTHX_
"Foo needs to be recompiled against the newly installed PDL"
);
If you add a new field to the F<Core> struct you should:
=over 5
=item *
discuss it on the pdl porters email list (pdl-devel
@lists
.sourceforge.net)
[
with
the possibility of making your changes to a separate
branch of the CVS tree
if
it's a change that will take
time
to complete]
=item *
increase by 1 the value of the C<
$pdl_core_version
> variable in
F<pdlcore.h.PL>. This sets the value of the
C<PDL_CORE_VERSION> C macro used to populate the Version field
=item *
add documentation (e.g. to L<PDL::API>)
if
it's a
"useful"
function
for
external module writers (as well as
ensuring the code is as well documented as the rest of PDL
;)
=back
=head1 BUGS
This description is far from perfect. If you need more details
or something is still unclear please ask on the pdl-devel
mailing list (pdl-devel
@lists
.sourceforge.net).
=head1 AUTHOR
Copyright(C) 1997 Tuomas J. Lukka (lukka
@fas
.harvard.edu),
2000 Doug Burke (djburke
@cpan
.org), 2002 Christian Soeller & Doug Burke,
2013 Chris Marshall.
Redistribution in the same form is allowed but reprinting requires
a permission from the author.