39

As an extension to moyner's answer, the on-chip sqrt is usually an rsqrt, i.e. a reciprocal square root that computes $a \rightarrow 1/\sqrt{a}$. So if in your code you're only going to use $1/r$ (if you're doing molecular dynamics, you are), you can compute r = rsqrt(r2) directly and save yourself the division. The reason why rsqrt is computed instead of ...


20

$n \log\log n$ is between $n$ and $n \log n$, and is a relatively common one to find in the wild.


18

Short answer, you want to have the leftmost index on the innermost loop. In your example, the loop indices would go k, j, i and the array indices would be i, j, k. This has to do with how MATLAB stores the different dimensions in memory. For more, see #13 of this reddit post.


16

The square root is implemented in hardware on most processors, that is, there are specific assembly instructions and the performance should be comparable in most languages because it is very hard to muck up the implementation. You will probably never be able to beat the FSQRT instruction, since it was designed by some smart hardware designer. How it is ...


15

My rule of thumb is that if you can compute some quantity efficiently (good utilization of the FPU) in less than 50 flops per double precision value, it's better to recompute than to store. The trend, which has been steady for decades, is for floating point capability to improve faster than memory performance, and is not likely to relent due to physical ...


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


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


11

A somewhat longer answer that explains why it's more efficient to have the left most index varying most rapidly. There are two key things that you need to understand. First, MATLAB (and Fortran, but not C and most other programming languages) stores arrays in memory in "column major order." e.g. if A is a 2 by 3 by 10 matrix, then the entries will be ...


10

Particle and domain decomposition are directly connected to the two main methods of speeding up force calculations for systems with limited-range interactions - Verlet neighbour lists and cell linked lists. If you'd like to get into details, there is a pretty nice book from Allen and Tildesley, called Computer Simulation of Liquids, considered by many to be ...


9

The integral in question is also known as the Boys function, after the British chemist Samuel Francis Boys who introduced its use in the early 1950s. A few years ago, I needed to compute this function in double precision, as fast as possible but accurately. I managed to achieve a relative error on the order of $10^{-15}$ across the entire input domain. It ...


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


7

I've had to deal with similar problems before, and my favourite solution is to use Memory-mapped I/O, albeit in C... The principle behind it is quite simple: instead of opening a file and reading from it, you load it directly to the memory and access it as if it were a huge array. The trick that makes it efficient is that the operating system doesn't ...


7

With your choice of fuzzy language in your comment, you are carefully skirting the answer to your question. The number of nonzeros in a matrix is a decent proxy for cost to apply it to a vector, since the cost is dominated by the time to load its entries from cache. Roughly speaking, if $$\mathrm{nnz}(AB) > \mathrm{nnz}(A) + \mathrm{nnz}(B),$$ then you ...


7

On top of $O(n\log(\log(n)))$, there's also $O(n \log^*(n))$ in which $\log^*$ is the number of times the logarithm function must be applied in order for the result to be less than or equal to 1. For instance, if you already know an Euclidean minimum spanning tree, the Delaunay triangulation may be discovered in $O(n\log^*(n))$ time. More extremely, one ...


7

Generally there isn't a way to compute the inverse of a sum of Kronecker products. However, suppose there is a factor in common, let's say $I_T$ here and your sum is $$ A = K \otimes I_T + I_T \otimes \Sigma $$ Then solving linear systems with $A$ becomes equivalent to solving a Sylvester equation. Thus perhaps one way to approach your problem would be ...


6

There are so many variables that play into the speed with which a program runs that it is impossible to tell from just your description. For example: Did other programs run at the same time? What is the clock speed of your processor? How was the program compiled and which processor was it optimized for? How long does the program run to begin with -- for ...


6

As Jed already remarked, it all depends on the matrices. But you can get an idea by considering, for example, the 5-point stencil of the Laplace equation in 2d. There, each row of the stiffness or mass matrix $A$ has five entries. But the product of two such matrices has 13 entries (apply the five point stencil to every location of the five point stencil). ...


6

I'm assuming your question comes from the observation that the I/O causes a significant overhead in your whole analysis. In that case, you can try to overlap I/O with computation. A successful approach depends on how you access the data, and the computation you perform on that data. If you can identify a pattern, or the access to different regions of the ...


6

My limited understanding of AD parallels what Matt has said. To speed up the computation of derivatives, the structure of the expression graph must exploit sparsity and scarcity in the set of Jacobian matrices. (See this paper by Griewank for more insight.) The software engineering tricks would likely be in the AD code itself to restructure the expression ...


6

You need to declare a LDLT object like this: LDLT<MatrixXd> ldlt(A); x = ldlt.solve(b); y = ldlt.solve(x); ...


6

$L_0$ leads to a combinatorial problem and is theoretically intractable, and practically slow. To a large degree, intractability of $L_0$ is often the reason we use $L_1$ instead. $L_1$ and $L_{\infty}$ leads to linear programming which is tractable, and you can exploit structure to develop even faster algorithms $L_2$ leads to second-order cone ...


6

There are many different ways to do this. One of the standard is a work-precision plot where you plot the amount of time or function calls that it takes in order to achieve a certain level of accuracy. You can find tons of examples at DiffEqBenchmarks.jl. Generally you slide a timestep or adaptivity tolerances along a window and plot all of the (time,error) ...


6

This can be done using a technique called the $\eta$ (eta) factorization. This method is commonly used in implementations of the simplex method for linear programming and can be found in many textbooks on linear programming. The procedure is as follows: Find an LU factorization of $A_{1}$, $PA_{1}=LU$ Using this factorization, you can easily solve $...


5

There is no magic way to speed up an arbitrary operation by using more memory. Similarly, there is no way to compress arbitrary data. It so happens that a lot of data has exploitable redundancy, but even then, changing the algorithm (e.g. changing a spatial discretization, applying formal model reduction methods, or using an adjoint model) tends to offer ...


5

You might take a look at Numerical Methods for Special Functions by Amparo Gil, Javier Segura, and Nico M. Temme.


5

First of all, note that this forum is for Computational Science, not Computer Science. These are different fields, with Computational Science being scientific computing, more like computational mathematics and scientific simulations. That said, even though the example in question is not relevant to this forum, there are things that should be discussed here. ...


4

I'd take a look at Abramowicz & Stegun's book, or the newer revision that NIST has published a couple of years ago and that's available online I believe. They also discuss ways to implement things in a stable way.


4

deal.II uses variants of your option (2). This is of course also essentially the way you store sparse matrices and is very efficient since you only have two arrays (the array that for each cell stores the beginning index of the data that corresponds to this cell, and the array that contiguously stores all the data of all cells). You traverse both arrays ...


4

By "Domain decomposition is a better choice only when linear system size considerably exceeds the range of interaction, which is seldom the case in molecular dynamics" the authors of that (very old) GROMACS paper mean that if the spatial size of the neighbour list is of the order of 1 nm, and the simulation cell is only several nanometers, then the overhead ...


4

Whether saving information makes something faster depends on how much time it takes to retrieve from storage and how much time it takes to recalculate the variables on the fly. As a general rule, it takes something on the order of fractions of nanoseconds to get something stored in cache or to perform a CPU instruction. Going to RAM can take anywhere from ~...


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