46

Julia, at this point (May 2019, Julia v1.1 with v1.2 about to come out) is quite mature for scientific computing. The v1.0 release signified an end to yearly code breakage. With that, a lot of scientific computing libraries have had the time to simply grow without disruption. A broad overview of Julia packages can be found at pkg.julialang.org. For core ...


35

I'm not aware of any recent overview articles, but I am actively involved in the development of the PFASST algorithm so can share some thoughts. There are three broad classes of time-parallel techniques that I am aware of: across the method — independent stages of RK or extrapolation integrators can be evaluated in parallel; see also the RIDC (revisionist ...


32

The central issue is the length of the critical path $C$ relative to the total amount of computation $T$. If $C$ is proportional to $T$, then parallelism offers at best a constant speed-up. If $C$ is asymptotically smaller than $T$, there is room for more parallelism as the problem size increases. For algorithms in which $T$ is polynomial in the input size $...


28

If not, is it possible to give a rough order-of-magnitude estimate for how long I should wait before considering it again? My rough, order-of-magnitude estimate of how long it takes computational science languages to mature is around a decade. Example 1: SciPy started in 2001 or so. In 2009, Scipy 0.7.0 was released, and the ODE integrator had an interface ...


27

There is a wide variety of algorithms; Barnes Hut is a popular $\mathcal{O}(N \log N)$ method, and the Fast Multipole Method is a much more sophisticated $\mathcal{O}(N)$ alternative. Both methods make use of a tree data structure where nodes essentially only interact with their nearest neighbors at each level of the tree; you can think of splitting the ...


26

First of all I wish to thanks Aron Ahmadia for pointing me to this thread. As for OpenCL in scientific code: OpenCL is meant to be a low-level API, thus it is crucial to wrap this functionality in some way in order to reach a reasonable productivity. Moreover, as soon as several compute kernels are involved, code can get VERY dirty if OpenCL kernel and ...


21

I believe Julia is worth learning. I have used it to produce a few research finite element codes, and produce them very quickly. I have been over all very pleased with my experience. Julia has enabled a workflow for me that I have found difficult to achieve with other languages. You may use it as a prototyping language like MATLAB, but unlike MATLAB when ...


20

Joblib does what you want. The basic usage pattern is: from joblib import Parallel, delayed def myfun(arg): do_stuff return result results = Parallel(n_jobs=-1, verbose=verbosity_level, backend="threading")( map(delayed(myfun), arg_instances)) where arg_instances is list of values for which myfun is computed in parallel. The main ...


19

Your matrix is of size 15,000 x 15,000, so you have 225M elements in the matrix. This makes for roughly 2GB of memory. This is much more than the cache size of your processor, so it has to be loaded completely from main memory in every matrix multiplication, making for approximately 100GB of data transfers, plus what you need for the source and destination ...


17

There are basically two major, commercial choices out there: DDT from Allinea (which is what we use at TACC) and Totalview (as mentioned in the other comment). They have comparable features, are both actively developed, and are direct competitors. Eclipse has their Parallel Tools Platform, which should include MPI and OpenMP programming support and a ...


17

The first thing is to recognize that you can do this using BLAS. If you data matrix is $X = [x_1 x_2 x_3 ...] \in \mathbb{R}^{m\times n}$ (each $x$ is a column vector corresponding to one measurement; rows are trials), then you can write the covariance as: $$ C_{ij} = E[x_i,x_j] - E[x_i] E[x_j] = \frac{1}{n} \sum_k x_{ik} x_{jk} - \frac{1}{n^2} \left(\sum_k ...


16

I must give the curmudgeon answer. My productivity has never been improved by any of the suggestions above. They are slow and expensive compared to my preferred option in parallel: one gdb session per process. Each gdb can connect to an MPI process and sit in an xterm (this happens automatically in PETSc using -start_in_debugger). I have used this for 15 ...


16

We are currently writing a paper that contains a number of comparable plots, and we more or less had the same problem. The paper is about comparing the scaling of different algorithms over the number of cores, which ranges between 1 and up to 100k on a BlueGene. The reason for using loglog-plots in this situation is the number of orders of magnitude involved....


16

As Paul states, without more information, it is hard to give advice without assumptions. With 10-20 variables and expensive function evaluations, the tendency is to recommend derivative-free optimization algorithms. I am going to disagree strongly with Paul's advice: you generally need a machine-precision gradient unless you're using some sort of special ...


15

A reduction implemented using MPI_Allreduce() is reproducible as long as you use the same number of processors, provided the implementation observed the following note appearing in Section 5.9.1 of the MPI-2.2 standard. Advice to implementors. It is strongly recommended that MPI_REDUCE be implemented so that the same result be obtained whenever the ...


15

Although this post is now two years old, in case someone stumbles across it, let me give a brief update: Martin Gander recently wrote a nice review article, that gives a historical perspective on the field and discusses many different PINT methods: http://www.unige.ch/~gander/Preprints/50YearsTimeParallel.pdf There is now also a community website which ...


14

Parallel geometric multigrid is straightforward to implement on structured grids. Algebraic and unstructured multigrid are more technical, see this answer for links to implementations. In a multiplicative method (e.g. $V$-cycles), only one level can be computed on at a time. Since the number of levels is $\log_{c} N$ where $N$ is the number of degrees of ...


14

One library to consider is BoxLib. Its key features (from the website) are: Support for block-structured AMR with optional subcycling in time Support for cell-centered, face-centered and node-centered data Support for hyperbolic, parabolic and elliptic solves on hierarchical grid structure C++ and Fortran90 versions Supports hybrid programming model with ...


14

Georg Hager wrote about this in Fooling the Masses - Stunt 3: The log scale is your friend. While it is true that log-log plots of strong scaling are not very discerning on the high end, they allow for showing scaling across many more orders of magnitude. To see why this is useful, consider a 3D problem with regular refinement. On a linear scale, you can ...


13

Look at the fast multipole method. It is highly scalable and $O(n)$. It allows trading off between precision and cost. Here's an example where it is run at 42 Tflops on a GPU cluster.


13

To give a theoretical aspect to this, $NC$ is defined as the complexity class that is solvable in $O(log^c n)$ time on a system with $O(n^k)$ parallel processors. It is still unknown whether $P=NC$ (although most people suspect it's not) where $P$ is the set of problems solvable in polynomial time. The "hardest" problems to parallelize are known as $P$-...


13

To the best of my knowledge, Numpy does not support independent streams. Indeed, getting independent streams from the Mersenne Twister (Pythons RNG) is notoriously difficult although it can be done. Consider using the RandomGen package. It is fully compatible with Numpy, and provides you with the PCG64 generator, supporting up to $2^{63}$ independent ...


12

HPC is a blend of mathematics, computation, computer science, and applications. You need to be able to understand all of them to be truly successful in the long run. However, you do not necessarily need to achieve the same level of proficiency in all of them. In the computation versus mathematics argument, for an engineer, I'd argue the numerical ...


12

Good is a relative term, and it will depend on the nature of the problem, the nature of the algorithm, and properties of the hardware involved. The only absolute reference point is ideal scaling (100% efficiency). You can claim your scaling is good if it is better than what anyone else has achieved for the same problem, or if it's "close" to ideal for ...


11

I'm a heavy user of OpenFOAM, so naturally I'd recommend it. It has a large amount of features including combustion models (though not necessarily precisely what you need) and has been used together with Canterra by other people. If you need a solver for a specific equation that hasn't been implemented yet, you can pretty much literally write your equations. ...


11

I am extremely surprised that there is no mention of conditioning or the shape of the spectrum in your discussion, as it will be the decisive property in whether or not iterative methods can beat dense methods. As an extreme example, suppose that your dense matrix is some small perturbation of the identity matrix. Then most iterative methods will converge ...


11

Boost Graph Library and LEMON As Daniel mentions in his comprehensive answer, the most full-featured general C++ library is the Boost Graph Library. There is a new distributed-memory extension capable of doing some basic algorithms such as breadth-first and depth-first search, minimum spanning trees, and connected components search, but I am not very ...


11

I would say that there are a number of reasons why there are no computational science contests besides the potentially massive computational resources required. Time limits: Writing scientific computing code is usually not something that you want to rush. A lot of emphasis is on making sure it is correct, and thorough consideration of test/corner cases. ...


11

Both the standard cluster and custom supercomputer (Anton) versions of molecular dynamics at D. E. Shaw Research are both deterministic and parallel invariant. That is, a test run on a single core generates the same bits as a massively parallel run. The techniques include Integer summation: Although each force term is computed in floating point, the total ...


10

I like aeismail's answer, but I'm going to provide an alternative perspective. In optimization, it's impossible to really learn the field without understanding real analysis. Even before you tackle numerical issues, you need to understand notions of convergence of sequences, because you are going to prove in classes that algorithms converge. You're going to ...


Only top voted, non community-wiki answers of a minimum length are eligible