Hot answers tagged

46

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


32

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


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


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


20

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


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


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

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

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

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

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

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

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

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

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


13

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


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

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


11

OPT++ is used internally by Dakota (Sandia), which is much more than an optimization library and is released under the LGPL. You should also take a look at TAO (ANL), released under a BSD-like license. An introduction to both OPT++ and TAO can be found here. Other alternatives are MOOCHO, NOMAD, and HOPSPACK, which are, as far as I know, also LGPL licensed. ...


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

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


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


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