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

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

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


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


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

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

I was able to get an answer on the math stackexchange site. You could find the answer here! It supposedly has to do with Hermite Normal Forms! Thank you, Heterotic for the answer.


3

This formula seems to assume that one operation counted by $\omega(n)$ takes as long as one I/O operation counted by $h(n,p)$. So to apply it, you have to be able to count operations in your code and scale them to the same time scale. The further elaboration dividing the numerator and denominator by $\omega$ probably scales this out effectively, but you ...


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


3

My first question is, when it comes to writing the software what is the best practice for the separation of these steps? Should a program do all of these steps or is it better to divide the steps between different scripts and save the data in between? How you separate the steps depends on many things, like how long it takes to solve the ODEs, what software ...


3

From Numerical Linear Algebra by Trefethen and Bau, page 175, it seems that the formula is \begin{align} \sum_{k=1}^{n}\sum_{j = k + 1}^{n} (2(n - j + 1) + 1). \end{align} Eyeballing it, it seems to agree with the formula given by Boyd in his convex optimization notes: $(1/3)n^{3} + 2n^{2}$.


3

Common subexpression elimination and inlining code are like the bread and butter of compiler optimizations, the first and easiest optimizations that a compiler can do. There is every reason to think the performance of the code will be exactly the same (while losing the comments), and there is no particular reason to think that a compiler would have trouble ...


3

The Fortran standard doesn't specify what precisions and ranges of floating point numbers need be supported by a given compiler. All it says is that at least two different kinds of real numbers need be supported and that larger of them takes exactly twice the memory to represent a number as the other (smaller) kind. What is actually supported depends upon ...


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