49

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 ...


30

There is one real and one practical reason. First, MPI was developed at a time when machines had exactly one processor core and when we wanted to couple different machines. It is today used on clusters of tens of thousands of machines, each of which happens to have many cores but the point is that it's still separate machines. Now, a processor core on ...


27

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 ...


24

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 ...


22

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 ...


18

The simple example from electromagnetics (EM) would be performing a parallel frequency sweep for a frequency-domain simulation, say, full-wave extraction of network parameters (S, Y, Z, etc) for a device. Since the simulation for each frequency point is highly independent from another, the simulation can be embarrassingly parallelized across different cores, ...


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

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 ...


15

What you're looking for is Numba, which can auto parallelize a for loop. From their documentation from numba import jit, prange @jit def parallel_sum(A): sum = 0.0 for i in prange(A.shape[0]): sum += A[i] return sum


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

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

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

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 ...


11

Without assuming something special on my_function choosing multiprocessing.Pool().map() is a good guess for parallelizing such simple loops. joblib, dask, mpi computations or numba like proposed in other answers looks not bringing any advantage for such use cases and add useless dependencies (to sum up they are overkill). Using threading as proposed in ...


11

Defining the auxiliary variable $y=Bx$ yields the following algebraically equivalent expanded system, $$\underbrace{\begin{bmatrix} 0 & A \\ B & -I \end{bmatrix}}_{K} \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_{u} = \underbrace{\begin{bmatrix} b \\ 0 \end{bmatrix}}_{f},$$ which you could solve with GMRES or another nonsymmetric Krylov method. ...


11

High-quality video encoding is something like this. The search space is so huge that it requires branching to prune it rapidly, but GPUs are terrible at that. Modern CPU short-vector SIMD works well for this, working on contiguous chunks of 16 to 64 bytes of data. And while still being tightly coupled to the CPU core which can branch efficiently on SIMD ...


10

The Thomas algorithm is very efficient because its operation count is very low and because data accesses are very likely to be cache hits once data is initially read from memory. There are two loops. The first loop traverses the data forward. Each element of the lower, main and upper triangle, along with the right-hand-side vector (which is typically ...


9

The MUMPS sparse direct solver can handle symmetric indefinite systems and is freely available (http://graal.ens-lyon.fr/MUMPS/). Ian Duff was one of the authors of both MUMPS and MA57 so the algorithms have many similarities. MUMPS was designed for distributed-memory parallel computers but it also works well on single-processor machines. If you link it ...


9

I'm a happy user of GoogleTest with a C++ MPI code in a CMake/CTest build environment: CMake automatically installs/links googletest from svn! adding tests is a one-liner! writing the tests is easy! (and google mock is very powerful!) CTest can pass command-line parameters to your tests, and exports data to CDash! This is how it works. A batch of unit-...


9

256 equations is a relatively small number. All of the usual integrators, such as those included in Matlab, Maple or Mathematica should have no real problem with equations of this size and should be able to return answers in a fraction of the time it would take an algorithm you would implement yourself, because they use sophisticated explicit/implicit and ...


9

One reason DG methods may receive more attention as a parallel method is that it is readily seen that the method is inherently local to an element. The coupling in DG methods is weak, as it only occurs through adjacent edges (or faces in 3d). So, for triangles or quads DG will communicate to three or four processers at most, respectively. Whereas CG methods ...


9

Like you say, using the Mersenne Twister for parallel computations is almost always done incorrectly, as the correct method is tricky to implement. By far the easiest and best answer would be to move away from the Mersenne Twister entirely, and use something like the PCG family, which provides multiple streams out of the box. The Mersenne Twister is known ...


9

Wolfgang Bangerth's answer is totally correct, and I only want to add one practical aspect. Portability across hardware Let's say you write a research code from scratch. You have a powerful multi-core shared memory machine at your department which can do the job. If you start out with a thread-based implementation, you can reach good performance on this ...


9

GPUs work with the model SIMD (single instruction multiple data) i.e. they execute an instruction over multiple data. To give an idea: under CUDA technology when you have got an if-then-else condition the two branches are executed in sequence over the respective data. In your question, the condition to favor a CPU suggests a MISD or MIMD model, i.e. ...


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