Hot answers tagged

62

As so often, the choice depends on (1) the problem you are trying to solve, (2) the skills you have, and (3) the people you work with (unless it's a solo project). I'll leave (3) aside for the moment because it depends on everyone's individual situation. Problem dependence: Fortran excels at array processing. If your problem can be described in terms of ...


44

It's a bad idea because vector needs to allocate as many objects in space as there are rows in your matrix. Allocation is expensive, but primarily it is a bad idea because the data of your matrix now exists in a number of arrays scattered around memory, rather than all in one place where the processor cache can easily access it. It's also a wasteful storage ...


37

I think that both C++ and Fortran are good enough and work well. However I think that Fortran is better for numeric scientific computing, for algorithms that can be expressed using arrays and don't need other sophisticated data structures, so in fields like finite differences/elements, PDE solvers, electronic structure calculations. Fortran is a domain ...


32

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


31

I'm also throwing my two cents in kind of late, but I've only just seen this thread and I feel that, for posterity, there are a few points that desperately need to be made. Note in the following that I will talk about C and not C++. Why? Well, otherwise it's apples and oranges to compare a full-fledged dynamically typed object-oriented language with ...


22

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 try and break down your requirements: Maintainability Reading/writing text data Strong interfaces/capability for LU factorizations Sparse linear solvers Performance and scalability to large data From this list, I would consider the following languages: C, C++, Fortran, Python, MATLAB, Java Julia is a promising new language, but the community is ...


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

In addition to the reasons Wolfgang mentioned, if you use a vector<vector<double> >, you'll have to dereference it twice every time you want to retrieve an element, which is more computationally costly than a single dereferencing operation. One typical approach is to allocate a single array (a vector<double> or a double *) instead. I've ...


17

My problem with expression templates is that they are a very leaky abstraction. You spend a lot of work writing very complicated code to do a simple task with nicer syntax. But if you want to change the algorithm, you have to mess with the dirty code and if you slip up with types or syntax, you get completely unintelligible error messages. If your ...


17

If you're looking for a good bound on your rounding error, you don't necessarily need an aribtrary-precision library. You can use running error analysis instead. I wasn't able to find a good online reference, but it's all described in Section 3.3 of Nick Higham's book "Accuracy and Stability of Numerical Algorithms". The idea is quite simple: Re-factor ...


17

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


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

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


16

From my 15 years of thinking about scientific software: If your code runs 25% faster because you write it in Fortran, but it takes you 4 times as long to write it (no STL, difficulty implementing complex data structures, etc), then Fortran only wins if you spend a significant fraction of your day twiddling thumbs and waiting for your computations to finish. ...


15

You should probably start with the LAPACK implementation, ?gtsv, e.g., dgtsv. If you want a distributed-memory version, then you might want to start with ScaLAPACK's p?gtsv. EDIT: Since your matrix does not change very often, you can avoid redundantly factoring the tridiagonal matrix by breaking up the LAPACK routine ?gtsv into the factorization step, ?...


15

Let me give an example based on experience. Most libraries I use from a day to day basis use OOP in some way. OOP is able to hide the complexity required for many domains, it is not a mechanism that really helps with performance. What can happen is that a library is able to use specific optimizations based upon the object hierarchy, but for the most part ...


15

I believe this discussion has come up a number of times on the PETSc list. My main reasons are: The C++ standard states that std::complex is only defined for the float, double, and long double datatypes. Thus it cannot be used for other datatypes, such as quad-precision. The standard makes no guarantees about the stability of the complex arithmetic. The ...


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

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


14

I think by and large, template metaprogramming has been found to be unusable in practice -- it compiles too slow, and the error messages we get are just impossible to decipher. The barrier to entry for newcomers is also just too high when using metaprogramming. Of course, generic programming is an entirely different issue, as witnessed by Trilinos, deal.II (...


13

My approach has been to use C++ for everything but computational kernels, which are usually best written in assembly; this buys you all of the performance of the traditional HPC approach but allows you to simplify the interface, e.g., by overloading computational kernels like SGEMM/DGEMM/CGEMM/ZGEMM into a single routine, say Gemm. Clearly the abstraction ...


13

I have always thought that we should use it in our own project, deal.II, because it is higher level than pure MPI and can save a few lines of code here and there. That said, what I learned over the years is that most high-level code doesn't actually have that much MPI code to begin with -- the 600,000 lines of code in deal.II have only ~50 calls to MPI. That'...


13

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

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


12

Others have commented on the issue of how difficult it is to write ET programs as well as the complexity of understanding error messages. Let me comment on the issue of compilers: It is true that a while back one of the big issues was finding a compiler that's compliant enough with the C++ standard to make everything work and make it work portably. As a ...


12

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


12

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

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


11

If you want to write something general-purpose, you can do it either with shell scripts if it is something very simple, as Pedro suggests, or aggregate in a higher-level mathematical programming language such as Python or MATLAB. I agree that plain text files are useful for smaller amounts of data, but you should probably switch to binary data for anything ...


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