  Standard full-shape clustering analyses in Fourier space rely on a fixed
power spectrum template, defined at the fiducial cosmology used to convert
redshifts into distances, and compress the cosmological information into the
Alcock-Paczynski parameters and the linear growth rate of structure. In this
paper, we propose an analysis method that operates directly in the cosmology
parameter space and varies the power spectrum template accordingly at each
tested point. Predictions for the power spectrum multipoles from the TNS model
are computed at different cosmologies in the framework of $\Lambda \rm{CDM}$.
Applied to the final eBOSS QSO and LRG samples together with the low-z DR12
BOSS galaxy sample, our analysis results in a set of constraints on the
cosmological parameters $\Omega_{\rm cdm}$, $H_0$, $\sigma_8$, $\Omega_{\rm b}$
and $n_s$. To reduce the number of computed models, we construct an iterative
process to sample the likelihood surface, where each iteration consists of a
Gaussian process regression. This method is validated with mocks from N-body
simulations. From the combined analysis of the (e)BOSS data, we obtain the
following constraints: $\sigma_8=0.877\pm 0.049$ and $\Omega_{\rm
m}=0.304^{+0.016}_{-0.010}$ without any external prior. The eBOSS quasar sample
alone shows a $3.1\sigma$ discrepancy compared to the Planck prediction.
