  We demonstrate how to compute the low energy spectrum for small ($N\le 50$),
but otherwise arbitrary, spin-glass instances using modern Graphics Processing
Units or similar heterogeneous architecture. Our algorithm performs an
exhaustive (i.e., brute-force) search of all possible configurations to select
$S\ll 2^N$ lowest ones together with their corresponding energies. We mainly
focus on the Ising model defined on an arbitrary graph. An open-source
implementation based on CUDA Fortran and a suitable Python wrapper are
provided. As opposed to heuristic approaches, ours is exact and thus can serve
as a references point to benchmark other algorithms and hardware, including
quantum and digital annealers. Our implementation offers unprecedented speed
and efficiency already visible on commodity hardware. At the same time, it can
be easily launched on professional, high-end graphics cards virtually at no
extra effort. As a practical application, we employ it to demonstrate that the
recent Matrix Product State based algorithm-despite its one-dimensional
nature-can still accurately approximate the low energy spectrum of fully
connected graphs of size $N$ approaching $50$.
