  Computing the Pfaffian of a skew-symmetric matrix is a problem that arises in
various fields of physics. Both computing the Pfaffian and a related problem,
computing the canonical form of a skew-symmetric matrix under unitary
congruence, can be solved easily once the skew-symmetric matrix has been
reduced to skew-symmetric tridiagonal form. We develop efficient numerical
methods for computing this tridiagonal form based on Gauss transformations,
using a skew-symmetric, blocked form of the Parlett-Reid algorithm, or based on
unitary transformations, using block Householder transformations and Givens
rotations, that are applicable to dense and banded matrices, respectively. We
also give a complete and fully optimized implementation of these algorithms in
Fortran, and also provide Python, Matlab and Mathematica implementations for
convenience. Finally, we apply these methods to compute the topological charge
of a class D nanowire, and show numerically the equivalence of definitions
based on the Hamiltonian and the scattering matrix.
