Hot answers tagged

35

First, if your undergraduates are like ours and had no prior introduction to computers, expect to spend some time teaching them how to use basic stuff like using a proper editor (i.e., not MS Word), the command line, etc. I think the answer somewhat depends on where you set the focus of your course (or what you are required to teach). For example: How ...


23

In 2014, I would've said Python. In 2017, I wholeheartedly believe that the language to teach undergraduates is Julia. Teaching is always about a tradeoff. On one hand, you want to choose something that is simple enough that it is easy to grasp. But secondly, you want to teach something that has staying power, i.e. something that can grow with you. The ...


19

Let me first answer why I think C++ interfaces to MPI have generally not been overly successful, having thought about the issue for a good long time when trying to decide whether we should just use the standard C bindings of MPI or building on something at higher level: When you look at real-world MPI codes (say, PETSc, or in my case deal.II), one finds ...


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


18

There seems to be quite a bit of confusion about how to apply multi-step (e.g. Runge-Kutta) methods to 2nd or higher order ODEs or systems of ODEs. The process is very simple once you understand it, but perhaps not obvious without a good explanation. The following method is the one I find simplest. In your case, the differential equation you would like to ...


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


15

A difficulty with any of these types of questions is that the answer is highly community-dependent. To answer some of your questions in haphazard order: MATLAB is used a lot both in academia and in industry. One of the reasons it's used quite a bit in industry is because it is taught in academia. I know for a fact that MATLAB is used at Lincoln Laboratory ...


15

Feature test macros: HPC is generally stuck on old compilers or compilers with partially conformant implementations. This can help ease the pain of working on the custom architectures common in HPC. Example: #ifdef __cpp_lib_source_location #include <source_location> #endif ... #ifdef auto sl = std::source_location(); std::cerr << "Error at line ...


14

In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=(). I think you're right that it's creating temporaries, which is too slow, but I think the reason for why it's doing that is wrong. Armadillo, like any good C++ linear ...


14

If you have constants that will not change before runs, declare them in a header file: //constants.hpp #ifndef PROJECT_NAME_constants_hpp #define PROJECT_NAME_constants_hpp namespace constants { constexpr double G = 6.67408e-11; constexpr double M_EARTH = 5.972e24; constexpr double GM_EARTH = G*M_EARTH; } #endif //main.cpp using namespace ...


13

I will address only the comparison of C to C++. While it is true that anything written in C can be ported to C++ with a few syntactic touch-ups, the communities have different values. The C library community, more than almost any other, values binary stability. Binary stability is critical for low-level libraries to avoid inflicting constant pain on those ...


13

Another alternative that may be in line with your train of thought is to use a namespace (or nested namespaces) to properly group constants. An example might be: namespace constants { namespace earth { constexpr double G = 6.67408e-11; constexpr double Mass_Earth = 5.972e24; constexpr double GM = G*Mass_Earth; }// constant properties ...


12

Two examples of libraries that use modern C++ constructs: Both the eigen and armadillo libraries (linear algebra) use several modern C++ constructs. For instance, they use both expression templates to simplify arithmetic expressions and can sometimes eliminate some temporaries: http://eigen.tuxfamily.org http://arma.sourceforge.net/ http://hpac.rwth-...


12

Consolidating the comments: No, you are very unlikely to beat a typical BLAS library such as Intel's MKL, AMD's Math Core Library, or OpenBLAS.1 These not only use vectorization, but also (at least for the major functions) use kernels that are hand-written in architecture-specific assembly language in order to optimally exploit available vector extensions (...


12

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: return y*sin(pi<Real>()*x)/pi<Real>(); The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and ...


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

As Misha and Geoff Oxberry pointed out, Mathematica really has a different focus (just because you can pound in a nail with a screwdriver doesn't mean you should). So I take your question as being "If I know Matlab, why should I learn Python?" [Edit: and so, apparently, did you.] For all intents and purposes, Matlab is the English of scientific computing -- ...


11

It is true that compilers are getting better and better at auto-vectorization, and for basic coefficient-wise operations like 2*A-4*B a library like Eigen cannot do much better than recent compilers. However, for slightly more complicated expressions like matrix products, reductions, transposition, powers, etc. the compiler cannot do much. On the other hand, ...


10

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'. Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix: most of the subroutines have a TRANS argument to accommodate for operation on the transpose matrix or on a matrix stored with a different memory layout. (...


10

I think in practice the impact is limited and will be limited. The reason why it is limited right now is that the big finite element packages are careful to write code that is portable, and so they do not yet use C++11 language constructs in their own codes. Of course, they will benefit from code like the one you show where, even without having to change ...


10

The Thomas algorithm is very efficient because its operation count is very low and because data accesses are very likely to be cache hits once data is initially read from memory. There are two loops. The first loop traverses the data forward. Each element of the lower, main and upper triangle, along with the right-hand-side vector (which is typically ...


10

I would suggest taking a look at Deal.II. It uses the STL, it's own iterators, shared pointers, etc. The various linear solvers can use the various matrices because of how it was designed. I haven't come across any use of move semantics, but that doesn't mean they aren't there. Here is a link.


10

I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you have more-or-less one-to-one correspondence between the bilinear form and a (sparse) matrix. The user can then use the resulting matrix as he or she sees fit. ...


9

If you can stand to use complex arithmetic, simultaneous iteration methods might be preferable for computing all the roots of your polynomial. The simplest simultaneous iteration method, the (Weierstrass-)Durand-Kerner method, is effectively equivalent to applying Newton-Raphson to the Vieta relations relating the coefficients and roots of a polynomial, ...


9

There are a couple subtleties to your question that I think are important: You're comparing an interpreted language (Python) to a compiled language (C++). Most scientific and engineering software is developed with a heavy Linux (and UNIX) bias, and is not usually known for cross-platform compatibility or great user support (big libraries, of course, ...


9

The Levenberg-Marquardt method can be used to minimize any problem of the form: $ \min f(x)=\sum_{i=1}^{m} f_{i}(x)^{2} $ However, if the objective functino to be minimized is not a sum of squares, then the method is no longer applicable.


9

Chances are that the evaluation of the functions is the most time consuming part of this computation. If that's the case, then you should focus on improving the speed of func() rather than trying to speed up the integration routine itself. Depending on the properties of func(), it's also likely that you could get a more precise evaluation of the integral ...


9

FTensor is a lightweight, header only, fully templated library that includes ergonomic summation notation. It has been tested extensively in 2, 3, and 4 dimensions, but should work fine for any number of dimensions.


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


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