  We develop a theory of the local density of states (LDOS) of disordered
superconductors, employing the non-linear sigma-model formalism and the
renormalization-group framework. The theory takes into account the interplay of
disorder and interaction couplings in all channels, treating the systems with
short-range and Coulomb interactions on equal footing. We explore 2D systems
that would be Anderson insulators in the absence of interaction and 2D or 3D
systems that undergo Anderson transition in the absence of interaction. We
evaluate both the average tunneling density of states and its mesoscopic
fluctuations which are related to the LDOS multifractality in normal disordered
systems. The obtained average LDOS shows a pronounced depletion around the
Fermi energy, both in the metallic phase (i.e., above the superconducting
critical temperature $T_c$) and in the insulating phase near the
superconductor-insulator transition (SIT). The fluctuations of the LDOS are
found to be particularly strong for the case of short-range interactions --
especially, in the regime when $T_c$ is enhanced by Anderson localization. On
the other hand, the long-range Coulomb repulsion reduces the mesoscopic LDOS
fluctuations. However, also in a model with Coulomb interaction, the
fluctuations become strong when the systems approaches the SIT.
