Hot answers tagged

19

Possibly one could start with the function $\mathtt{expm1}$ which is part of the C99 standard, and calculates $e^x-1$ accurately near $x=0$.


17

This is an instance of cancellation error. The C standard library (as of C99) includes a function called expm1 that avoids this problem. If you use expm1(x) / x instead of (exp(x) - 1.0) / x, you won't experience this issue (see graph below). The details and solution of this particular problem are discussed at length in Section 1.14.1 of Accuracy and ...


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

Since you're already using C++ and your matrices are symmetric positive definite, I would perform an unpivoted $LDL^T$ factorization of $Q$ and also of $12I-Q-J$. Here I'm assuming that $12I-Q-J$ is also positive definite, otherwise the $LDL^T$ will require pivoting for numerical stability (it's also possible that even though it's not positive definite, ...


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

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

In alphabetical order (disclaimer: I am the main author of Elemental): DPLASMA Distributed Parallel Linear Algebra Software for Multicore Architectures (DPLASMA) is a relatively recent and ongoing effort by Bosilca et al. to extend PLASMA to distributed-memory machines. Version 1.0.0 supports distributed Cholesky factorizations, among many other operations....


7

Have you had a look at VMD? I used it ages ago to produce movies from simulation snapshots. Way back then, it could read a sequence of PDB files, render them (or generate POV-Ray scripts to raytrace them), and store them as individual images. I then used mencoder to generate MPEG-4 files out of the stills. Those were the days. I haven't used VMD since, but ...


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


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

What have you tried so far? This completely naive implementation manages to compute 7 (maybe 7.5) digits in 2.5 seconds on my laptop: #include <iostream> #include <complex> #include <cmath> #include <iomanip> int main () { const double alpha = 1; std::cout.precision(16); std::complex<double> sum = 0; for (unsigned int ...


6

Bill answered the first part, so I'll only answer the second question. An MPI send is blocking if it does not return until it is safe to modify the send buffer and a receive is blocking if it does not return until the receive buffer contains the newly-received message. In practice, outside of buffered sends (thanks, Hristo Iliev), this implies that ...


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

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

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

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

MPI_Probe allows you to test for a message without actually receiving it. You must complete all non-blocking communications with an appropriate communications completion fuction like MPI_Wait and friends, otherwise the runtime will not free up internal resources associated with the communications leading to resource leaks and other problems. For example, you ...


5

Without some information about the construction of these $12\times 12$ positive definite real symmetric matrices, the suggestions to be made are of necessity fairly limited. I downloaded the Armadillo package from Sourceforge and took a look at the documentation. Try to improve performance of separately computing $\det(Q)$ and $\det(12I - Q - J)$, where $J$...


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

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

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


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


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