  We present algorithmic and implementation details for the fully
self-consistent finite-temperature $GW$ method in Gaussian Bloch orbitals for
solids. Our implementation is based on the finite-temperature Green's function
formalism in which all equations are solved on the imaginary axis, without
resorting to analytical continuation during the self-consistency. No
quasiparticle approximation is employed and all matrix elements of the
self-energy are explicitly evaluated. The method is tested by evaluating the
band gaps of selected semiconductors and insulators. We show agreement with
other, differently formulated finite-temperature sc$GW$ implementations when
finite-size corrections and basis set errors are taken into account. By
migrating computationally intensive calculations to GPUs, we obtain scalable
results on large supercomputers with nearly optimal performance. Our work
demonstrates the applicability of Gaussian orbital based sc$GW$ for $\emph{ab
initio}$ correlated materials simulations and provides a sound starting point
for embedding methods built on top of $GW$.
