How to call C and FORTRAN functions from PDL - PDLPorters/pdl GitHub Wiki
You have a FORTRAN subroutines and, possibly, also a C function that perform some numerical computation and you want to call them from within PDL.
Assume you have a FORTRAN F77 and a FORTRAN 90 subroutine like the following in a file called producto.f :
subroutine producto(n,a,b,c)
implicit none
integer n,j
real a(*)
real b(*)
real c(*)
do j=1,n
c(j)=a(j)*sin(b(j))
enddo
end
and prodexp.f90 :
subroutine prodexp(n, a, b, c)
implicit none
integer n,j
real a(n)
real b(n)
real c(n)
c=a*exp(b)
end subroutine prodexp
You also would like to include yet another function in the C language, 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);
extern int producto(long int * n,float *a, float *b, float *c);
extern int prodexp(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;
printf("Results for suma\n");
suma(n,a,b,c);
for (j=0;j<n;j++) {
printf("%d : %f * sin(%f) = %f\n",j,a[j],b[j],c[j]);
}
}
Note that at the beginning of the program we have declared the types of the C function and of the FORTRAN subroutines. The correspondence between FORTRAN and C variables are usually explained in the compiler manuals of our platform of choice. Normally they do not vary much from one architecture to another and the correspondence between FORTRAN and C types is relatively obvious, as seen in this example. Note also the call by reference of all the arguments in the C declaration of the FORTRAN subroutines ( all the variables are preceded by an asterisk). We should also prepend the word extern to the function declarations as in the example.
We can now compile and run the program to check that our assumptions on the C and FORTRAN variable correspondence are correct. We should also, in this step, be aware of the compilation switches we need for the subroutine or function naming convention. In the case of the SUN machine we would use:
calbet@lago:~/pdlpp>cc -c suma.c
calbet@lago:~/pdlpp>f77 -c -ext_names=plain producto.f
NOTICE: Invoking /opt/SUNWspro/bin/f90 -f77 -ftrap=%none -c -ext_names=plain producto.f
producto.f:
producto:
calbet@lago:~/pdlpp>f90 -c -ext_names=plain prodexp.f90
Note the FORTRAN compilation switch -ext_names=plain which tells the compiler to leave the FORTRAN subroutine name exactly as they are without adding any extra underscores to their names. In order to avoid future confusion, it is also normally preferable to always use lowercase for the subroutine names since the FORTRAN compiler usually translates them to lowercase internally (when called from C if you prefer).
Now we can compile and link the main program. On a SUN machine:
calbet@lago:~/pdlpp>cc -c mainsumprod.c
calbet@lago:~/pdlpp>f90 -o mainsumprod mainsumprod.o suma.o producto.o prodexp.o
Finally we can run the program to verify that everything went smoothly
calbet@lago:~/pdlpp>./mainsumprod
Results for 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
Results for producto
0 : 1.000000 * sin(5.000000) = -0.958924
1 : 2.000000 * sin(6.000000) = -0.558831
2 : 3.000000 * sin(7.000000) = 1.970960
3 : 4.000000 * sin(8.000000) = 3.957433
Results for prodexp
0 : 1.000000 * sin(5.000000) = 148.413162
1 : 2.000000 * sin(6.000000) = 806.857605
2 : 3.000000 * sin(7.000000) = 3289.899414
3 : 4.000000 * sin(8.000000) = 11923.832031
On a Linux PC we would compile and link like this:
calbet@laguna:~/pdlpp> gcc -c suma.c
calbet@laguna:~/pdlpp> g77 -c -fno-underscoring producto.f
calbet@laguna:~/pdlpp> gfortran-4.0 -c -fno-underscoring prodexp.f90
Note again the -fno-underscoring compiling switches to avoid the system naming the subroutines with additional underscores.
We can now compile and link the main program:
calbet@laguna:~/pdlpp> gcc -o mainsumprod mainsumprod.c suma.o producto.o prodexp.o -lm -lg2c
Note the need to link with two libraries. One of them is the usual one needed to use the math C libraries (-lm) and the other one to link FORTRAN and C code together (-lg2c).
We can now run the program and we should get
calbet@laguna:~/pdlpp> ./mainsumprod
Results for 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
Results for producto
0 : 1.000000 * sin(5.000000) = -0.958924
1 : 2.000000 * sin(6.000000) = -0.558831
2 : 3.000000 * sin(7.000000) = 1.970960
3 : 4.000000 * sin(8.000000) = 3.957433
Results for prodexp
0 : 1.000000 * sin(5.000000) = 148.413162
1 : 2.000000 * sin(6.000000) = 806.857605
2 : 3.000000 * sin(7.000000) = 3289.899414
3 : 4.000000 * sin(8.000000) = 11923.832031
Now we want to make a PDL module called Sumprod that will allow us to access the suma C function, the producto and the prodexp FROTRAN subroutines from within PDL. Note the important name change between the name of the C and FORTRAN functions and filenames (suma, producto, and prodexp) and the name of the PDL module or package to include the C function, Sumprod. The easiest way to do this is to use the PDL::PP module. It is documented PDL::PP documentation and in the PDL FAQs. We can generate the PDL::PP definition code in a file called Sumprod.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_def('producto',
Pars => 'float a(n); float b(n);float [o]c(n);',
GenericTypes => [F],
Code => 'producto(&$SIZE(n),$P(a),$P(b),$P(c));'
);
pp_def('prodexp',
Pars => 'float a(n); float b(n);float [o]c(n);',
GenericTypes => [F],
Code => 'prodexp(&$SIZE(n),$P(a),$P(b),$P(c));'
);
pp_done();
Most of the contents of the Sumprod.pd file should be self explanatory. Note the special sintax
$SIZE(n)
which tells PDL::PP to use the size of the piddles as input to the suma C function. Also note that in the FORTRAN functions we need to preprend an the ampersand reference operator to the size parameter
&$SIZE(n)
The reason for this is because in FORTRAN all calls are by reference.
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, producto and prodexp functions 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:
# Makefile.PL for a package defined by PP code.
use PDL::Core::Dev; # Pick up development utilities
use ExtUtils::MakeMaker;
use Config;
use ExtUtils::F77;
# -- Add new subroutines here! --
my @src = qw(suma.c producto.f prodexp.f90);
# -- Try not to mess with below stuff --
my $compiler_available = ExtUtils::F77->testcompiler;
my $myf77 = ExtUtils::F77->compiler();
my $underscore = _;_
my $trail_ = ExtUtils::F77->trail_; # fortran names need trailing underscore
if (($myf77 eq 'g77') and $trail_) { $underscore = '-fno-underscoring'; }
if (($myf77 eq 'f77') and $trail_ and ($Config{'osname'} eq 'dec_osf')) { $underscore = '-assume nounderscore'; }
print "Underscoring define = $underscore\n";
# Might want to add 'no underscore' rules for other architectures below...
my $runtime = ExtUtils::F77->runtime;
my @obj = @src;
map {s/\.f90/\.o/;} @obj; # swap .f90 for .o
map {s/\.[fc]/\.o/;} @obj; # swap .f, .c for .o
WriteMakefile(
'NAME' => 'PDL::Sumprod',
'VERSION_FROM' => 'Sumprod.pd',
'TYPEMAPS' => [&PDL_TYPEMAP()],
'OBJECT' => 'Sumprod.o ' . join (" ", @obj),
'DEFINE' => $define,
# 'PM' => { 'Sumprod.pm' => '$(INST_LIBDIR)/Sumprod.pm'},
'INC' => &PDL_INCLUDE(), # add include dirs as required by your lib
'LIBS' => "$runtime -lm", # add link directives as necessary
'clean' => {'FILES' => join (" ", @obj) },
'dist' => { COMPRESS => 'gzip', SUFFIX => 'gz' },
);
# Add genpp rule; this will invoke PDL::PP on our PP file
# the argument is an array reference where the array has three string elements:
# arg1: name of the source file that contains the PP code
# arg2: basename of the xs and pm files to be generated
# arg3: name of the package that is to be generated
sub MY::postamble {
$myfflags = ExtUtils::F77->cflags();
my $orig = pdlpp_postamble(["Sumprod.pd",Sumprod,PDL::Sumprod]);
my $added = _;_
foreach my $s (@src) {
my $o = $s;
$o =~ s/\.f90/\.o/;
$o =~ s/\.[cf]/\.o/;
if ($s =~ /\.c$/) { # C code
$added .= "$o: $s\n\t$Config{'cc'} -c -o $o -g $s\n\n";
} elsif ( $s =~ /\.f90$/) { # FORTRAN 90 code
$added .= "$o: $s\n\tf90 -ext_names=plain -c -o $o $s\n\n";
} else { # FORTRAN code
$added .= "$o: $s\n\t$myf77 -ext_names=plain -c -o $o -g $myfflags $underscore $s\n\n";
}
}
return $orig . $added;
}
This Makefile.PL file would be used in the SUN machine. Note that we can modify manually the compilation behaviour by changing explicitly the compilation flags and the linking ligraries. We can even make the modifications in the compilation line as we have done here
$added .= "$o: $s\n\t$myf77 -c -o $o -g -ext_names=plain $myfflags $underscore $s\n\n";
You should modify this file by listing all your C files separated by spaces in the
my @src=qw(suma.c)
line. You should also change the name of the current module, Sumprod, to the name you want to give to your newly created PDL module all throughout your Makefile.PL. Add also any necessary linking options in,
'LIBS' => "$runtime -lm", # add link directives as necessary
or in the compilation line as was shown above.
For example, for the PC Linux platform you would use the following Makefile.PL
# Makefile.PL for a package defined by PP code.
use PDL::Core::Dev; # Pick up development utilities
use ExtUtils::MakeMaker;
use Config;
use ExtUtils::F77;
# -- Add new subroutines here! --
my @src = qw(suma.c producto.f prodexp.f90);
# -- Try not to mess with below stuff --
my $compiler_available = ExtUtils::F77->testcompiler;
my $myf77 = ExtUtils::F77->compiler();
my $underscore = _;_
my $trail_ = ExtUtils::F77->trail_; # fortran names need trailing underscore
if (($myf77 eq 'g77') and $trail_) { $underscore = '-fno-underscoring'; }
if (($myf77 eq 'f77') and $trail_ and ($Config{'osname'} eq 'dec_osf')) { $underscore = '-assume nounderscore'; }
print "Underscoring define = $underscore\n";
# Might want to add 'no underscore' rules for other architectures below...
my $runtime = ExtUtils::F77->runtime;
my @obj = @src;
map {s/\.f90/\.o/;} @obj; # swap .f90 for .o
map {s/\.[fc]/\.o/;} @obj; # swap .f, .c for .o
WriteMakefile(
'NAME' => 'PDL::Sumprod',
'VERSION_FROM' => 'Sumprod.pd',
'TYPEMAPS' => [&PDL_TYPEMAP()],
'OBJECT' => 'Sumprod.o ' . join (" ", @obj),
'DEFINE' => $define,
# 'PM' => { 'Sumprod.pm' => '$(INST_LIBDIR)/Sumprod.pm'},
'INC' => &PDL_INCLUDE(), # add include dirs as required by your lib
'LIBS' => "$runtime -lm", # add link directives as necessary
'clean' => {'FILES' => join (" ", @obj) },
'dist' => { COMPRESS => 'gzip', SUFFIX => 'gz' },
);
# Add genpp rule; this will invoke PDL::PP on our PP file
# the argument is an array reference where the array has three string elements:
# arg1: name of the source file that contains the PP code
# arg2: basename of the xs and pm files to be generated
# arg3: name of the package that is to be generated
sub MY::postamble {
$myfflags = ExtUtils::F77->cflags();
my $orig = pdlpp_postamble(["Sumprod.pd",Sumprod,PDL::Sumprod]);
my $added = _;_
foreach my $s (@src) {
my $o = $s;
$o =~ s/\.f90/\.o/;
$o =~ s/\.[cf]/\.o/;
if ($s =~ /\.c$/) { # C code
$added .= "$o: $s\n\t$Config{'cc'} -c -o $o -g $s\n\n";
} elsif ( $s =~ /\.f90$/) { # FORTRAN 90 code
$added .= "$o: $s\n\tgfortran-4.0 -c -o $o -g $myfflags $underscore $s\n\n";
} else { # FORTRAN code
$added .= "$o: $s\n\t$myf77 -c -o $o -g $myfflags $underscore $s\n\n";
}
}
return $orig . $added;
}
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..3\n"; }
END {print "not ok 1\n" unless $loaded;}
use PDL;
use PDL::Sumprod;
$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";
$c=producto($a,$b);
print $c,"\n";
print "ok 2\n";
$c=prodexp($a,$b);
print $c,"\n";
print "ok 3\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,
xcalbet@laguna:~/pdlpp$ perl Makefile.PL
ExtUtils::F77: Version 1.15
Loaded ExtUtils::F77 version 1.15
Found compiler g77
ExtUtils::F77: Using system=Linux compiler=G77
Checking for gcc in disguise:
Compiler is gcc version 4.1.2 20060729 (prerelease) (Debian 4.1.1-10)
Runtime: -L/usr/lib/gcc/i486-linux-gnu/3.4.6/../../../../lib -L/usr/lib -lg2c -lm -L/usr/lib/gcc/i486-linux-gnu/4.1.2 -lgcc
ExtUtils::F77: Validating -L/usr/lib/gcc/i486-linux-gnu/3.4.6/../../../../lib -L/usr/lib -lg2c -lm -L/usr/lib/gcc/i486-linux-gnu/4.1.2 -lgcc [ok]
ExtUtils::F77: Compiler: g77
ExtUtils::F77: Cflags: -O
Compiling the test Fortran program...
Executing the test program...
Congratulations you seem to have a working f77!
Underscoring define = -fno-underscoring
Writing Makefile for PDL::Sumprod
and
xcalbet@laguna:~/pdlpp$ make
cp Suma.pm blib/lib/PDL/Suma.pm
cp Sumprod.pm blib/lib/PDL/Sumprod.pm
/usr/bin/perl /usr/share/perl/5.8/ExtUtils/xsubpp -typemap /usr/share/perl/5.8/ExtUtils/typemap -typemap /usr/local/lib/perl/5.8.4/PDL/Core/typemap.pdl Sumprod.xs > Sumprod.xsc && mv Sumprod.xsc Sumprod.c
cc -c -I/usr/local/lib/perl/5.8.4/PDL/Core -D_REENTRANT -D_GNU_SOURCE -DTHREADS_HAVE_PIDS -DDEBIAN -fno-strict-aliasing -pipe -I/usr/local/include -D_LARGEFILE_SOURCE -D_FILE_OFFSET_BITS=64 -O2 -DVERSION=\"0.1\" -DXS_VERSION=\"0.1\" -fPIC "-I/usr/lib/perl/5.8/CORE" Sumprod.c
In file included from /usr/local/lib/perl/5.8.4/PDL/Core/pdlcore.h:10,
from Sumprod.xs:9:
/usr/local/lib/perl/5.8.4/PDL/Core/ppport.h:227:1: warning: "PERL_UNUSED_DECL" redefined
In file included from Sumprod.xs:6:
/usr/lib/perl/5.8/CORE/perl.h:163:1: warning: this is the location of the previous definition
cc -c -o suma.o -g suma.c
g77 -c -o producto.o -g -O -fno-underscoring producto.f
gfortran-4.0 -c -o prodexp.o -g -O -fno-underscoring prodexp.f90
Running Mkbootstrap for PDL::Sumprod ()
chmod 644 Sumprod.bs
rm -f blib/arch/auto/PDL/Sumprod/Sumprod.so
LD_RUN_PATH="/usr/lib/gcc/i486-linux-gnu/3.4.6/../../../../lib" cc -shared -L/usr/local/lib Sumprod.o suma.o producto.o prodexp.o -o blib/arch/auto/PDL/Sumprod/Sumprod.so \
-L/usr/lib/gcc/i486-linux-gnu/3.4.6/../../../../lib -L/usr/lib -lg2c -lm -L/usr/lib/gcc/i486-linux-gnu/4.1.2 -lgcc -lm \
chmod 755 blib/arch/auto/PDL/Sumprod/Sumprod.so
cp Sumprod.bs blib/arch/auto/PDL/Sumprod/Sumprod.bs
chmod 644 blib/arch/auto/PDL/Sumprod/Sumprod.bs
Manifying blib/man3/PDL::Suma.3pm
Manifying blib/man3/PDL::Sumprod.3pm
We are now ready to test our newly generated module!!!!
xcalbet@laguna:~/pdlpp$ make test
PERL_DL_NONLAZY=1 /usr/bin/perl "-Iblib/lib" "-Iblib/arch" test.pl
1..3
[ 5.5403 5.58385 6.01001 7.34636]
ok 1
[-0.958924 -0.558831 1.97096 3.95743]
ok 2
[148.413 806.858 3289.9 11923.8]
ok 3
You can verify that the ouput 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),
xcalbet@laguna:~/pdlpp$ su
Password:
laguna:/home/xcalbet/pdlpp# make install
Installing /usr/local/lib/perl/5.8.8/auto/PDL/Sumprod/Sumprod.so
Installing /usr/local/lib/perl/5.8.8/auto/PDL/Sumprod/Sumprod.bs
Files found in blib/arch: installing files in blib/lib into architecture dependent library tree
Installing /usr/local/lib/perl/5.8.8/PDL/Suma.pm
Installing /usr/local/lib/perl/5.8.8/PDL/Sumprod.pm
Installing /usr/local/man/man3/PDL::Suma.3pm
Installing /usr/local/man/man3/PDL::Sumprod.3pm
Writing /usr/local/lib/perl/5.8.8/auto/PDL/Sumprod/.packlist
Appending installation info to /usr/local/lib/perl/5.8.8/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::Sumprod;
$a=pdl([1,2,3,4]);
$b=pdl([5,6,7,8]);
$c=suma($a,$b);
print $c,"\n";
print "ok 1\n";
$c=producto($a,$b);
print $c,"\n";
print "ok 2\n";
$c=prodexp($a,$b);
print $c,"\n";
print "ok 3\n";
If you run it, you get again the same results,
xcalbet@laguna:~/pdlpp$ ./installtest.pl
[ 5.5403 5.58385 6.01001 7.34636]
ok 1
[-0.958924 -0.558831 1.97096 3.95743]
ok 2
[148.413 806.858 3289.9 11923.8]
ok 3
Congratulations!!!