Hot answers tagged

48

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

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


21

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


18

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

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


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


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


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

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

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

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


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

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

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

No, use one of the free available linear algebra libraries. A discussion about different libraries can be found here: Recommendations for a usable, fast C++ matrix library?


10

You might want to look at Python with NumPy/SciPy. (There are many more packages available; for example, pandas is a dedicated data analysis library). The syntax is somewhere between Matlab and C and so should not be too hard to pick up. There's a helpful NumPy for Matlab users page, too. One of the advantages is that it is fairly easy to "outsource" parts ...


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

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


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.


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