# Tag Info

56

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

31

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

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

23

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

23

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

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

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

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

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

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

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

12

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

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

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

10

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

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

8

For [CU]BLAS, there is a wrapper called 'thunking' in the CUDA toolkit (src/fortran_thunking.{c,h}) that takes pointers from CPU memory and does all the GPU allocation/copying for you. You can plug it into your code with a preprocessor statements like #define ZGEMV CUBLAS_ZGEMV #define ZGEMM CUBLAS_ZGEMM ... For LAPACK, Magma has CPU-side interfaces for ...

8

From my experience Julia isn't ready for (scientific) everyday use yet (I'm talking about the stabilized version 0.4 of march 2016). The language itself is nice, feature rich and consistent; a refreshing step forward from both matlab and python. There are certainly uses cases where Julia is a good choice even at this early stage. But if you need reliable and ...

8

From my many years writing FEM software, I believe that the statement that DG schemes are better suited to parallelization than CG schemes is apocryphal. It is frequently used in introductions of DG papers as a justification for DG methods, but I have never seen it substantiated with a reference that actually investigated the question. It is similar to every ...

8

In fact, the precise total number of operations is very rarely used as a measure of computational cost. Instead, you will most often see the computational order (i.e. $\mathcal{O}(n^3)$). This "big O" notation roughly means that the number of operations is proportional to $n^3$ and tells you how the total number of operations scales as the number of unknowns ...

8

I think some of your issues are more important than others and some of your emphasis is misplaced. In pursuing overhead, you are in danger of making your program unmaintainable. It is easier to write a common program and direct surplus effort somewhere more interesting. I apologize for pontificating like this. If statements. From a strict programming ...

8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...

8

Like Brian said, the Xeon Phi cores are not at all comparable to the CUDA ones. The problem with the Phi is that it's somewhere between two horses. If you are doing highly parallel floating point calculations, NVIDIA will provide you with something like 3 times the performance at 1/4th of the price. For double precision the gap is smaller, but NVIDIA still ...

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