66

First of all, as skillman and Dan have pointed out, profiling is essential. I personally use Intel's VTune Amplifier on Linux as it gives me a very fine-grained overview of where time was spent doing what. If you're not going to change the algorithm (i.e. if there will be no major changes that will turn all your optimizations obsolete), then I'd suggest ...


50

Language designers face many choices. Ken Kennedy emphasized two: (1) better abstractions and (2) higher- or lower-level (less or more machine-like) code. While functional languages like Haskell and Scheme focus on the former, traditional scientific-computing languages like Fortran and C/C++ focused on the latter. Saying that one language is faster than ...


40

I'm going to break up my answer into three parts. Profiling, speeding up the python code via c, and speeding up python via python. It is my view that Python has some of the best tools for looking at what your code's performance is then drilling down to the actual bottle necks. Speeding up code without profiling is about like trying to kill a deer with an ...


38

Scientific software is not that much different from other software, as far as how to know what needs tuning. The method I use is random pausing. Here are some of the speedups it has found for me: If a large fraction of time is spent in functions like log and exp, I can see what the arguments to those functions are, as a function of the points they are ...


37

In Matlab, the ‘\’ command invokes an algorithm which depends upon the structure of the matrix A and includes checks (small overhead) on properties of A. If A is sparse and banded, employ a banded solver. If A is an upper or lower triangular matrix, employ a backward substitution algorithm. If A is symmetric and has real positive diagonal elements, ...


37

If you want to see what a\b does for your particular matrix you can set spparms('spumoni',1) and figure exactly what algorithm you were impressed by. For example: spparms('spumoni',1); A = delsq(numgrid('B',256)); b = rand(size(A,2),1); mldivide(A,b); % another way to write A\b will output sp\: bandwidth = 254+1+254. sp\: is A diagonal? no. sp\: is band ...


31

I think the (first order) right thing to do is look at the ratio of flops to bytes needed in the algorithm, which I call $\beta$. Let $F_{\mathrm{max}}$ be the maximum flop rate of the processor, and $B_{\mathrm{max}}$ the maximum bandwidth. If $\frac{F_{\mathrm{max}}}{\beta} > B_{\mathrm{max}}$, then the algorithm will be bandwidth limited. If $B_{\...


23

The design of Fortran allows the compiler to perform stronger optimizations in some cases, optimizations that are not generally available to C. One famous example is the handling of aliasing. In Fortran, you can access a specific memory area only though the specific symbol associated with that memory area. This knowledge allows the compiler to employ smart ...


23

Not only do you have to have intimate knowledge of your compiler, you also have intimate knowledge of your target architecture and operating system. What can affect performance? If you want to squeeze every last ounce of performance, then every time you change your target architecture, you will have to tweak and re-optimise your code. Something which was ...


22

I don't see why one has to be the "winner"; this isn't a zero-sum game, where flop counts and memory access have to drown the other out. You can teach both of them, and I think they both have their uses. After all, it's hard to say that your $O(N^4)$ algorithm with $O(N)$ memory accesses is necessarily going to be faster than your $O(N \log N)$ algorithm ...


19

The closest positive answers to your question that I could find is for sparse diagonal perturbations (see below). With that said, I do not know of any algorithms for the general case, though there is a generalization of the technique you mentioned for scalar shifts from SPD matrices to all square matrices: Given any square matrix $A$, there exists a Schur ...


18

MKL (from Intel) is optimized for Intel processors, and probably has the "upper hand" there in many cases. But it is also "famous" for choosing the "worst" code-paths for AMD processors, as described here.


17

My problem with expression templates is that they are a very leaky abstraction. You spend a lot of work writing very complicated code to do a simple task with nicer syntax. But if you want to change the algorithm, you have to mess with the dirty code and if you slip up with types or syntax, you get completely unintelligible error messages. If your ...


17

Matlab interprets sequences of multiplications and/or divisions from left to right. Hence $A*B*C*v$ is much more expensive than $A*(B*(C*v))$, as you have two matrix products and one matrix-vecor product in place of three matrix-vector products. On the other hand, $A*(B*(C*v))$ should be slightly faster than if you save the intermediates in separate ...


16

For sparse matrices, Matlab uses UMFPACK for the "\" operation, which, in your example, basically runs through the values of a, inverts them, and multiplies them with the values of b. For this example, though, you should use b./diag(a), which is a lot faster. For dense systems, the backslash-operator is a bit more complicated. A brief description of what is ...


15

Let me give an example based on experience. Most libraries I use from a day to day basis use OOP in some way. OOP is able to hide the complexity required for many domains, it is not a mechanism that really helps with performance. What can happen is that a library is able to use specific optimizations based upon the object hierarchy, but for the most part ...


14

I think by and large, template metaprogramming has been found to be unusable in practice -- it compiles too slow, and the error messages we get are just impossible to decipher. The barrier to entry for newcomers is also just too high when using metaprogramming. Of course, generic programming is an entirely different issue, as witnessed by Trilinos, deal.II (...


14

For starters, I wouldn't use intermediate variables, but brackets. Unless, of course, you're interested in the intermediate results, but I'm guessing not. I tried the following in Matlab: >> N = 500; >> A = rand(N); B = rand(N); C = rand(N); v = rand(N,1); >> tic, for k=1:100, A*B*C*v; end; ...


14

Since the matrices are so small, all of the cost is going to be in call overhead. If you will do the transformation many times, it will be faster to precompute D=A*B*C once and then for each vector apply v_f=D*v_i. You could also consider bringing this out to a mex file.


14

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


13

BLAS is not monlithic. BLAS1 and BLAS2 are memory bandwidth limited, and there is not much you can do to speed them up beyond the obvious (loop unrolling, cache blocking for level 2). BLAS3 is more interesting and the prototypical benchmark here is matrix-matrix multiplication. To my knowledge GOTOBlas has always been the clear winner here, see for example ...


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 don't think Fortran is that close to the metal (see other answer) but it tends to optimize very easily. Loops are simple, and the language readily supports vectorization extensions (okay when I used it in my first job we were targeting a wide range of vector big iron). There is also the large factor of inertia. A lot of numeric code is in Fortran, so ...


12

If $(D_{i} + A)$ is diagonally dominant for each $i$, then recent work by Koutis, Miller, and Peng (see Koutis' website for work on symmetric diagonally dominant matrices) could be used to solve each system in $\mathcal{O}(n^2 \log(n))$ time (actually $\mathcal{O}(m\log(n))$ time, where $m$ is the maximum number of nonzero entries in $(D_{i} + A)$ over all $...


12

Others have commented on the issue of how difficult it is to write ET programs as well as the complexity of understanding error messages. Let me comment on the issue of compilers: It is true that a while back one of the big issues was finding a compiler that's compliant enough with the C++ standard to make everything work and make it work portably. As a ...


12

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


12

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


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