15

For the first part of my question, I found this very useful comparison for performance of different linear interpolation methods using python libraries: http://nbviewer.ipython.org/github/pierre-haessig/stodynprog/blob/master/stodynprog/linear_interp_benchmark.ipynb Below is list of methods collected so far. Standart interpolation, structured grid: http:/...


14

To make more robust comparisons (on linux), you can : 1) On Intel CPUs the turbo overclocks your CPU. This is controlled by the temperature of the CPU, so it can behave differently from one run to the other. On Linux, you can block the frequency of the CPU as follows. For example, for 2.4GHz: echo 1 > /sys/module/processor/parameters/ignore_ppc for ...


14

Here's the deal with GPUs. On a GPU, every single core is slow. Really slow. However, you have thousands of cores. If you can effectively use the thousands of cores at a time, then your algorithm will run better on the GPU. If you cannot, then it will run much slower on the GPU. Linear algebra is one domain where parallelism is really well established. ...


14

There has been considerable recent interest in numerical linear algebra using mixed precision with some combination of 16, 32, 64, and 128 bit floating point arithmetic. For example, a low precision factorization of a matrix can be used used to precondition a higher precision iterative solution. Since the factorization takes $O(N^{3})$ operations and the ...


13

In general, both methods of performance comparisons have their place. Comparing cpu time is in a sense the most interesting metric, because at the end of the day you are really interested in which of the methods is faster. (But make sure that the termination criteria are comparable; e.g., that both methods yield an approximation with the same accuracy). ...


13

There are some differences, however they aren't necessarily in hardware or specs. Note that this is all information I have gained from forums or news releases, so take it all with a grain of salt. The first is the "scalability and reliability" (source). The K20 was designed to sit in a cluster system and run at full tilt 24/7. The Titan is more designed ...


13

This list is nowhere near complete, but hopefully the size of it will give a hint as to the scale of possible factors. I am assuming you are compiling the code from source on your platform of choice. Software Standard Library Performance Lin. Alg. Library Performance (if the software links to outside libraries) Compiler Choice Compiler Optimization ...


12

I think it depends on your purpose. If you are trying to assess the overall performance of the code or environment, then I'd encourage you to run it however you think most people will run it on a desktop environment: leave things open but make sure nothing is crunching in the background or hogging all the memory. The biggest culprits, in my experience, are ...


12

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: return y*sin(pi<Real>()*x)/pi<Real>(); The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and ...


11

Fortran's allocatable variables are automatically deallocated when the variable goes out of scope (see http://www.fortran90.org/src/best-practices.html#allocatable-arrays). This means that it is not possible to create a memory leak by failing to deallocate an allocatable array. This is one of the big benefits of using allocatable arrays rather than pointers. ...


11

MATLAB's \ (aka mldivide) command does not blindly compute the inverse of the matrix. Instead, it uses one of several algorithms based on the type of matrix (see the "Algorithms" section of http://www.mathworks.com/help/matlab/ref/mldivide.html). In the case of a triangular matrix, MATLAB will use a triangular solver which is at least as good as yours in ...


11

In the 1980's era Intel 80x86 architecture, there was a scalar floating point unit that had instructions like FSIN, FCOS, etc. for computing functions like sin and cos. These functions were implemented in microcode and might take 100's of CPU cycles to execute. Later, Intel added Streaming SIMD Extensions (SSE) which gave the processor parallel floating ...


10

Doing one-off best rational approximations can often be accomplished by "manual" iterations of the Remez algorithm: interpolate a rational approximation with (relative or absolute) alternating sign errors at an initial guess for interpolation points, locate one (or more) points where the actual error exceeds that of the interpolation points and pivot (...


10

I think in practice the impact is limited and will be limited. The reason why it is limited right now is that the big finite element packages are careful to write code that is portable, and so they do not yet use C++11 language constructs in their own codes. Of course, they will benefit from code like the one you show where, even without having to change ...


9

This is usually caused by trying to use a threaded MKL combined with MPI, resulting in over-subscription. Either explicitly configure PETSc to use non-threaded MKL or add MKL_NUM_THREADS=1 to your environment.


9

OpenBlas is quite fast, so you can link it to LAPACK. Have you tried precompiled version of LAPACK/BLAS from your CPU vendor? For example AMD ACML (free) or Intel MKL (free on linux for non-commercial and non-academic use)? You simply need to unpack and run install file. In my opinion the only advantage of using ATLAS is then when you use some unusual CPU. ...


9

tl;dr Use loops My numbers indicate that ifort is smart enough to recognize the loop, forall, and do concurrent identically and achieves what I'd expect to be about 'peak' in each of those cases. gfortran, on the other hand, does a bad job (10x or more slower) with forall and do concurrent, especially as N gets large. Both ifort and gfortran seem to ...


9

Chances are that the evaluation of the functions is the most time consuming part of this computation. If that's the case, then you should focus on improving the speed of func() rather than trying to speed up the integration routine itself. Depending on the properties of func(), it's also likely that you could get a more precise evaluation of the integral ...


9

Consolidating the comments: No, you are very unlikely to beat a typical BLAS library such as Intel's MKL, AMD's Math Core Library, or OpenBLAS.1 These not only use vectorization, but also (at least for the major functions) use kernels that are hand-written in architecture-specific assembly language in order to optimally exploit available vector extensions (...


8

Be careful of urban myths here. If two compilers generate the same assembly code, then the resulting programs will have the same performance. For any given piece of logic, there is a program in assembly language that minimizes its execution time. That program doesn't care which compiler generated it. That said, compiled languages exist to make life easier ...


8

If this actually works, and it seems to, that'd be awesome because we can get a lot of undocumented data about the cache of a GPU. A frustrating aspect of high performance computing research is digging through all the undocumented instruction sets and architecture features when trying to tune code. In HPC, the proof is in the benchmark, whether it is High ...


8

The benchmarks on the Julia website 1 2 include R and Matlab as competitors. Note that these are benchmarks focusing on testing the pure speed of the language, not the quality of the underlying linear algebra or FFT libraries. The speed for operations that are outsourced to these libraries (such as a large matrix multiplication) can vary a lot depending on ...


8

I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS array, otherwise the interface will create an internal copy with the correct storage order. Python 3.6.2 (default, Jul 22 2017, 21:19:22) [GCC 7.1.1 20170516] on linux Type "help", "copyright", "...


7

It is important to preallocate correctly. This is almost certainly the reason why your assembly was slow. If you are starting with a dense matrix representation, just scan through it once counting the number of nonzeros per row, then call MatSeqAIJSetPreallocation(). See this FAQ. The option MAT_IGNORE_ZERO_ENTRIES is really intended to be used when there is ...


7

Prior answers to this question have covered most of the salient points, but I want to add one comment with respect to this: does MKL have the upper hand for some tasks? The MKL team is in a unique position to know about future Intel instruction sets and their implementations in specific processors. Furthermore, they have access to proprietary processor ...


7

The repository package is not safe to use with threading due to the way it was compiled. I reported the bug on the Lapack forum, but it will take a long time for workarounds or solutions to trickle down into the repository. If you compile it yourself, be sure to add the "-frecursive" to gfortran.


7

Possible? Yes. Useful? No. The optimizations I'm going to list here are unlikely to make more than a tiny fraction of a percent difference in the runtime. A good compiler may already do these for you. Anyway, looking at your inner loop: for (s=0,j=a;j<b;j+=h){ func2 = func(j+h); s = s + 0.5*(func1+func2)*h; func1 = func2; ...


7

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


7

You can convert $\mathbf{b}-\mathbf{x}$ into polar coordinates, and do the dot product in this system. This changes $((\mathbf{b}-\mathbf{x})\cdot\nabla)^n\frac{1}{r}$ to $$\left((\mathbf{b}-\mathbf{x})_r\frac{\partial}{\partial r} + (\mathbf{b}-\mathbf{x})_{\theta}\frac{1}{r}\frac{\partial}{\partial \theta}\right)^n\frac{1}{r}$$ Here, I am using the ...


7

The appropriate and fastest library depends on several things. Which Bessel functions (only J, Y & Hankel or modified Bessel functions I & K too), for which types of arguments (real or complex, integer, fractional or general order)? Amos's libraries are written in Fortran-77 (there are Fortran-90 coverted versions of TOMS 644 on a mirror of Alan ...


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