The output resp. the prediction is measured with a loss function. Then what? How is the error used in the gradient descent approach?
If the output of your neural network is $\hat{y}_\theta$, where $\theta$ refers to the parameters of the net, and if your loss function is $\mathcal{L}(y,\hat{y}_\theta)$, where $y$ is the target, then the gradient descent update rule is
$$
\theta_{k+1} \leftarrow \theta_k - \alpha \cdot \nabla_\theta \mathcal{L}(y,\hat{y}_\theta),
$$
where $\alpha$ is the learning rate. The "error" and "loss" are the same in this context. You use the backpropagation algorithm to compute $\nabla_\theta \mathcal{L}(y,\hat{y}_\theta)$ exactly, since $\hat{y}_\theta$ is a complicated function of $\theta$.
So I guess the final loss function takes into account all these parameters.
Yes these parameters are in $\theta$.
Will gradient descent then optimize this function (with all parameters)?
Yes. You can think of the gradient $\nabla_\theta\mathcal{L}(y,\hat{y}_\theta)$ as
$$
\nabla_\theta\mathcal{L}(y,\hat{y}_\theta) = \frac{\partial \mathcal{L}(y,\hat{y}_\theta)}{\partial \theta}.
$$
Then, the gradient can be computed using matrix calculus. Note that if the parameters are stored in matrices and vectors, the gradient can be computed individually for each of these, but the update has to happen to all parameters simultaneously.
So it would be a very high-dimensional function to minimize?
Yes, very.