32

Although LAPACK has some incredibly optimized code, it can still be worth it to write your own version in a few cases. The most important reason (and the reason they make you do it in your course). It's a great learning experience. You will never fully understand something like the QR algorithm (which is actually incredibly complicated to get the details ...


19

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


15

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

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

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

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

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


12

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


12

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


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

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


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


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

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


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

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

There is a nice discussion on StackOverflow regarding floating point vs integer operations. In short, the performance of the operations depends a lot on processor architecture how the data is stored in memory and in which order it is accessed if (and which) SSE/AVX/AVX2/etc instructions are used (and how efficiently) This probably provides some insight ...


8

@ThiysSteel covers a lot, here is another perspective I find important: Even if you have available excellent implementations of any algorithm you might need, you still need to understand some of the internals in order to make good use of your library. Some examples: Any linear-algebra library can compute an inverse matrix really fast and as accuratly as ...


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

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

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


7

In general, I agree with Chris's comment that using a compiled language with the allocation of the matrices on the stack can help significantly. Several possibilities if we are limited to Python and numpy: consider np.array vs np.matrix, it might happen that np.matrix is faster than np.array matrix-matrix product (it is unclear what you are using now, and ...


7

You can calculate GFLOP rates this way, but the numbers are pretty meaningless on today's hardware: Floating point operations require a variable number of clock cycles. An addition is generally cheaper than a multiplication, but each generally takes more than one clock cycle of the 2.8 billion cycles you quite. When you have hyperthreading, you have two ...


7

To clarify @ThiysSteel's good answer: The point is not to attempt to surpass the optimization of code written by very experienced people, who'd wrangled with it for decades. The point is to acquaint yourself with the issue that were addressed in that code, to appreciate the "niceties", and not be naive about related matters in the future. (It's ...


7

I used to try to optimize code via using assembly language (as opposed to C). I had some clear success, where a real-time microphone array worked with assembly language, but it would not work at all in the C language. However, since that time, more than 20 years have passed. It's very difficult to achieve this sort of improvement nowadays. Compilers have ...


6

Traditionally, I'd say that people more or less understood the number of floating-point operations required to solve a problem ($O(n^3)$ and the like), and so reporting performance as a rate had meaning. Readers could then obtain the run time via the method you describe, but they could also compare that to the theoretical peak performance of the hardware ...


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