  Matrix diagonalization is almost always involved in computing the density
matrix needed in quantum chemistry calculations. In the case of modest matrix
sizes ($\lesssim$ 5000), performance of traditional dense diagonalization
algorithms on modern GPUs is underwhelming compared to the peak performance of
these devices. This motivates the exploration of alternative algorithms better
suited to these types of architectures. We newly derive, and present in detail,
an existing Chebyshev expansion algorithm [W. Liang et al, J. Chem. Phys. 2003]
whose number of required matrix multiplications scales with the square root of
the number of terms in the expansion. Focusing on dense matrices of modest
size, our implementation on GPUs results in large speed ups when compared to
diagonalization. Additionally, we improve upon this existing method by
capitalizing on the inherent task parallelism and concurrency in the algorithm.
This improvement is implemented on GPUs by using CUDA and HIP streams via the
MAGMA library and leads to a significant speed up over the serial-only approach
for smaller ($\lesssim$ 1000) matrix sizes. Lastly, we apply our technique to a
model system with a high density of states around the Fermi level which
typically presents significant challenges.
