A linear response function (LRF) determines the mean-response of a nonlinear climate system to weak imposed forcings, and an eddy flux matrix (EFM) determines the eddy momentum and heat flux responses to mean-flow changes. Neither LRF nor EFM can be calculated from first principles due to the lack of a complete theory for turbulent eddies. Here the LRF and EFM for an idealized dry atmosphere are computed by applying numerous localized weak forcings, one at a time, to a GCM with Held-Suarez physics and calculating the mean-responses. The LRF and EFM for zonally-averaged responses are then constructed using these forcings and responses through matrix inversion. Tests demonstrate that LRF and EFM are fairly accurate. Spectral analysis of the LRF shows that the most excitable dynamical mode, the neutral vector, strongly resembles the model's Annular Mode. The framework described here can be employed to compute the LRF/EFM for zonally-asymmetric responses and more complex GCMs. The potential applications of the LRF/EFM constructed here are i) forcing a specified mean-flow for hypothesis-testing, ii) isolating/quantifying the eddyfeedbacks in complex eddy-mean flow interaction problems, and iii) evaluating/improving more generallyapplicable methods currently used to construct LRFs or diagnose eddy-feedbacks in comprehensive GCMs or observations. As an example for iii, in Part 2, the LRF is also computed using the fluctuation-dissipation theorem (FDT), and the previously-calculated LRF is exploited to investigate why FDT performs poorly in some cases. It is shown that dimension-reduction using leading EOFs, which is commonly used to construct LRFs from the FDT, can significantly degrade the accuracy due to the non-normality of the operator.