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.


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


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


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


8

The inverse Langevin function $\mathcal{L}^{-1}(x)$ is an odd function. Therefore one needs to consider only approximation on the interval $[0, 1]$; the negative half-plane is treated via symmetry about the origin. Because of the singularity at unity, a polynomial approximation is not a fruitful direction of the investigation, although some authors have ...


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

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

So, you are comparing a generally slower Matlab implementation of algorithm A to a generally faster C++ implementation of algorithm B, and still getting the advantage for A. I would say, congratulations, you certainly have a stronger point now, since the "competing" algorithm is given an advantage. It's worth no note though, that implementation in C++ is ...


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

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


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

$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

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

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

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

There are infinitely many, since $O(n(\log n)^\alpha) \subsetneq O(n(\log n)^\beta)$ for any $\alpha<\beta$. So, in particular, $O(n) = O(n(\log n)^0) \subsetneq O(n(\log n)^\alpha) \subsetneq O(n\log n)$ for any $\alpha\in (0,1)$.


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


4

Your loop 2 does not preallocate $C$. Hence it gets resized $n$ times, which requires many reallocations. Try adding C = zeros(n, n) before the for loops. [EDIT: so does the first one, so probably this isn't the true reason.] (That said, I second ChrisRackauckas's comment --- there is little point in testing timings under Matlab at this level.)


4

Expanding on my previous comment. There are efficient algorithms to solve linear systems of the form $$ (K \otimes M+I_T\otimes \Sigma)\operatorname{vec}(X) = \operatorname{vec}(B), $$ so in practice you can just use those instead of an explicit inverse. Explicit inverses are overrated, in some fields. :) The trick is that these linear systems are ...


4

You might want to use the fact that: $$ ||A||_2=\sigma_\max(A) $$ where $\sigma_\max$ is the largest singular value. If you are interested in details, this Math SO question should be interesting. Thus, $$ ||A^{-1}||_2=\frac{1}{\sigma_\min(A)} $$ where $\sigma_\min$ is the smallest singular value. You certainly want to avoid the actual calculation of the ...


3

The factor you miss is that when you do matrix-matrix multiplications, you have to reload data from memory. Since your matrices are big enough that they won't fit into the cache, you access main memory which is incredibly slow compared to a floating point operation. To give you an idea how slow: a floating point multiplication on a typical processor today ...


3

For dense matrices or order $n$, the theoretical number of FLOPS required for matrix product is $2n^3$. In your case, since you are multiplying two matrices of order 1000, the number of FLOPS should be, in MFLOPS, 2000. Therefore, the 848.3 MFLOPS performance of your code is less than 50% of the theoretical performance for matrix multiplication. Scaling $n$,...


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