# Tag Info

### 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 ...
• 11.6k
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, ...
• 4,521

### 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. ...
• 29.8k
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: ...
• 3,111
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. ...
• 2,001
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 ...
• 11.4k

### 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: ...
• 3,673
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 ...
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: ...
• 1,030

### 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 ...
• 2,185
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 ...
• 666
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 ...
• 8,382

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

### 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 ...
• 1,000

### 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 ...
• 1,877
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 ...
• 2,139

### 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 ...
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 ...
• 2,473
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....
• 435

### 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 ...
• 211

### 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. ...
• 11.6k

### 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 ...
• 50.7k

### 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: ...
• 50.7k
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 ...
• 11.4k

### 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 ...
• 763

### 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 ...
• 697
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 ...
• 1,522

### 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++ ...
• 171
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 ...
• 17.6k
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$...