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 ...
Chris Rackauckas's 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: ...
Richard's user avatar
  • 3,606
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 ...
Kirill's user avatar
  • 11.4k
14 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 ...
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: ...
spektr's user avatar
  • 3,818
14 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. ...
user14717's user avatar
  • 2,085
13 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 ...
BrunoLevy's user avatar
  • 2,285
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: ...
njuffa's user avatar
  • 1,320
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 ...
ggael's 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 ...
Anton Menshov's user avatar
  • 8,572
11 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 ...
Tolga Birdal's user avatar
  • 2,179
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 ...
fred's user avatar
  • 1,000
10 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 ...
SGJ's user avatar
  • 221
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 ...
knl's user avatar
  • 2,041
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 ...
Damascus Steel's user avatar
9 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 ...
wyer33's user avatar
  • 717
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. ...
Chris Rackauckas's user avatar
9 votes
Accepted

Searching for recent code source for "Parallel scientific computing in C++ and MPI "

First of all, the book you mention is very old. In fact, it misses the last two MPI standards (3 and 4), and every C++ standard from C++11 on. Secondly, know that MPI has officially no C++ bindings, ...
Victor Eijkhout's 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: ...
Wolfgang Bangerth's 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 ...
Kirill's 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 ...
horchler's user avatar
  • 763
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++ ...
SteveL's 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 ...
Brian Borchers's 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$...
knl's user avatar
  • 2,041
7 votes

Physics Simulation in C++

I think you are missing a very important and crucial step that lies exactly between the physics and simulation: the mathematical model. In order to model any physics, one has to formulate the ...
Anton Menshov's user avatar
  • 8,572
7 votes
Accepted

How to justify using available code (in different language) for comparing algorithms

So, you are comparing a generally slower Matlab implementation of algorithm A to a generally faster C++ implementation of algorithm B, and still getting the advantage for A. I would say, ...
Anton Menshov's user avatar
  • 8,572
7 votes
Accepted

Compiled c++ code runs much faster with double than float. Explanation?

I have looked at your simple code example, and my suspicion is that what you observe in loss of speed is due to the C-heritage requirement that right-hand side expressions be evaluated using ...
hardmath's user avatar
  • 3,249
7 votes
Accepted

The Formula of Explicit Runge-Kutta Fourteen order

The 14th order methods due to Feagin can be found in DifferentialEquations.jl. An example using them with 128-bit floating point arithmetic is as follows: ...
Chris Rackauckas's user avatar
6 votes

In C++, how do you calculate the analytical value of $\int_a^b \left|\sin x \right|\,dx$?

You could just go for the direct final answer of the integral: $$\int_{a}^{b}\vert {\sin(x)}\vert dx = 2\left(\left\lfloor\frac{b}{\pi}\right\rfloor - \left\lfloor\frac{a}{\pi}\right\rfloor\right) + \...
SassyQuatch's user avatar
6 votes

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

XTensor is a modern approach and is getting more and more popular. https://github.com/QuantStack/xtensor
Jingpeng Wu's user avatar

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