Most Popular
1500 questions
19
votes
6 answers
To what extent is generic and meta-programming using C++ templates useful in computational science?
The C++ language provides generic programming and metaprogramming through templates. These techniques have found their way into many large-scale scientific computing packages (e.g., MPQC, LAMMPS, CGAL, Trilinos). But what have they actually…
Deathbreath
- 1,042
- 7
- 20
19
votes
2 answers
How to determine if a numerical solution to a PDE is converging to a continuum solution?
The Lax equivalence theorem states that consistency and stability of a numerical scheme for a linear initial value problem is a necessary and sufficient condition for convergence. But for nonlinear problems, numerical methods can converge very…
Jed Brown
- 25,650
- 3
- 72
- 130
19
votes
4 answers
Is Fortuna or Mersenne Twister preferable as an algorithmic RNG?
A recent answer mentioned the use of Fortuna or Mersenne Twister Random Number Generators (RNGs) to seed a Monte Carlo simulation. I hadn't heard of Fortuna before so I looked it up - looks like it is mainly intended for cryptographic use.
I…
winwaed
- 658
- 5
- 13
19
votes
3 answers
Is it well known that some optimization problems are equivalent to time-stepping?
Given a desired state $y_0$ and a regularization parameter $\beta \in \mathbb R$, consider the problem of finding a state $y$ and a control $u$ to minimize a functional
\begin{equation}
\frac{1}{2} \| y - y_0 \|^2 + \frac{\beta}{2} \| u…
Andrew T. Barker
- 677
- 3
- 10
18
votes
6 answers
How to reorder variables to produce a banded matrix of minimum bandwidth?
I'm trying to solve a 2D Poisson equation by finite differences. In the process, I obtain a sparse matrix with only $5$ variables in each equation. For example, if the variables were $U$, then the discretization would yield:
$$U_{i-1,j} + U_{i+1,j}…
Paul
- 12,045
- 7
- 56
- 129
18
votes
1 answer
Intuitive motivation for BFGS update
I am teaching a numerical analysis survey class and am seeking motivation for the BFGS method for students with limited background/intuition in optimization!
While I don't have time to prove rigorously that everything converges, I'm looking to give…
Justin Solomon
- 1,341
- 7
- 24
18
votes
3 answers
Solving unconstrained nonlinear optimization problems on GPU
I am trying to solve some unconstrained nonlinear optimization problems on GPU (CUDA).
The objective function is a smooth nonlinear function, and its gradient is relatively cheap to compute analytically, so I don't need to bother with numerical…
user0002128
- 341
- 2
- 4
18
votes
4 answers
Which methods can ensure that physical quantities remain positive throughout a PDE simulation?
Physical quantities like pressure, density, energy, temperature, and concentration should always be positive, but numerical methods sometimes compute negative values during the solution process. This is not okay because the equations will compute…
Jed Brown
- 25,650
- 3
- 72
- 130
18
votes
4 answers
What are some good strategies to test a floating point arithmetic implementation for double numbers?
For IEEE, the single representation is 1-bit sign, 8-bit exponent and 23-bit mantissa. This means that at each exponent value, you can test all 2^23-1 (roughly 9mil cases) possible combination of binary representation (give or take). Then you do it…
Quang Thinh Ha
- 475
- 1
- 4
- 8
18
votes
5 answers
What is the best way to determine the number of non zeros in sparse matrix multiplication?
I was wondering whether there is a fast and efficient method to find the number of non zeros in advance for sparse matrix multiplication operation assuming both matrices are in CSC or CSR format.
I know there is one in smmp package but I need…
Recker
- 333
- 2
- 7
18
votes
2 answers
Is there an efficient algorithm for matrix-valued continued fractions?
Suppose I have a matrix equation recursively defined as
A[n] = inverse([1 - b[n]A[n+1]]) * a[n]
Then the equation for A[1] looks similar to a continued fraction, for which there are some highly efficient methods that avoid tedious recalculation…
Lagerbaer
- 487
- 4
- 9
18
votes
2 answers
Complexity of matrix inversion in numpy
I am solving differential equations that require to invert dense square matrices. This matrix inversion consumes the most of my computation time, so I was wondering if I am using the fastest algorithm available.
My current choice is…
physicsGuy
- 409
- 1
- 4
- 12
18
votes
1 answer
How to Run MPI-3.0 in shared memory mode like OpenMP
I am parallelizing code to numerically solve a 5 Dimensional population balance model. Currently I have a very good MPICH2 parallelized code in FORTRAN but as we increase parameter values the arrays become too large to run in distributed memory…
Franklin Betten
- 191
- 4
18
votes
4 answers
Portable multicore/NUMA memory allocation/initialization best practices
When memory bandwidth limited computations are performed in shared memory environments (e.g. threaded via OpenMP, Pthreads, or TBB), there is a dilemma of how to ensure that the memory is correctly distributed across physical memory, such that each…
Jed Brown
- 25,650
- 3
- 72
- 130
18
votes
3 answers
What programming strategies can I take for easily modifying algorithm parameters?
Developing scientific algorithms is a highly iterative process often involving changing lots of parameters that I will want to vary either as part of my experimental design or as part of tweaking algorithm performance. What strategies can I take for…
Scottie T
- 289
- 1
- 4