Porting Python Code to Fortran for Parallel Computing

[SevenThunders]

I have a simulation that involves a lot of dense linear algebra with a lot of complex arrays, that was written in python's numpy and then sped up with cython. Unfortunately I'm stuck with python's GIL (global interpreter lock), which prevents the use of numpy with multi-threaded code.

I'm looking at rewriting sizable chunks of this code into C, or C++ or Fortran. Fortran is actually looking pretty good right now given it's built in array support and support for parallel operations. The newer compilers support array slicing, operator overloading and even classes.

The global mindset in the computer science world is that Fortran is a dead language and that we should all be programming in Java or C++, or the latest toy language. However it's looking mighty attractive to me right now for the kind of high performance computing I need. So someone please talk me out of it LOL.

One thing though that surprised me is that I have not found any higher level wrappers that are maintained for lapack and blas operations in Fortran. These things abound in C++ and in many scripting languages.

 

[StoneTemplePython]

on the Python side of things, a couple thoughts.

One is that Cython generally doesn't play well with builtin numpy functions (a very simple example: numpy's dot products vs writing your own for loop) as each time you call a numpy function/method, it's a slow call to the interpreter in Cython. You may have speed issues long before threading comes up as an issue.

Two, while a lot less mature, Numba was made exactly for these situations. It plays extremely well with numpy (but not much else). In general if you use the "nogil" feature, you can bypass GIL. Here's a link that may be of interest:

https://roman-kh.github.io/numba-2/- - - - -

The global mindset in the computer science world is that Fortran is a dead language and that we should all be programming in Java or C++, or the latest toy language. However it's looking mighty attractive to me right now for the kind of high performance computing I need. So someone please talk me out of it LOL.

If you are extremely concerned about performance I don't think anyone would suggest Java. For example: just about every granular thing for Machine Learning is done in C++ these days. Fortran is fine and perhaps simpler than C++, but C++ is the standard in Machine Learning world.
- - - - -
Personally, I'd probably reach for Julia in situations like this, but a lot people would concur with C++ or Fortran. Good luck.

 

[SevenThunders]

Super appreciate your reply StoneTemplePython. Over the years I have become leery of complete rewrites, so I'm somewhat partial to my cython code since I have reasonable confidence it is correct. Although it would be reasonable to consider converting it back to python and then optimizing hot spots in C++ or Fortran if that makes sense.

While I have done work in machine learning, this application is all about iteratively solving hundreds of medium sized, complex, dense linear least squares problems. I have a love hate relationship with C++, since I find you are always wrestling with the compiler, and I just want to get results. As such I try to stick to a reasonable subset of the language (e.g., I avoid complex class templates and try to exploit return value optimization.)

Julia is probably the go-to scripting language of the future for scientific programming, but this is legacy code in python/cython and there will be no rewrite of that in the near future, especially since Julia is a moving target. I was just wondering if any Fortran users could weigh in on this. I'd hate to start in and then have to change course because the language gets in the way.

Finally I'll take another look at Numba (it's been years since I considered it), but I suspect the legacy issue might preclude it.

 

[Dr Transport]

...The global mindset in the computer science world is that Fortran is a dead language and that we should all be programming in Java or C++, or the latest toy language. However it's looking mighty attractive to me right now for the kind of high performance computing I need. So someone please talk me out of it LOL..

Fortran isn't dead, I know of at least 4 major us govt programs written in it and there are no plans to re-write the codes in any newer language. I also know that there are multiple academic codes that won't die and the aerospace industry has numerous codes in it also which won't die anytime soon.

I thought I left it behind when I got my degree, but have programmed in it off and on for the past 15 years and my new job has all the codes associated with it written in a combination of Fortran 90/95/98/2003 etc...This is in addition to awk and sed...

it won't hurt you to learn it if you already haven't...

 

[SevenThunders]

Fortran isn't dead, I know of at least 4 major us govt programs written in it and there are no plans to re-write the codes in any newer language. I also know that there are multiple academic codes that won't die and the aerospace industry has numerous codes in it also which won't die anytime soon.

I thought I left it behind when I got my degree, but have programmed in it off and on for the past 15 years and my new job has all the codes associated with it written in a combination of Fortran 90/95/98/2003 etc...This is in addition to awk and sed...

it won't hurt you to learn it if you already haven't...

It all depends on what you call "learn". Like many engineers/scientists back in the day, Fortran was the language they taught you in your intro programming course. It was not a pleasant experience, but punchcards, limited computer access, and a hostile, student hating professor were part of that process.

I'm not worried about "learning" the language. I just want to make sure that I can be productive in it and that I can get the coarse grained parallel performance I'm looking for.

If you've been using it, I'm curious about whether that was productive or not and I also wonder what you do for linear algebra (e.g. least squares problems, matrix inversions, eigenvalue problems etc.). Do you just call raw blas or lapack libraries? Another issue I've had with the linear algebra is that libraries like openblas are multithreaded, but possibly have been throttling my coarse grained parallelization efforts due to some odd effects when I call blas routines from different threads.

 

[jedishrfu]

If you're thinking of a rewrite you should also look into using Julia over Fortran.

While Fortran has evolved over the decades its still Fortran and isn't using the most cutting strategies for computation.

For that, you need to consider other languages. Julia is a rather young language that is making inroads into computational modeling. It works well with Fortran, C/C++, R, and Python with a syntax similar to Matlab.

https://julialang.org/

 

[SevenThunders]

If you're thinking of a rewrite you should also look into using Julia over Fortran.

While Fortran has evolved over the decades its still Fortran and isn't using the most cutting strategies for computation.

For that, you need to consider other languages. Julia is a rather young language that is making inroads into computational modeling. It works well with Fortran, C/C++, R, and Python with a syntax similar to Matlab.

https://julialang.org/

Yes I've already completed a good sized project in Julia and it was pretty nice. Fortran is faster though isn't it? Why can't Fortran be cutting edge? It seems that Fortran 2008 incorporates a bunch of modern idioms. Maybe the only thing they aren't doing is making functions first class objects, but they do have class objects, a type system, and at least their operator overloading is polymorphic.

I haven't written anything serious in Fortran yet, but I think my biggest complaint so far is all the redundant declarations. At any rate I would need something I can easily call from Python, which are easily satisfied by C++ and Fortran. There is the PyCall.jl module and I think there is a module that does allow julia calls within python, but I'm a bit cautious for production code that will be used by my client for the long haul.

 

[Dr Transport]

I can't say that my programs were not productive, each and every project has its pluses and minuses. As for linear algebra, etc, if you need Lapack or whatever, that is what you use. I used IMSL for a long time and was happy with the results. I didn't do any multi-threading for the things I did, it wasn't necessary.

 

[FactChecker]

Yes, FORTRAN can be faster, even than C (this surprised me). In a speed comparison, Fortran and C are the fastest and some other languages are shockingly slow, depending on which features are used. If you are converting a lot of matrix code (or any structures with > 1 dimension), you should consider the difficulty of converting the Python index order (row-major) into FORTRAN index order (column-major), and the Python index-base 0 to the FORTRAN index-base 1. Python also allows negative indices, which would be a problem to convert if it occurs a lot in your code.

Other than those issues, I would be very surprised if most languages did not have libraries that do what you need.

 

[jim mcnamara]

If you are running on a desktop or have cloud access consider pyCUDA, which allows you to parallelize python code on NVIDIA GPU's or on a cloud service, many provide NVIDIA GPU's.

See: https://developer.nvidia.com/pycuda Free to download.

 

[SevenThunders]

If you are running on a desktop or have cloud access consider pyCUDA, which allows you to parallelize python code on NVIDIA GPU's or on a cloud service, many provide NVIDIA GPU's.

See: https://developer.nvidia.com/pycuda Free to download.

Back in the day I was using cython to hook into Nvidia's CUDA. It's great for fine grained parallelization, and thus works best for the same code replicated over a lot of data. I might be able to shoehorn my problem into that, but the lowest hanging fruit is to simply be able to do something like a simple parallel for loop, with each inner loop process performing some least squares iterations.

I've used OpenMP and Matlabs parallel library to good effect for this sort of thing. I never thought that numpy would preclude this. But numpy depends on Python's GIL, and there is no easy way around it. Even the Numba documentation alludes to this. This creates the problem of replacing all the numpy code inside that inner loop with something that can be parallelized. Hopefully this is where Fortran can come to play I think.

 

[newjerseyrunner]

Doesn't python have the ability to write extensions for it? I write a lot of server code in PHP, which is a single threaded language, but it allows me to create "phpized" libraries in C/C++ and interface with them directly from the PHP, I'm sure python has something similar.

 

[harborsparrow]

One caveat for porting numerical code "with a lot of arrays" is to look into how each language stores the arrays internally. My husband was porting a lot of FORTRAN code into C and ran into this. One language stores the rows across and the columns down, or vice versa, or something to that effect. And it can matter, sometimes, in obtaining an accurate port.

The other is, that IMO concerns about performance in Java are no longer reasonable. Java is pretty fast and optimizing these days. But back in the day, it was a concern.

 

[FactChecker]

One caveat for porting numerical code "with a lot of arrays" is to look into how each language stores the arrays internally. My husband was porting a lot of FORTRAN code into C and ran into this. One language stores the rows across and the columns down, or vice versa, or something to that effect. And it can matter, sometimes, in obtaining an accurate port.

The other is, that IMO concerns about performance in Java are no longer reasonable. Java is pretty fast and optimizing these days. But back in the day, it was a concern.

I agree with your warning about the row-priority versus column priority conversion. But I am not so sure about the speed comment. The OP was concerned with speed. When speed is the main issue, there can be a lot of difference between languages.

 

[SevenThunders]

After a week or two writing and speed profiling some Fortran code, I can honestly say that I'm extremely happy with it's performance and it's expressive power. As a matter of fact my fairly simple least squares solver appears to be competitive with Lapack's gels subroutine. Especially nice for my application is the ability to use array slicing notation.

This is exemplified in this little iterative subroutine I wrote to actually apply a QR factorization to an overdetermined rectangular complex matrix,

! generate a QR decomposition implicitly using Householder xforms
pure subroutine qhr(A,b)
    implicit none
    complex*16, intent(inout) :: A(:,:)
    complex*16, intent(inout) :: b(:,:)
! local temporaries,  first a vector for a householder reflector
    integer :: mrows, ncols
    complex*16 hv(size(A,1),1)
    integer k
    mrows = size(A,1)
    ncols = size(A,2)
! loop over each sub-matrix starting along diagonal
    do k=1,ncols
        hv(k:mrows,:) = household(A(k:mrows,k))
! apply householder reflector to A
         call applyhouse(hv(k:mrows,:), A(k:mrows,k:ncols))
! apply householder reflector to b
        call applyhouse(hv(k:mrows,:), b(k:mrows,:))
    end do

end subroutine qhr

Fortran gets a few things right here. One is that vanilla arrays appear to be flat contiguous blobs of memory. Also the array slicing designates strides for a subarray, which can be written to, to modify the original array. This reduces copying overhead. Fortran appears to support recursive subroutines, which I assume are tail-call-optimized.

With native support for classes, user defined types (soon to support parametric types), and intrinsics such matrix transpose and matrix multiplication, there is a lot to like for applications that need to do a lot of linear algebra. Having programmed in virtually every language out there, from Haskell to C++, D, Julia, python and java, I'd say I could be more productive in Fortran, than any other low level language, for general scientific programming, in particular I'd rather use Fortran than C++ for numerics.

I'm not saying here that, Python or Julia or Matlab is less productive than Fortran, quite the contrary, but if you need the performance, Fortran is looking pretty good to me. The only catch so far has been that the built in parallel do constructs such as do concurrent, or for all don't actually seem to spawn multiple threads in gfortran. I suspect the support is better with Intel's ifortran. Also fortran automatically lowers the rank (aka dimension) of an array if you slice it with a single value. This forces you to program all the edge cases wherein you might have different dimension arrays as inputs to a subroutine.

I'm waiting for functions as first class objects, better support for co-arrays and map reduce constructs run concurrently over multiple threads to complete the trifecta. Nevertheless I'm sufficiently impressed to continue fully with rewriting my core cython code into fortran, to gain added speed.