22 April 2009

Clawpack Update - Inter-Language Communication I

Yesterday I finished up a nice toy application that does the following:

  1. compiles a Fortran dynamic library containing a subroutine that accepts a double and a function pointer and simply calls the pointed function with the double as an argument,

  2. creates a test C program that passes a double and simple squaring function to the Fortran library,

  3. compiles a Cython class containing a method that passes a double and a simple squaring function to the Fortran library,

  4. executes a Python script that imports the Cython class and executes various methods.



First, let's take a look at the Fortran library. Fortran 95's way of defining function pointers is a bit different from C's method. Here's the code:

! clawpack.f95
subroutine evolve(q, rp)
! define inputs
double precision, intent(inout) :: q
interface
subroutine rp(q)
double precision, intent(inout) :: q
end subroutine
end interface

! execute input function
call rp(q)
end subroutine


We compile this using the following commands (Mac OS X):

$ gfortran -g -c -fPIC -arch i386 -arch ppc clawpack.f95 -o clawpack.o
$ gcc -arch i386 -arch ppc -dynamiclib clawpack.o -o libclawpack.dylib

I'll discuss why I explicitly include the above architecture flags when we get to the Cython portion of this discussion. In order to test this functionality, I wrote a little C program. But first, for C to interface with the library, we need a header file. Taking a look at the library symbol table will show that gfortran prefixes each subroutine with an underscore "_". Therefore, the header should look a little something like this:


// clawpack.h
extern void evolve_(double *q,
void (*rp)(double *q));


Note that Fortran subroutines actually accepts pointers to their inputs rather than the contents at the addresses, themselves; at least this is so from the C perspective. Therefore, we pass in a pointer to the double precision variable. Additionally, the C equivalent of a Fortran subroutine is
a void function. This website provides some nice examples of C/Fortran equivalences and communication.

Now for the program itself:


// test.c
#include
#include
#include "clawpack.h"

// an example function that we will pass to evolve_
void c_square(double *q);
void c_square(double *q) {
*q = (*q)*(*q);
}

int main(int argc, char **argv) {
double *val;
val = (double *) malloc( sizeof(double) );
if (argc == 1)
*val = 0.0;
else
*val = atof( argv[1] ); // obtain val as command line arg.

// square once
c_square(val);
printf("val ** 2 = %f\n", *val);
// square again from Fortran library
evolve_(val, c_square);
printf("val ** 4 = %f\n", *val);

free(val);
return 0;
}


Compile against the library with the header and you'll get something like this:

$ ./test_prog 2.0
val ** 2 = 4.000000
val ** 4 = 16.000000
$ ./test_prog 5.0
val ** 2 = 25.000000
val ** 4 = 625.000000
$


This means that function pointer passing from C to Fortran works! In particular, the value of val is retained throughout the process - important for when we begin passing around important information like adaptive meshes and boundary value data.

There's a lot that goes on in the Cython module; both in it's implementation and compilation. Therefore, I'll reserve an entire post for all of that information. Expect to see it tomorrow!

Clawpack Update - Thoughts on 5.0 Design

After spending some time trying to create a Sage wrapper for the Clawpack package, I decided to dial it back a bit and talk to Kyle about Clawpack 5.0. The goal is to provide an easy to use Python interface - as opposed to the "copy and compile" structure of the current iteration - while maintaining the option to provide fast Fortran or C subroutines for use in the computationally heavy portions of the calculation. Enter Cython: a language for easily writing C extensions for Python.

The way we see it thus far, and this is definitely subject to change, there will be three primary components:

  • Fortran Library - the main Fortran library will contain the engine of the computations. In order to define the particular functionality, say for example when the user needs to provide an appropriate Riemann solver or grid for their particular problem, they pass in pointers to their own functions. Thus, the user can potentially solve complicated problems by writing the appropriate subroutines in a single Fortran or C file, pass pointers to these functions to the Clawpack library routines, and link and compile against that library.

  • Python Interface - much like Pyclaw, there will be a Python interface to all of this functionality. The user will have access to user-friendly classes, such as Solver, that will let them define and customize the behavior of the program so that they can solve their particular problem.

  • Cython Bridge - finally, Cython will act much like a bridge between the Python and Fortran. Cython provides a lot of functionality for communication between Python and compiled libraries, making it possible to create highly-optimized algorithms with a nice Python interface. Also, it will hopefully allow users to pass in functions written in C or Fortran in order to provide really fast computation.



In my next post, I'll describe what I've done thus far to get the languages talking to one another.

11 April 2009

Level Up!

Every year, whenever I attend the three-hour-long Easter Vigil Mass and am attentive throughout the entire service I always feel like I gain two levels in Catholicism at the end.

Sorry for the cheesy Dungeons and Dragons joke. Happy Easter!

09 April 2009

Clawpack - Design Strategy

Over the past two weeks I have been experimenting with different ways to incorporate Clawpack into Sage. It's been a challenge since, as I've discovered, the two systems have quite distinct design philosophies. Though, part of the issue is the language difference. Fortran doesn't have the built-in modularity provided by Python.

My investigation has led me to several design strategies so far. I share these not only so I can organize my thoughts but also to see if anyone has suggestions on how to proceed.

Strategy #1: Insert Pyclaw
Pyclaw is a Clawpack Python wrapper designed by Kyle Mandli. I've already partially documented my attempts at incorporating this code into Sage in a Google Notebook which you can view here. The idea is to simply connect the Pyclaw code to Sage as a collection of additional Python modules.
Pros:

  • Conceptually straightforward

  • Python to Fortran interface is already done


Cons:

  • Requires augmentation of high-level Pyclaw code

  • Need to rework plotting routines to use Sage's plotting functionality

  • Pyclaw is constantly changing



Strategy #2: Wrap the Core Fortran Code: Cython
Another approach is to create a Cython wrapper of the core Fortran code found in claw/clawpack. Then, the user can write Cython code that will interface with the library. However, since most of the Subroutines are stubs this would require some compiling against the library thus overriding the symbol table.
Pros:

  • Fast.

  • Allows the user to write the step1, src1, and other user-supplied Clawpack functions in Cython in addition to Fortran.


Cons:

  • Recompiling routines might be a bit tricky, especially if the user is working in the Sage Notebook.

  • Need to figure out a way to handle claw.data files.

  • What happens when you Cython extern and then change the library contents / symbol table? Needs to be examined.



Clawpack and Sage are two different beasts and it'll be a challenge to harmoniously combine the two. However, like I said before, this will be a great way to introduce Sage to numerics.