I work on Kokkos Kernels, which is performance-portable math libraries for the Kokkos programming system. We implement BLAS, sparse matrix operations, a variety of solvers, some ODE stuff, and various utilities for a variety of serial and parallel execution modes on x86, ARM, and the three major GPU vendors.
The simplest and highest-impact tests are all the edge cases - if an input matrix/vector/scalar is 0/1/-1/NaN, that usually tells you a lot about what the outputs should be.
One of the things that's been thorniest for us is determining a sensible numerical limit for error in our algorithms. The simplest example is a dot product - summing floats is not associative, so doing it in parallel is not bitwise the same as serial. For dot in particular it's relatively easy to come up with an error bound, but for anything more complicated it takes a particular expertise that is not always available. This is certainly a work in progress for us, and sometimes (usually) we just pick a magic tolerance out of thin air that seems to work.
We do test solvers using analytical solutions.
We also test solvers by inverting them, e.g. if we're solving Ax = y, for x, then Ax should come out "close" to the original y (see error tolerance discussion above).
We also test many things against specific regressions we've goofed up before.
One of the most interesting things to me is that our test suite has identified many bugs in vendor math libraries (OpenBLAS, MKL, cuSparse, rocSparse, etc.) - a major component of what we do is wrap up these vendor libraries in a common interface so our users don't have to do any work when they switch supercomputers, so in practice we test them all pretty thoroughly as well.
All good stuff, I’d add though that for many numerical routines you end up testing on simple well known systems that have well defined analytic answers.
So for e.g. if you’re writing solvers for ODE systems, then you tend to test it on things like simple harmonic oscillators.
If you’re very lucky, some algorithms have bounded errors. For e.g. the fast multipole method for evaluating long range forces does, so you can check this for very simple systems.
If you have a numerical routine that converts from one representation to another, you can test it by calling the routine implementing the inverse conversion.
I work on Kokkos Kernels, which is performance-portable math libraries for the Kokkos programming system. We implement BLAS, sparse matrix operations, a variety of solvers, some ODE stuff, and various utilities for a variety of serial and parallel execution modes on x86, ARM, and the three major GPU vendors.
The simplest and highest-impact tests are all the edge cases - if an input matrix/vector/scalar is 0/1/-1/NaN, that usually tells you a lot about what the outputs should be.
One of the things that's been thorniest for us is determining a sensible numerical limit for error in our algorithms. The simplest example is a dot product - summing floats is not associative, so doing it in parallel is not bitwise the same as serial. For dot in particular it's relatively easy to come up with an error bound, but for anything more complicated it takes a particular expertise that is not always available. This is certainly a work in progress for us, and sometimes (usually) we just pick a magic tolerance out of thin air that seems to work.
We do test solvers using analytical solutions.
We also test solvers by inverting them, e.g. if we're solving Ax = y, for x, then Ax should come out "close" to the original y (see error tolerance discussion above).
We also test many things against specific regressions we've goofed up before.
One of the most interesting things to me is that our test suite has identified many bugs in vendor math libraries (OpenBLAS, MKL, cuSparse, rocSparse, etc.) - a major component of what we do is wrap up these vendor libraries in a common interface so our users don't have to do any work when they switch supercomputers, so in practice we test them all pretty thoroughly as well.
All good stuff, I’d add though that for many numerical routines you end up testing on simple well known systems that have well defined analytic answers.
So for e.g. if you’re writing solvers for ODE systems, then you tend to test it on things like simple harmonic oscillators.
If you’re very lucky, some algorithms have bounded errors. For e.g. the fast multipole method for evaluating long range forces does, so you can check this for very simple systems.
If you have a numerical routine that converts from one representation to another, you can test it by calling the routine implementing the inverse conversion.