How to call C functions from PDL

You have a C function that performs some numerical computation and you want to call it from within PDL.

Assume you have a function like the following in a file called suma.c :

#include <math.h>
 
void suma(long int n,float *a, float* b, float *c) {
 
  long int j;
 
  for (j=0;j<stdio.h>

void suma(long int n,float *a, float* b, float *c);

int main() {

  float a[4];
  float b[4];
  float c[4];
  long int n;
  long int j;

  a[0]=1.;
  a[1]=2.;
  a[2]=3.;
  a[3]=4.;

  b[0]=5.;
  b[1]=6.;
  b[2]=7.;
  b[3]=8.;

  n=4;

  suma(n,a,b,c);

  for (j=0;j<n;j++) {
     printf("%d : cos(%f) + %f = %f\n",j,a[j],b[j],c[j]);
  }

}

We can now compile and run the program to see the results:

calbet@lago:~/pdlpp>cc -o mainsuma mainsuma.c suma.c -lm
mainsuma.c:
suma.c:


calbet@lago:~/pdlpp>./main_suma
0 : cos(1.000000) + 5.000000 = 5.540302
1 : cos(2.000000) + 6.000000 = 5.583853
2 : cos(3.000000) + 7.000000 = 6.010007
3 : cos(4.000000) + 8.000000 = 7.346356

Now we want to make a PDL module called Suma that will allow us to access the suma function from within PDL. Note the important name change (uppercase) between the name of the C function and filename, suma and suma.c, and the name of the PDL module or package to include the C function, Suma. The easiest way to do this is to use the PDL::PP module. It is documented in PDL::PP documentation and in the PDL FAQs. We can generate the PDL::PP definition code in a file called Suma.pd,

$VERSION = '0.1';

#-----------------------------------------------------------------------------------

pp_def('suma',
                Pars => 'float a(n); float b(n);float [o]c(n);',
                GenericTypes => [F],
                Code => 'suma($SIZE(n),$P(a),$P(b),$P(c));'
        );

pp_done();

Most of the contents of the Suma.pd file should be self explanatory. Note the special syntax

 $SIZE(n)

which tells PDL::PP to use the size of the piddles as input to the suma C function.

Also worth noting is the

[o]c(n)

prependend to what we want to be the output parameters of the C function. This is important because it will allow us to call the suma function within PDL with an output parameter as will be shown below.

Now comes the tricky part. We want to compile everything together and make PDL be aware of our newly created suma function inside PDL. For this the best thing is to copy the Makefile.PL file shown here:

# Use this as a template for the Makefile.PL for
# any external PDL module.

use ExtUtils::MakeMaker;

use PDL::Core::Dev;
@pack = (["Suma.pd",Suma,PDL::Suma]);

%hash = pdlpp_stdargs(@::pack);
# $hash{'OPTIMIZE'}  = '-g'; # If you want to debug, uncomment this.
# $hash{INC} .= " -I/usr/local/include";  # uncomment as required
# $hash{LIBS}[0] .= " -L/usr/local/lib -lmylib "; # uncomment as required

for $file ( qw( suma ) ) {
   my $n = "$file\$(OBJ_EXT)";
   $hash{OBJECT} .= " $n";
   $hash{clean}{FILES} .= " $n";
}
@prereq = (
          PDL       => "2.4.0",
         );
WriteMakefile(
       %hash,
       'NAME'  => 'PDL::Suma',
       'VERSION'       => 0.01,
       'PREREQ_PM'     => {@prereq}, # e.g., Module::Name => 1.1
     ($] ge '5.005') ? (
                        'AUTHOR' => 'My name , ([email protected])',
                        ) : (),
       );

# Add genpp rule
# add other makefile additions as required (see also ExtUtils::MakeMaker)
sub MY::postamble {
       pdlpp_postamble(@::pack);
}

You should modify this file by listing all your C files separated by spaces in the

 for $file ( qw( suma ) ) {

line. Note that the files are listed without the c extension. You should also change the name of the current module, Suma, to the name you want to give to your newly created PDL module all throughout your Makefile.PL. This Makefile.PL is designed to work automatically, if it does not work for you, try with a more manual approach of Makefile.PL as the one in How to call C and FORTRAN functions from PDL.

We are now ready to compile the module, but it is better to finish up one more thing before doing that in order for the Makefile generator to work properly. To call the function from within PDL we have to run perl/PDL allowing it access to the module. The easiest way to do this it to write a short program called test.pl in the following way,

# Before `make install' is performed this script should be runnable with
# `make test'. After `make install' it should work as `perl test.pl'

######################### We start with some black magic to print on failure.

# Change 1..1 below to 1..last_test_to_print .
# (It may become useful if the test is moved to ./t subdirectory.)

BEGIN { $| = 1; print "1..1\n"; }
END {print "not ok 1\n" unless $loaded;}
use PDL;
use PDL::Suma;

$loaded=1;

######################### End of black magic.

# Insert your test code below (better if it prints "ok 13"
# (correspondingly "not ok 13") depending on the success of chunk 13
# of the test code):


$a=pdl([1,2,3,4]);
$b=pdl([5,6,7,8]);

$c=suma($a,$b);

print $c,"\n";

print "ok 1\n";

You should substitute the last few lines of code by whatever you feel appropriate to test your subroutine inside PDL.

We are now really ready to compile the module by typing,

calbet@lago:~/pdlpp>perl Makefile.PL
Writing Makefile for PDL::Suma

and

calbet@lago:~/pdlpp>make
/opt/SUNWspro/bin/cc -c   -I/homespace/calbet/perl/lib/site_perl/5.8.6/sun4-solaris/PDL/Core  -I/usr/local/include -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64 -O    -DVERSION=\"0.1\" -DXS_VERSION=\"0.1\" -KPIC  "-I/homespace/calbet/perl/lib/5.8.6/sun4-solaris/CORE"   Suma.c
"Suma.xs", line 237: warning: implicit function declaration: suma
/opt/SUNWspro/bin/cc -c -o suma.o -g suma.c
Running Mkbootstrap for PDL::Suma ()
chmod 644 Suma.bs
rm -f blib/arch/auto/PDL/Suma/Suma.so
LD_RUN_PATH="" /opt/SUNWspro/bin/cc  -G -L/usr/lib -L/usr/ccs/lib  -L/opt/SUNWspro/prod/lib -L/usr/local/lib -L/usr/ucblib Suma.o suma.o  -o   blib/arch/auto/PDL/Suma/Suma.so   -lm 
chmod 755 blib/arch/auto/PDL/Suma/Suma.so
cp Suma.bs blib/arch/auto/PDL/Suma/Suma.bs
chmod 644 blib/arch/auto/PDL/Suma/Suma.bs
Manifying blib/man3/PDL::Suma.3

We are now ready to test our newly generated module!!!!

calbet@lago:~/pdlpp>make test
PERL_DL_NONLAZY=1 /homespace/calbet/perl/bin/perl "-Iblib/lib" "-Iblib/arch"   test.pl 
1..1
[ 5.5403 5.58385 6.01001 7.34636]
ok 1

You can verify that the output is the same as the one above.

Now you are ready to install this module in your PDL installation on your computer by typing (obtaining root privileges if necessary),

calbet@lago:~/pdlpp>su
Password:
root@lago:~/pdlpp#make install
rm -f blib/arch/auto/PDL/Suma/Suma.so
LD_RUN_PATH="" /opt/SUNWspro/bin/cc  -G -L/usr/lib -L/usr/ccs/lib  -L/opt/SUNWspro/prod/lib -L/usr/local/lib -L/usr/ucblib Suma.o suma.o  -o  blib/arch/auto/PDL/Suma/Suma.so   -lm
chmod 755 blib/arch/auto/PDL/Suma/Suma.so
Installing  /homespace/calbet/perl/lib/site_perl/5.8.6/sun4-solaris/auto/PDL/Suma/Suma.so
Installing /homespace/calbet/perl/lib/site_perl/5.8.6/sun4-solaris/auto/PDL/Suma/Suma.bs
Files found in blib/arch: installing files in blib/lib into architecture  dependent library tree
Installing /homespace/calbet/perl/lib/site_perl/5.8.6/sun4-solaris/PDL/Suma.pm
Installing /homespace/calbet/perl/man/man3/PDL::Suma.3
Writing /homespace/calbet/perl/lib/site_perl/5.8.6/sun4-solaris/auto/PDL/Suma/.packlist
Appending installation info to /homespace/calbet/perl/lib/5.8.6/sun4-solaris/perllocal.pod

You will now be able to access your module and function from a regular PDL program like the following called installtest.pl,

#!/usr/bin/env perl

use PDL;
use PDL::Suma;

$a=pdl([1,2,3,4]);
$b=pdl([5,6,7,8]);
 
$c=suma($a,$b);
 
print $c,"\n";

If you run it, you get again the same results,

calbet@lago:~/pdlpp>./installtest.pl
[ 5.5403 5.58385 6.01001 7.34636]

Congratulations!!!