# Parallel Scientific Computation Software Development Language?

I want to develop a parallel scientific computation software from scratch. I want some thoughts on which language to start. The program involves reading/writing data to txt files and doing heavy computations in parallel, with many LU factorizations and the use of sparse linear solvers. The candidate solutions I was thinking are Fortran 2003/2008 with OpenMP or co-array, C++ with openmp cilk+ or TBB, python. Any other, documented, suggestions are welcome! I know very well C, Fortran and Java (in that order). I've done some scripting in python but basic stuff.

I know fortran is very fast, but, hard to maintain and parallelize. C++ is said to be slow unless you use external libraries etc Python I like, but is it realistic to write a full scale, industrial level software upon?

The software needs to be able to handle big amounts of data and be effective with scientific computations. The performance is of the essence.

For the background, I already have a working software written in Fortran. Many people were involved in development over many years and the code is really dirty. Maintaining and parallelizing the code has proved a nightmare and I'm thinking of alternatives.

• As a C++ wonk, I wouldn't call Fortran hard to maintain. Maintainability is tied to good practices for the most part, not language choice. The slowness of C++ is oversold. Also, I would recommend that you augment this post to describe your data size and turnaround time requirements. I've seen "big" vary by 9 or 10 orders of magnitude depending on who I'm talking to. Jun 10 '12 at 15:38
• @BillBarth is too modest in selling his C++ skills, but he's also too generous in his statement that the "slowness of C++ is oversold". There have been a number of C++ vs Fortran threads in scicomp.stackexchange.com that have discussed this very question and the general conclusion was that it's simply not true any more than C++ is slower than Fortran for almost all cases. I personally think that today it could be considered an urban myth. What is very much true is that if you take into account maintainability of the code, then Fortran doesn't fare very well today. Jun 10 '12 at 22:48
• @BillBarth I find Fortran unmaintainable due to cumbersome handling of pointers and function pointers. Also, it starts at 1. ;-) Jun 10 '12 at 23:59
• @BillBarth and others, if you'd like to continue to discuss the general merits of the Fortran, C++, and other languages, please take it to the scicomp chat room and @ anybody you'd like to specifically address. Jun 11 '12 at 6:26
• @AronAhmadia: ah, come on, I have so much to say to Jed ;-) (Jed: some other time. In our case, no STL for sparse matrices, but lots of it in the adaptive mesh data structures.) Jun 11 '12 at 15:59

Let me try and break down your requirements:

• Maintainability
• Strong interfaces/capability for LU factorizations
• Sparse linear solvers
• Performance and scalability to large data

From this list, I would consider the following languages:

C, C++, Fortran, Python, MATLAB, Java

Julia is a promising new language, but the community is still forming around it and it has not been deployed in any major new codes.

This is easy to get right in any programming language. Make sure you are appropriately buffering and coalescing your I/O access, and you will get good performance from any of the languages you should consider. Avoid the stream objects in C++ unless you know how to use them performantly.

### Strong interfaces/capability for LU factorizations

If you are performing dense LU factorizations, you will want to use LAPACK, or ScaLAPACK/Elemental for parallel functionality. LAPACK and ScaLAPACK are written in Fortran, Elemental is written in C++. All three libraries are performant and well-supported and documented. You can interface into them from any of the languages you should consider.

### Sparse linear solvers

The premier freely available sparse linear solvers are almost all available through PETSc, written in C, which is well-documented and supported. You can interface into PETSc from any of the languages you should consider.

### Performance and scalability to large data

The only parallel programming paradigms you mention are shared memory based, which means you are not considering an MPI-based (message-passing), distributed-memory computing approach. In my experience, it is much easier to write code that scales well beyond a dozen cores using a distributed-memory solution. Almost all University "clusters" are MPI-based these days, large shared-memory machines are expensive, and correspondingly rare. You should consider MPI for your approach, but my advice will apply regardless of the programming paradigm you choose.

With regards to on-node performance, if you are writing numerical routines yourself, it is easiest to get good serial performance in Fortran. If you have a little bit of experience in C, C++, or Python, you can get very comparable performance (C and C++ are dead-even with Fortran, Python and MATLAB come within about a 25% time overhead without much effort). MATLAB does this through a JIT compiler and very good linear algebra expressivity. You will likely need to use either Cython, numpy, numexpr, or embed numerical kernels to get the claimed performance from Python. I can't comment on Java's performance, because I don't know the language very well, but I suspect it is not far from Python's if written by an expert.

### A note on interfaces

I hope I've convinced you that you are going to be able to do everything you want in any of the programming languages you are considering. If you are using Java, the C interfaces will be a little challenging. Python has excellent C and Fortran interface support through ctypes, Cython, and f2py. LAPACK is already wrapped and available through scipy. MATLAB has all of the functionality you need in its native libraries, but is not readily scalable or particularly easy to run on clusters. Java can support C and Fortran interfaces with the JNI, but is not commonly found on clusters and in parallel software for scientific computing.

### Maintainability

A lot of this is going to come down to personal flavor, but the general consensus on maintainability is that you want to minimize the number of lines of code in your software, write modular code with well-defined interfaces, and for computational software, provide tests that verify the correctness and functionality of the implementation.

### Recommendation

I personally have had a lot of luck with Python and I recommend it for many computational projects. I think you should strongly consider it for your project. Python and MATLAB are probably the most expressive of the languages available for scientific computing. You can easily interface Python to any other programming language, you can use f2py to wrap your current Fortran implementation and piece-by-piece rewrite whichever parts you wish in Python while verifying that you are maintaining functionality. At this time, I would recommend a combination of the official Python 2.7 implementation with scipy. You can get very easily started with this stack from the freely available Enthought Python Distribution.

You could also do most of this in C, C++, or Fortran. C and C++ are very appealing languages for professional developers with a lot of experience, but frequently trip new developers and are in this sense probably not a great idea for a more academic code. Fortran and MATLAB are popular in academic computation, but are weak at the advanced data structures and expressivity Python offers (think of a Python dict object, for example).

### Related Questions:

• A very well documented, all inclusive answer. Under Fortran I use a lot Lapack. I will take a look at python and try to wrap my Fortran code to begin with and slowly slowly move to Python. The only thing that scares me is the 25% time overhead I might have. But if it comes with the benefit of more expressive code and better parallel computing handling, I'll go for it. I mentioned shared memory only because the software currently runs in an interactive way (make a change on the data and rerun) on 2,4,8,24,48-core shared memory computers of researchers in the Uni under Windows and Linux. Jun 10 '12 at 18:23
• I don't know how you can make the claim of 25% overhead for numerical kernels written in Python. Pure Python numerical kernels are often order of 100x slower than C. Numpy and numexpr can do a decent job with certain expressions, but that's hardly writing new numerical kernels in Python. Cython can make some things fast, but usually not within 25% of C. Python is a fine "glue" language, but I think Aron is overselling it as a general purpose solution for performance-sensitive tasks. Jun 10 '12 at 20:31
• I/O is the weak point of Fortran, because Fortran requires a lot of structure in I/O. My second-hand experience from talking with colleagues in my lab who work with Cython matches what Jed says about Cython; at least one of them writes hand-tuned C to replace the Cython for performance-intensive tasks, and then I believe performance of Python calling the resulting C code is closer to Aron's claim. Also, if you're going to mention PETSc and Python, you might as well mention petsc4py. Last I saw (this was a few years ago), there were no good MPI interfaces for Java. Has that changed? Jun 10 '12 at 21:39
• @GeoffOxberry: The Java MPI bindings exist but have not been updated in nearly a decade. I consider their status dubious. Fortran has numerous I/O options that can be made to go very quickly. I'd recommend exploring Parallel HDF5 (and HDF5, generally). If I/O is truly dominant (more than 50% of run time), more serious measures might be in order, but otherwise, the quality and portability of and HDF-like interface is probably worth it. Jun 10 '12 at 23:04
• @BillBarth: I'll have to check that out. My comment about Fortran I/O comes from the viewpoint of someone once recommending that I write an input file parser in Fortran. It's possible, by enforcing a great deal of structure, but I just haven't seen easily and widely used regex parser or XML parser libraries in Fortran (to give some examples). There's good reason for that: we're the only people using Fortran anymore. Perhaps we're thinking of different use cases. Jun 10 '12 at 23:12

In addition to Aron's very comprehensive answer, I'd take a look at the various threads on scicomp.stackexchange that dealt with the question which programming language to take -- both regarding the speed of programs as well as the question of how easy or hard it is to write and maintain software in these languages.

That said, in addition to what has been written there, let me make a few observations:

(i) You include co-array Fortran in your list. To my knowledge, the number of compilers that actually support it is very small -- and my, in fact, be zero. The most widely available Fortran compiler is GNU gfortran, and while the current development sources parse a subset of co-arrays, I believe that it doesn't actually support any of it (i.e., it accepts the syntax but implements none of the semantics). This is of course a general observation about newer Fortran standards: that the lag with which compilers actually support new standards is measured in several years -- compilers have only fully implemented Fortran 2003 in the last couple of years, and only partially support Fortran 2008. This shouldn't stop you from using any of it if you have a compiler that happens to support what you use, but you must known that you put yourself on a portability island.

(ii) The same is certainly true with C++/Cilk+: Yes, Intel is developing this on a branch of GCC but it's not available in any of the GCC releases and will, likely, not be for a while. You can expect it to take another 2-3 years at least till you will find Cilk+ with the GCC versions installed on typical linux machines.

(iii) C++/TBB is a different story: The TBB has been around for a while, has a very stable interface and is compilable with most any C++ compiler that has existed for the past several years (on linux as well as on windows). We have been using it in deal.II for several years already with good results. There is also a very good book on it.

(iv) I have my very own opinion on OpenMP, namely that it's a solution in search of a problem. It works well for parallelizing the inner loops which is what might be of interest if you have very regular data structures. But it's rarely what you want to do if you need to parallelize something -- because what you really want to do is to parallelize the outer loops. And for that, solutions such as the TBB are much better solutions because they use the mechanisms of the programming language rather than trying to describe what happens outside the language (via #pragmas) and in such a way that you have no access to thread handles, result status indicators, etc, from within your program.

(v) If you're experimental, you might also take a look at the new programming languages that are designed for parallel programming and, in particular, for tasks like the ones you describe. There are essentially two I'd take a look at: X10 and Chapel. I've seen nice tutorials on Chapel, and it seems well designed, though both of course today are insular solutions as well.

• For the record, Intel claims to have a parallelized (distributed memory) co-array Fortran built into their current compilers. We're looking into it at TACC, but I have nothing to report yet. Cray also has an implementation in their compiler, but this is only available on a small integer number of machines throughout the world. I don't think anyone implements the full Fortran 2008 standard with respect to co-arrays, yet, but there's more than nascent support in a few compilers. Cilk+ is, of course, also available with the Intel compilers, but being reliant is probably not yet wise. Jun 10 '12 at 23:10
• The Fortran 2008 standard wasnt approved until late 2010 so it will be a few years before CAF will be widely available. G95 actually had a (non-free) implementation but is not developed anymore (the developer had joined PathScale). Jun 10 '12 at 23:50
• Most of g95 has ultimately ended up in gfortran but it may be that CAF isn't part of that. Jun 11 '12 at 15:57
• I believe the Intel compiler provides good support of co-array. They have built it using mpiexec. It will not be my first choice. The nice thing is that the same implementation can run on shared and distributed memory (I ran a few tests). With opteron shared-memory processors reaching 60-cores at really reasonable prices, I want to see my shared-memory options first. Jun 11 '12 at 22:39

Generally, if you are really serious about this software project, I would suggest a complete re-write in whatever language you yourself feel most comfortable with. It sounds like you will be doing the work alone, and therefore you will get the best results in the language you feel most at home with.

More specifically, though, regarding parallelism, I would encourage you to try to think a bit outside of the box. OpenMP has its strengths, but is stuck in a mindset of taking a sequential code and slapping-on parallelism here and there. The same goes, in essence, for Intels TBB.

Cilk is definitely a step in the right direction, i.e. it forces you to re-think your problem/solution in an inherently parallel setup. What I don't like about it, though, is that it is yet another language. Also, since it can only roughly infer relations between parallel tasks, the scheduler can be quite conservative and may not scale well for certain problems.

The good news is, though, that, again, if you're serious about your implementation, you can do what Cilk does, e.g. re-write your problem as a set of inter-dependent tasks and distribute them over a number of processors/cores, all on your own either using pthreads or misusing OpenMP to spawn processes. A nice example of how this can be done is the QUARK scheduler used in the PLASMA library. A nice comparison of its performance vs. Cilk is given here.

• I will look at the links suggested. The comparison paper is very nice! Thanks! I've been thinking about pthreads but I want the program to be cross-platform. From what I know pthreads have problems under windows (wrong?). Jun 11 '12 at 22:34
• @p3tris: The "p" in pthreads is for POSIX, so it's about as portable as it can be. There are some compliant Windows implementations such as pthreads-win32 or within the cygwin project. Jun 12 '12 at 7:40
• Based on stackoverflow.com/q/2797690/801468 I see there are a lot of stuff needed to sort out to use it. Given that I'm not a programmer, I'd prefer to stick with something more tested. Jun 12 '12 at 9:40

There's been little discussion of coarray fortran in the above comments. At this time, and to my limited knowledge, coarray support in compilers are roughly as follows:

• Cray has a compiler which supports at least the basic coarray features. I've used it to write code that was intended to be "educational", but I'd say that you could write real code in coarray fortran. The syntax and concepts are mostly much simpler than MPI, but as always, there are lotsa traps, and the traps are different from MPI.
• Intel fortran has coarray support built on top of their MPI library. Supposedly this limits their theoretical peak performance, but I haven't seen any metrics.
• Gfortran supports coarrays, but only for a single image (or single rank, in MPI speak). Hence, no real parallelization is available until gfortran 4.8 or 4.9 is out.

Generally, I'd be careful if starting a coarray based code. The syntax is simple and much more convenient than Fortran/C/C++ with MPI, but then, it's just not as full-featured. For instance, MPI supports a lot of reduction operations etc. which could be very convenient for you. It would really depend on your need for a lot of communication. If you want an example, let me know and I can provide you with a few, if I can dig up the files.

Have a look at Spark it's a distributed framework for computations in memory which takes advantage of functional programming. Structure of a program in Spark is very different when compared to MPI, basically you write a code like for single computer, which is automatically distributed as functions to data located in memory. It supports Scala, Java and Python.

Logistic Regression (scala):

//load data to distributed memory
val points = spark.textFile(...).map(parsePoint).cache()
var w = Vector.random(D) // current separating plane
for (i <- 1 to ITERATIONS) {