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


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


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


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


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


22

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


20

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

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


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

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

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

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


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

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


10

Particle and domain decomposition are directly connected to the two main methods of speeding up force calculations for systems with limited-range interactions - Verlet neighbour lists and cell linked lists. If you'd like to get into details, there is a pretty nice book from Allen and Tildesley, called Computer Simulation of Liquids, considered by many to be ...


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 trick is to interleave each process's LCG stream: for $p$ processes, we modify the LCG $$ x_{n+1} := a x_n + c\;\;\; (\bmod m),$$ to be $$ x_{n+p} := A_p x_n + C_p\;\;\; (\bmod m),$$ where $A_p$ and $C_p$ effectively step forward $p$ steps. We can quickly derive them by expanding the original LCG step: $$ x_{n+2} = a (a x_n + c) + c\;\;\; (\bmod m)...


9

You should also look at libMesh. It's targeted at finite element methods, but other than that, I think it checks most of your boxes. Unlike BoxLib, it's a fully unstructured, mixed element type library, which is to stay that it supports tets, pyramids, prisms, and hexahedra in the same mesh. It also has one of the largest sets of integration rules for high-...


9

Start by grocking Amdahl's Law. Basically anything with a large number of serial steps will benefit insignificantly from parallelism. A few examples include parsing, regex, and most high-ratio compression. Aside from that, the key issue is often a bottleneck in memory bandwidth. In particular with most GPU's your theoretical flops vastly outstrip the ...


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


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