7

Using an Eigen matrix type where the number of rows and columns is encoded into the type at compile time gives you an edge over LAPACK, where the matrix size is known only at runtime. This extra information allows the compiler to do full or partial loop unrolling, eliminating lots of branch instructions. If you're looking at using an existing library rather ...


7

A little playing with the sequence of numbers generated by the C code shows that the sequence is $z_{i+1}=5z_{i}+273 \mod 2^{16}$ This is a linear congruential generator (LCG). It's easy to show that this LCG has full period (See theorem 7.1 in Law's Simulation Modeling and Analysis, 5th ed. and check the three conditions.) I can't find the generator ...


5

OK, you have a very nice problem, I tried to run some benchmarks. First, I don't have your parameters so I used your small example. Second, since you do not specify the language, I used C + GSL (since I'm not that familiar with C++) #if __STDC_VERSION__ >= 199901L #define _XOPEN_SOURCE 700 #else #define _XOPEN_SOURCE 600 #endif /* __STDC_VERSION__ */ #...


4

Another idea could be to use a generative approach (a program writing a program). Author a (meta)program that spits out the sequence of C/C++ instructions to perform unpivoted** LU on a 10x10 system.. basically taking the k/i/j loop nest and flattening it into O(1000) or so lines of scalar arithmetic. Then feed that generated program into whichever ...


3

Once you're at the level of a long list of expressions, there is little you can still do other than hope that the compiler finds opportunities for optimization at the assembly level. This may bring you 10 or 20% but not something that's going to make a meaningful difference. What you ought to do is go back to the original formula you're trying to implement....


2

Your question leads to two different considerations. First, you need to pick the right algorithm. Hence, the question if the matrices have any structure, should be considered. E.g., when the matrices are symmetric, a Cholesky decomposition is more efficient than LU. When you only need a limited amount of accuracy an iterative method can be faster. Second, ...


1

Would it be possible to do the 1 billion runs method-parallel? Even if there is no way to parallelize the interior of the method body, depending on your problem, you might reach some form of concurrency by evoking the method in parallel. If you do not re-use your output y[100] as an input x[10] to the function then nobody stops you from executing it within ...


1

I would try blockwise inversion. https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion Eigen uses an optimized routine to calculate the inverse of a 4x4 matrix, which is probably the best you're going to get. Try using that as much as possible. http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html Top left: 8x8. Top right: 8x2. ...


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