  In order to study structure of proto-neutron stars and those in subsequent
cooling stages, it is of great interest to calculate inhomogeneous hot and cold
nuclear matter in a variety of phases. The finite-temperature
Hartree-Fock-Bogoliubov (FT-HFB) theory is a primary choice for this purpose,
however, its numerical calculation for superfluid (superconducting)
many-fermion systems in three dimensions requires enormous computational costs.
To study a variety of phases in the crust of hot and cold neutron stars, we
propose an efficient method to perform the FT-HFB calculation with the
three-dimensional (3D) coordinate-space representation. Recently, an efficient
method based on the contour integral of Green's function with the shifted
conjugate-orthogonal conjugate-gradient method has been proposed [Phys. Rev. C
95, 044302 (2017)]. We extend the method to the finite temperature, using the
shifted conjugate-orthogonal conjugate-residual method. We benchmark the 3D
coordinate-space solver of the FT-HFB calculation for hot isolated nuclei and
fcc phase in the inner crust of neutron stars at finite temperature. The
computational performance of the present method is demonstrated. Different
critical temperatures of the quadrupole and the octupole deformations are
confirmed for $^{146}$Ba. The robustness of the shape coexistence feature in
$^{184}$Hg is examined. For the neutron-star crust, the deformed neutron-rich
Se nuclei embedded in the sea of superfluid low-density neutrons appear in the
fcc phase at the nucleon density of 0.045 fm$^{-3}$ and the temperature of $k_B
T=200$ keV. The efficiency of the developed solver is demonstrated for nuclei
and inhomogeneous nuclear matter at finite temperature. It may provide a
standard tool for nuclear physics, especially for the structure of the hot and
cold neutron-star matters.
