Post a Comment On: Charlie's Second Blog

I recently tried to write a numerical computation software that needed to perform a lot of vector operations. My first idea was to use python, and the numpy library.
But then I found that this simple operation :

y = a * x + y

Where x and y are large numpy arrays, and a is a scalar,
was much slower than the same thing coded in C, like this :

void func(float x[N], float y[N], float a) {
for(size_t i = 0; i < N; ++i) {
y[i] += a * x[i];
}
}

The reason is that numpy has to create a temporary object to store the result of the multiplication, and then also create a new object to store the result of the addition. When the size of the arrays becomes large - I was dealing with several Megabytes data arrays - it can make a lot of difference.

But I don't like the C code either, because of the very ugly loop.

So what can I do ? Well, first I could use the C blas library, then the operation is replaced by a function call, like this :

cblas_daxpy(N, a, (double*)x, 1, (double*)y, 1);

But it still looks over complicated.

OK maybe some other languages can be more suitable for this task. What about C++ ? With the boost::ublas library, it looks already much better :

using namespace boost::numeric::ublas;
void func(const vector& x, vector& y, double a) {
y += a * x;
}

Well, almost what I wanted, I could even use template parameters to make the function accepts any kind of input data.

But since at the end I want to be able to call this function from python, C++ is not the easiest choice (beside, I don't like C++ that much). Finally I used a language I really don't like, but that appeared to be the best for this particular task : Fortran90

Fortran90 is already much better than the horrible Fortran77 (That nobody should ever use). and most important, it has vector operations included in the language.

So now my simple operation looks like this :

subroutine func(x, y, a)
implicit none
REAL, dimension(:), intent(in) :: x
REAL, dimension(:), intent(inout) :: y
REAL, intent(in) :: a

y = a * x + y

end subroutine func

(Note : in Fortran90, you DON'T need to start every lines at the 7th character !)

I can call this fortran function from python quite easily, thanks to F2PY, and voilà. Then I can extend my function to perform all sort of numerical operations.

Conclusion : Fortran90 is not as bad as I though it was, combined with python it can even be quite useful. But I still don't like this language, I am looking forward for the D language to have good numerical operations library, as well as python wrapping (it already exists, but there are still a few issues.)