32 votes
Accepted

Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

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). ...
user avatar
  • 1,189
20 votes
Accepted

The real myth of GPU (specifically CUDA) really speed up FEM/CFD

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 ...
user avatar
19 votes

What is the preferred and efficient approach for interpolating multidimensional data?

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-...
user avatar
14 votes
Accepted

How to compare runtimes of two algorithms in a reproducible way

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 ...
user avatar
14 votes
Accepted

Increasing computational performance by using 16 bit numbers

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 ...
user avatar
13 votes
Accepted

Beating typical BLAS libraries matrix multiplication performance

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 ...
13 votes
Accepted

Whittaker-Shannon interpolation: Accuracy dies with speedup; can it be fixed?

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: ...
user avatar
  • 1,030
13 votes
Accepted

Is there an algorithm or graph theory that allows me to not need to store an intermediate matrix when calculating AT*Y1*A + BT*Y2*B?

BLAS may not have a function to compute what you are asking for, but the product $$ Y_N = A^TY_AA + B^T Y_B B $$ means that the entries $(Y_N)_{ij}$ are defined by $$ (Y_N)_{ij} = \sum_{k,l} (A^T)...
user avatar
12 votes
Accepted

Fortran 90/95: Deallocating variables

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 ...
user avatar
12 votes
Accepted

Hardware performance, floating point functions

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 ...
user avatar
11 votes

Why is my MATLAB code for back-substitution slower than the backslash operator?

MATLAB's \ (aka mldivide) command does not blindly compute the inverse of the matrix. Instead, it uses one of several algorithms ...
user avatar
10 votes

Open source implementation of rational approximation to a function

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 ...
user avatar
  • 3,189
9 votes
Accepted

Loop optimization with f2py, Cython and Numba

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 ...
user avatar
  • 3,789
8 votes
Accepted

Is R or Matlab currently faster?

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 ...
user avatar
8 votes
Accepted

Integer operations vs floating point operations

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 ...
user avatar
  • 8,382
8 votes
Accepted

Fastest Way to Mutiply $10^4$ 2x2 Matrices

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 ...
user avatar
  • 8,382
8 votes

How to properly calculate CPU and GPU FLOPS performance?

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 ...
user avatar
8 votes

Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

@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 ...
user avatar
  • 181
7 votes
Accepted

What is the fastest opensource implementation of Bessel functions computation?

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,...
user avatar
7 votes
Accepted

How to evaluate a series of derivatives?

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{...
user avatar
7 votes

Performance differences between ATLAS and MKL?

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 ...
user avatar
  • 2,001
7 votes

Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

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 ...
user avatar
7 votes

Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

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 ...
user avatar
6 votes

Lua and Fortran (or Python) I/O

This is probably not the answer you are looking for, but I wanted to state it anyway: Your choice of programming languages introduces two difficulties you will encounter as your program grows. First,...
user avatar
6 votes
Accepted

Performance of adding eight numbers sequentially vs. in a tree

I think your analysis is basically right. Some notes. 1. Pipelining is the wrong word here; what you're looking at here is data dependency. A CPU pipeline splits an individual instruction into ...
user avatar
  • 11.4k
6 votes
Accepted

What is the state of the art algorithm for diagonalizing real symmetric matrices?

The QR/Francis algorithm is the go-to choice for dense eigenproblems, but there are a few competitors around: The Jacobi algorithm (like QR, another algorithm with an unfortunate name, which can be ...
user avatar
6 votes
Accepted

Why is the speed of the parts of the LU-decomposition so different?

First, don't forget to also time the LU decomposition in a loop! Otherwise it's not really a fair comparison. If I do that, I get the following timings: ...
user avatar
6 votes
Accepted

Implementation of Jacobi iteration

If you really want to speed it up, I'd suggest changing your iteration to a SSOR iteration. This can converge much faster than jacobi, and if you are already OpenMP parallelized, there's not much you ...
user avatar
  • 1,902
6 votes
Accepted

Is operation count a reliable predictor of performance when comparing two formulations?

No. There a number of aspects of modern computer architecture that make operation counts unreliable as a way to compare the performance of algorithms. These include memory caches, vector ...
user avatar
6 votes

Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

This might be tangent but i feel it worth adding. A common way for user code to beat libraries is to exploit structure the libraries are not aware of. Block diagonal matrices are easier to work in ...
user avatar

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