Jon’s Project Blog

I developed many finite elements during my Ph.D. research. Many finite element programs are built upon a mix of C/C++/Fortran. I had felt more comfortable with C/C++ over Fortran, and I have become fairly familiar with the C/Fortran calling conventions for struct/common blocks and etc.

Now that I’ve graduated, I no longer have access to Matlab for rapid prototyping and rapid checking of solutions to differential equations that I may be interested in. I recently ran across Julia which is a newer language. I decided to give it a try as I’ve done with R, Octave, and a slew of other “Matlab” clones. Julia has some nice features such as parallelism that I have not touched. It’s definitely not a Matlab clone, but it shares many similarities, and has some degree of built-in parallelism.

Recently, I’ve been involved in another material modeling project, and I’ve decided to use Julia during this project as a “rapid prototyping” replacement for Matlab. I found a gnuplot interface that seems to work pretty well called Gaston that seems to work pretty well. However, ultimately my code must run in Fortran for this project where it is compiled into the finite element code as a user material routine. Therefore, I decided to try to use Julia’s C/Fortran calling interface.

To do this in Julia, you must compile your subroutines and functions into a shared library. To illustrate this, we will look at the following fortran code that will help us investigate the underlying structure of arrays, matrices, and 3rd-order tensors in Julia.

c matrix test file
 
	subroutine vec(n,vector)
 
	integer n
	real*8 vector(*)
	write(*,*) (vector(i), i=1,n)
 
	end
 
	subroutine mat(m,n, matrix)
 
	integer m,n
	real*8 matrix(m,n)
	write(*,*) ((matrix(i,j),j=1,n),i=1,m)
 
	end
 
	subroutine tten(m,n,o,matrix)
 
	integer m,n,o
	real*8 matrix(m,n,o)
	write(*,*) (((matrix(i,j,k),j=1,n),i=1,m),k=1,o)
 
	end

We can then compile the fortran code into a shared object by the following:

gfortran matrix.f -o matrix.o -shared -fPIC

Then in Julia we can paste in the following code:

matrix = dlopen("matrix.o")
vec_ = dlsym(matrix,:vec_)
mat_ = dlsym(matrix,:mat_)
tten_ = dlsym(matrix,:tten_)

a = linspace(1.,4.,4)
ccall(vec_,Void,(Ptr{Int32},Ptr{Float64}),&int32(4),a)
a = linspace(1.,4.,4)'
ccall(vec_,Void,(Ptr{Int32},Ptr{Float64}),&int32(4),a)

b = ones(Float64,4,1)
c = b * a
ccall(mat_,Void,(Ptr{Int32},Ptr{Int32},Ptr{Float64}),&int32(4),&int32(4),c)

d = ones(Float64,4,4,2)
d[:,:,1] = c
ccall(tten_,Void,(Ptr{Int32},Ptr{Int32},Ptr{Int32},Ptr{Float64}),&int32(4),&int32(4),&int32(2),d)

The conclusion is that the underlying representation in Julia matches the ordering you would expect in Fortran for arrays, matrices, and 3rd-order tensors. With a bit of syntactic sugar, one can mask all the pointer casting, and create faster C/Fortran versions of the same Julia code. Hopefully, this languages provides a good way of rapid prototyping codes in a higher-level language, and then porting pieces to more traditional compiled object codes into high performance computing frameworks.

Julia isn’t the easiest language to use at the moment, but it seems to hold much promise. I had to write my ODE solver, which I guess is not too difficult. Now that I have this C/Fortran interfacing figured out, it should be trivial to add in ODE support from libraries from the national labs that have been well tested.