Parallel Scientific Computation Software Development Language?

Question

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.

Answer 1

Let me try and break down your requirements:

• Maintainability
• Reading/writing text data
• 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.

Reading/writing text data

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:

• Recommendations for a usable, fast Java matrix library?

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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) {
  val gradient = points.map(p =>
    (1 / (1 + exp(-p.y*(w dot p.x))) - 1) * p.y * p.x
  ).reduce(_ + _)
  w -= gradient
}
println("Final separating plane: " + w)

There's an extension to called MLib (Machine Learning library) which uses a Fortran library for some low level computations (for Python I guess numpy is used). So, the idea is simple, concentrate on your algorithm and leave optimizations to lower levels (order of processing, data distribution, etc.).