Hot answers tagged

24 votes

What language should I use when teaching an undergraduate course in computer programming?

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 ...
user avatar
18 votes
Accepted

Applying the Runge-Kutta method to second order ODEs

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, ...
user avatar
17 votes

Looking for Runge-Kutta 8th order in C/C++

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. ...
user avatar
15 votes
Accepted

C++ Best practices for dealing with many constants, variables in scientific codes

If you have constants that will not change before runs, declare them in a header file: ...
user avatar
  • 3,111
15 votes
Accepted

What are new c++20 features that are relevant to scientific computation?

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. ...
user avatar
  • 2,001
14 votes
Accepted

Super C++ optimization of matrix multiplication with Armadillo

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 ...
user avatar
  • 11.4k
14 votes

C++ Best practices for dealing with many constants, variables in scientific codes

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: ...
user avatar
  • 3,673
13 votes
Accepted

Beating typical BLAS libraries matrix multiplication performance

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 ...
13 votes
Accepted

Whittaker-Shannon interpolation: Accuracy dies with speedup; can it be fixed?

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: ...
user avatar
  • 1,030
12 votes

Modern C++ in scientific computing?

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 ...
user avatar
  • 2,185
12 votes
Accepted

Why do libraries need hand-vectorized code instead of compiler auto vectorization

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 ...
user avatar
  • 666
12 votes
Accepted

How to document math formulations in scientific computing codes?

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 ...
user avatar
  • 8,382
11 votes

Parallel vs Serial Thomas Algorithm

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 ...
user avatar
10 votes

Modern C++ in scientific computing?

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 ...
user avatar
  • 1,000
10 votes

Are there any "light-weight" FEM packages around?

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 ...
user avatar
  • 1,877
10 votes
Accepted

Method to check for positive definite matrices

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 ...
user avatar
  • 2,139
9 votes

Fast, lightweight C++ tensor library for dimension-agnostic code

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 ...
user avatar
9 votes
Accepted

Looking for Runge-Kutta 8th order in C/C++

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 ...
user avatar
  • 2,473
9 votes
Accepted

Programming Finite Element Methods in C++

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....
user avatar
  • 435
9 votes

I am searching for a C++ code implementing the complex polygamma function

I implemented a complex polygamma function in the Julia standard library (https://github.com/JuliaLang/julia/pull/7125). It is possible to call this from C++ by linking to the libjulia library. The ...
user avatar
  • 211
9 votes

Recommended language/environment for large scale semi-continuous biological models

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. ...
user avatar
8 votes

Programming Finite Element Methods in C++

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 ...
user avatar
8 votes

Super C++ optimization of matrix multiplication with Armadillo

@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: ...
user avatar
8 votes
Accepted

Stabilizing a 3x3 real symmetric matrix eigenvalue calculation

This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a ...
user avatar
  • 11.4k
8 votes

Fastest method forfinding a solution to x*log(x)

This equation is not polynomial. Assuming both $K$ and $C$ are positive (as in your linked problem), then the solution of $C x \ln_2(x) - K = 0$ can be found in terms of the Lambert $W$ function or ...
user avatar
  • 763
8 votes

Method to check for positive definite matrices

Mostly, I'm leaving this answer here as a cautionary tale to not use a Choleski factorization. Most of the time, this is a fine answer. However, very specifically, it can and will fail, so if this ...
user avatar
  • 697
7 votes
Accepted

How to efficiently assemble global stiffness matrix in sparse storage format (c++)

When I do assembly, I use a Coordinate sparse matrix format. This is basically just a list of (row, col, value) tuples. This is especially useful for sparse matrix construction when the exact pattern ...
user avatar
  • 1,522
7 votes

Modern C++ in scientific computing?

The HPX library makes heavy use of a range of C++11 features such as move constructors and is also aiming to be a complete implementation of N4409 (Working Draft, Technical Specification for C++ ...
user avatar
  • 171
7 votes
Accepted

A good, simple book/resource on Parallel Programming in C++ for scientific computing

One of the first things that you need to understand about parallel programming is the difference between shared memory multiprocessor computer systems and distributed memory clusters. A shared ...
user avatar
7 votes

Matlab, Mathematica & LAPACK returning 3 different eigenvectors

Seems that you have a duplicate eigenvalue. Thus, you have two eigenpairs $(\lambda_1, x_1)$ and $(\lambda_2, x_2)$ where $\lambda_1 = \lambda_2$. Denote $\lambda = \lambda_1 = \lambda_2$. Let $\alpha$...
user avatar
  • 1,877

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