Hot answers tagged

17

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods,...


12

I prefer using doxygen that supports C++ and LaTeX comments, both inline and as separate equations. This way, you will keep your comments, including, say, the rigorous mathematical formulation of the algorithm, very close to the source code. The generation of the documentation can be included in the overall workflow (say, a Makefile or CMake target, ...


11

Such an effect happens because of how the data of the int** a is stored in memory (as per C/C++). This question on StackOverflow has answers with some more details (in particular, a difference between int** and int[][] that many users noted in the comments), and how it looks like an array of arrays - it's just laid out contiguously in memory. It's worth ...


10

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...


9

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods. Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL. GSL is published under the GPL license, and Boost Odeint under the ...


9

Java has been around for almost 20 years now as a major programming language, but it hasn't caught on in scientific computing so far. I think that's a good indicator for what's going to happen in the future. My take is that the issue isn't speed. Most people are probably willing to give up 20% of performance (or even a factor of 2) if they would be vastly ...


8

I think some of your issues are more important than others and some of your emphasis is misplaced. In pursuing overhead, you are in danger of making your program unmaintainable. It is easier to write a common program and direct surplus effort somewhere more interesting. I apologize for pontificating like this. If statements. From a strict programming ...


8

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


8

In theory, as the original authors, you're free to pick and name a standard, then expect others to follow it. In practise, if you're supporting an HPC system, then your choice is likely to be restricted among the compiler standards that the given system you're using for testing has for its tool stack (you are testing your software, aren't you?). This might ...


7

This might be an accuracy problem in computing the second term, because of those large exponentials when $x \gg 1$. I would first work on that term: gather $e^x$ out from numerator and denominator and simplify it out (leaving a summand $e^{-x}$ in the numerator). However, rather than guessing it might be necessary to understand better where the accuracy ...


6

I would recommend looking at a time-dependent example because PETSc can provide a lot more diagnostics if you formulate at that level. For example, you could use a Rosenbrock method, an additive Runge-Kutta IMEX method, or others. Some involve Newton iteration, but that is not required and is not always the best approach. To use those methods, you can ...


6

You're probably better served writing a C wrapper for the Fortran implementation you linked to: Colavecchia, F. D., Gasaneo, G., "f1: a code to compute Appell's F1 hypergeometric function", Computer Physics Communications, Volume 157, Issue 1, p. 32-38 (2004), found at http://cpc.cs.qub.ac.uk/summaries/ADSJ. The R package appell wraps that implementation of ...


6

Whether your code is efficient or not, it will not work for any numbers over 32767 as written. This is because the int data type is a signed type of 16 bit length. One bit is used for the sign and 15 are used for the value of the integer, making the largest storable integer 215-1=32767. If you wish to support numbers larger than this, you will need to use a ...


6

I think your analysis is basically right. Some notes. 1. Pipelining is the wrong word here; what you're looking at here is data dependency. A CPU pipeline splits an individual instruction into multiple steps, and different steps of consecutive instructions can then be executed concurrently. A data dependency, on the other hand, is a situation where an ...


6

It seems unlikely to me. The Java MPI APIs haven't been worked on in years (so you're wrong about #4), and the JVM's floating-point performance is notoriously poor. Java may out perform C/C++ or Fortran in some areas due to rapid thread creation and easy memory management, but these aren't the bottlenecks in typical scientific programs. As to your #5, the ...


6

I would argue that Java will in fact REDUCE productivity when compared with modern c++, or even with modern Fortran for the purpose of scientific computing. Writing A = B*C+2*D is just so much more readable than A = B.mult(C).add(D.mult(2)) Assuming the code above deals with arrays, both C++ and fortran will also produce significantly more efficient ...


6

BLAS routines do not typically use stable summation algorithms. In the case of gsl, you can look up its source code online - the source of gsl's sdot is contained in gsl/cblas/source_dot_r.h, and contains this loop: for (i = 0; i < N; i++) { r += X[ix] * Y[iy]; ix += incX; iy += incY; } It's just a straightforward sum. The corresponding ...


6

Of course it makes sense to use the GSL (or another library for that matter) for several reasons: Don't reinvent the wheel. The work has been done, you can spend your time on more useful things. If you do decide to implement these basic things yourself, the risk that your code will probably contain some bugs and will be slower, less memory efficient etc ...


6

I would like to hear comments from users that have some practical models (e.g. black-box hyperparameter optimization) which are still needed to be solved acceptably - whether this method works or not for their models, possibly with the description of the model. Looks like you want somebody to invest what may be considerable time and energy in trying out ...


5

summarizing some points: If it's a long-term integration of a non-dissapative model, a symplectic integrator is what you're looking for. Otherwise, since it's an equation of motion, Runge-Kutta Nystrom methods will be more efficient than a transformation to a first order system. There are high order RKN methods due to DP. There are some implementations, ...


5

Some thoughts from someone who has worked a fair amount in compiled languages, and has done a tiny bit of FVM: Typically, if you have experience programming in C, you sketch out a high-level description (pseudocode) of what you would like to do. Then you look for libraries that might implement the data structures and capabilities you need for your high-...


5

Yes, you want to call the BLAS routine DGEMM. The place to start for how to call it from C is to look at the documentation for DGEMM, which you can find online. Then you want to understand how to call FORTRAN routines from C (DGEMM, like all the standard BLAS routines has a FORTRAN calling convention). For example, this document https://computing.llnl.gov/...


5

I think there's a simple way to do this. You have a rational function of identical cosh/sinh terms, where every expression is a homogeneous polynomial in cosh/sinh, and the only problem is that these exponential terms overflow. The function does not diverge as these terms approach infinity, so if you divide every numerator and denominator by the same power ...


5

You should definitely jump to C99, or newer(!). The C99 standard introduced the restrict keyword. Loosely speaking, with this keyword you can inform the compiler that A[i] and B[j] do not access the same memory location. In that case the compiler can generate better optimized code. For example, it makes it easier for the compiler to auto-vectorize code. ...


4

Use $\exp\left(-\tfrac1{x^2}\right) =\frac1{\exp\left(\tfrac1{x^2}\right)}$ and expand in the denominator, for example as $\exp\left(-\tfrac1{x^2}\right) \approx \frac1{1+1/x^2+1/(2x^4)}=\frac{4x^4}{(1+2x^2)^2+1}$ which is correct for $1/x^6$ below machine precision, i.e. for $x>1000$ and has the qualitatively correct behavior for smaller x. You can ...


4

I'd imagine that most users of random number generators are ultimately interested in floating-point values. This is why the Double precision SIMD-oriented Fast Mersenne Twister (dSFMT) exists. However, there is newer C code for the WELL RNG that returns unsigned long values. Looking at the code, it appears that the earlier version was casting unsigned long ...


4

Since Doug has pointed out the need for a bigger integer datatype to do the problem you want to tackle, let's talk about the logic of your program and improvements. Your approach is (A) input alpha as the number to factor, (B) make a list of all prime numbers less than alpha, and (C) see which of these divide alpha and add the ones that do to an array ...


4

Note: I haven't run your code. Perhaps this is a problem with the row-major/column-major conventions at play here: http://docs.nvidia.com/cuda/cublas/index.html#data-layout. You seem to use row-major matrices, but BLAS uses the column-major convention. When you pass row-major $A$ and $B$ to dgemm, it implicitly interprets them as $A^t$ and $B^t$, giving ...


4

You haven't told us whether your matrix is positive definite or not- this is a very important factor in selecting a solver. You also haven't told us whether the matrix is well conditioned or poorly conditioned. A third important question is whether you need a fairly accurate solution (e.g. accurate to 10 digits) or whether a less accurate solution (e.g. ...


4

One possibility is to wrap the arbitrary precision interval implementation in Arb. This will not be as fast as a dedicated double precision implementation, but it might still be fast enough. Note 1: this code requires Arb version 2.8.0 or later. #include "acb_hypgeom.h" void hyp1f1ix(double * re, double * im, double a, double b, double x) { long prec; ...


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