Our general interest is in self-consistent-field (scf) theories of disordered fermions. They generate physically relevant subensembles ("scf ensembles") within a given Altland-Zirnbauer class. We are motivated to investigate such ensembles (i) by the possibility to discover new fixed points due to (long-range) interactions; (ii) by analytical scf theories that rely on partial self-consistency approximations awaiting a numerical validation; and (iii) by the overall importance of scf theories for the understanding of complex interaction-mediated phenomena in terms of effective single-particle pictures. In this paper we present an efficient, parallelized implementation solving scf problems with spatially local fields by applying a kernel-polynomial approach. Our first application is the Boguliubov-deGennes theory of the attractive-U Hubbard model in the presence of on-site disorder; the sc fields are the particle density n(r) and the gap function Delta(r). For this case, we reach system sizes unprecedented in earlier work. They allow us to study phenomena emerging at scales substantially larger than the lattice constant, such as the interplay of multifractality and interactions or the formation of superconducting islands. For example, we observe that the coherence length exhibits a nonmonotonic behavior with increasing disorder strength already at moderate U. With respect to methodology our results are important because we establish that partial self-consistency ("energy-only") schemes as typically employed in analytical approaches tend to miss qualitative physics such as island formation.