We provide an end-to-end exploration of a distinct modified gravitational theory in Jordan-Brans-Dicke (JBD) gravity, from an analytical and numerical description of the background expansion and linear perturbations, to the nonlinear regime captured with a hybrid suite of N-body simulations, to the cosmological constraints from existing probes of the expansion history, the large-scale structure, and the cosmic microwave background (CMB). We have focused on JBD gravity as it both approximates a wider class of Horndeski scalar-tensor theories on cosmological scales and allows us to adequately model the nonlinear corrections to the matter power spectrum. In a combined analysis of the Planck 2018 CMB temperature, polarization, and lensing reconstruction, together with Pantheon supernova distances and the Baryon Oscillation Spectroscopic Survey (BOSS) measurements of baryon acoustic oscillation distances, the Alcock-Paczynski effect, and the growth rate, we constrain the JBD coupling constant to ωBD>970 (95% confidence level; C.L.) in agreement with the General Relativistic expectation given by ωBD→∞. In the unrestricted JBD model, where the effective gravitational constant at present, Gmatter/G, is additionally varied, increased dataset concordance (e.g., within 1σ agreement in S8=σ8Ωm/0.3) enables us to further include the combined ("3×2pt") dataset of cosmic shear, galaxy-galaxy lensing, and overlapping redshift-space galaxy clustering from the Kilo Degree Survey and the 2-degree Field Lensing Survey (KiDS×2dFLenS). In analyzing the weak lensing measurements, the nonlinear corrections due to baryons, massive neutrinos, and modified gravity are simultaneously modeled and propagated in the cosmological analysis for the first time. In the joint analysis of all datasets, we constrain ωBD>1540 (95% C.L.), Gmatter/G=0.997±0.029, the sum of neutrino masses, Σmν<0.12 eV (95% C.L.), and the baryonic feedback amplitude, B < 2.8 (95% CL), all in agreement with the standard model expectation. In fixing the sum of neutrino masses, the lower bound on the coupling constant strengthens to ωBD>1460 and ωBD>2230 (both at 95% C.L.) in the restricted and unrestricted JBD models, respectively. We explore the impact of the JBD modeling choices, and show that a more restrictive parametrization of the coupling constant degrades the neutrino mass bound by up to a factor of three. In addition to the improved concordance between KiDS×2dFLenS and Planck, the tension in the Hubble constant between Planck and the direct measurement of Riess et al. (2019) is reduced to ∼3σ; however, we find no substantial model selection preference for JBD gravity relative to ΛCDM. We further show that a positive shift in the effective gravitational constant suppresses the CMB damping tail, which might complicate future inferences of small-scale physics, given its degeneracy with the primordial helium abundance, the effective number of neutrinos, and the running of the spectral index.