Hot answers tagged

17 votes

Looking for Runge-Kutta 8th order in C/C++

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. ...
user avatar
12 votes
Accepted

How to document math formulations in scientific computing codes?

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 ...
user avatar
  • 8,452
11 votes
Accepted

Memory Management - Why certain initializations order are faster?

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 (...
user avatar
  • 8,452
10 votes
Accepted

Numerical evaluation of the Exponential Integral Ei by rational Chebyshev approximations fails

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)}{\...
user avatar
  • 11.4k
9 votes

Use of Java or Scala in HPC?

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 ...
user avatar
9 votes
Accepted

Looking for Runge-Kutta 8th order in C/C++

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 ...
user avatar
  • 2,473
8 votes
Accepted

C - OpenMP, MPI, Serial Program

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 ...
user avatar
  • 11.4k
8 votes
Accepted

What PRNG function is this?

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 ...
user avatar
8 votes
Accepted

C standard for computational science

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 ...
user avatar
  • 2,199
8 votes
Accepted

Suming up the series $1+1/2+...+1/n$ in C

I see several problems with your code: type error in scanf (n is an integer, not a float) i=0 and not i==0 in the for loop ...
user avatar
  • 196
7 votes

Numerical integration giving incorrect sign

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 ...
user avatar
6 votes
Accepted

Performance of adding eight numbers sequentially vs. in a tree

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 ...
user avatar
  • 11.4k
6 votes

Use of Java or Scala in HPC?

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 ...
user avatar
  • 10.8k
6 votes

Use of Java or Scala in HPC?

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 ...
user avatar
  • 2,473
6 votes

Do BLAS routines compute their respective operations with minimum error?

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 ...
user avatar
  • 11.4k
6 votes
Accepted

Using GSL for basic operations

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 ...
user avatar
  • 6,091
6 votes

I've developed a derivative-free optimization method, looking for comments

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 ...
user avatar
  • 2,116
5 votes

Looking for Runge-Kutta 8th order in C/C++

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 ...
user avatar
5 votes

BLAS, LAPACK or ATLAS for Matrix Multiplication in C

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 ...
user avatar
  • 5,744
5 votes
Accepted

How to deal with big numbers in intermediate calculations?

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 ...
user avatar
  • 11.4k
5 votes

C standard for computational science

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 ...
user avatar
  • 551
5 votes
Accepted

Recommendations for ODE solvers for stiff equations

This is a huge open ended question, but I'll copy the recommendation section from the current release of DifferentialEquations.jl: Stiff Problems For stiff problems at high tolerances (>1e-2?) it ...
user avatar
5 votes

float in C: Size of the exponent

$1.2\times 10^{-38}$ is a decimal value. You need to convert it to binary exponential notation to see how the bits are set. Another small mistake: $2^6+...2^1+2^0$ does not equal $2^{27}$, it equals $...
user avatar
4 votes
Accepted

What is wrong with this matrix multiplication?

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-...
user avatar
  • 11.4k
4 votes

solver linear system equation

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 ...
user avatar
4 votes

Kummer's confluent hypergeometric for complex arguments in C/C++?

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: ...
user avatar
4 votes
Accepted

What do C, C++ and Java have that Fortran 2003 don't?

One thing that C++ includes that Fortran doesn't have is extensive support for generic programming and compile-time code ...
user avatar
4 votes

Jacobi iteration for finite difference: when to stop?

For solving systems you shouldn't be comparing the results between each iteration but rather computing the residual. If you consider the matrix representation to be in the standard form: $ A\cdot y ...
user avatar
4 votes
Accepted

C or fortran library to solve linear 2D/3D elliptic PDE

Most of the widely used finite element libraries are written in C++. If all you really care for -- and if all you will ever care for -- is solving an elliptic PDE on a rectangle, then it's probably ...
user avatar
4 votes
Accepted

Number of decimal of float and double in C

Prof. Bangerth's comment is completely correct. To add more detail, we can refer to IEEE754 standard which defines the floats as ...
user avatar

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