# Tag Info

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

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

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

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

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

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

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

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

9

One of the authors of fenics, A. Logg, have written a very good paper on datastructures of storing meshes. The paper is A. Logg (2009). Efficient Representation of Computational Meshes http://arxiv.org/abs/1205.3081 In fact it's always a tradeoff between storing all the topological informations (nodes around nodes, faces around nodes, etc...) OR having to ...

9

The standard way of approaching this is really to attempt a Cholesky factorization and check to see if such a factorization exists. This is both fast and realiable. Here it is in MATLAB notation: A = zeros(3); % some matrix [~,p] = chol(A) If the input matrix is not positive definite, then $p$ will be a positive integer, e.g. $p=1$ and MATLAB will ...

9

You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq. JuliaDiffEq and DifferentialEquations.jl has a large feature set dedicated to efficiently integrating computationally-difficult differential equations. It has a simple ...

8

I am assuming that you are setting the error tolerance at 1e-12. You are correct that when an adaptive scheme accepts the current step size, it assumes the 5th order scheme was, for all intents and purposes, the "correct" answer. However this is only when it accepts the current step. If the difference between the 4th and 5th order steps are too large, it ...

8

CPLEX is commercial grade and solves very large LPs and IPs. If CPLEX didn't pan out for you, then switching to a different solver may not be the answer. Here's one SO Question on good solvers. Instead, here are three suggestions for you, rather than going for the solver speed or the specific algorithm. Focus on your formulation. How many constraints and ...

8

In deal.II, we basically only use vectors. Maps are too slow and scatter data all around memory, so we typically don't use them if the keys are integers and within a given range. For example, for the connectivity between cells, you can do arrays (STL vectors) in which you store neighbor indices and so that neighbor indices $4i\ldots 4i+3$ correspond to cell \$...

8

@BillGreene points to the "return value optimization" as a way around the fundamental problem, but this actually only helps for one half of it. Assume you have code of this form: struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); }; ExpensiveObject operator+ (ExpensiveObject &obj1, ExpensiveObject &obj2) { ...

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