So the only way to find uniform estimate of partial sum of Fourier coefficients is to go through Hecke eigenvalue relation with Fourier coefficient? Is there any way to find the estimate directly, without using Ramanujan on average?

I am sorry, but I am really confused. Equation (8.15) of Iwaniec's 'Spectral Methods of Automorphic forms' does not have $r^\epsilon$ term. It is, in our language,$$N^2\mathrm{res}_{s=1} L(s,\phi\times\phi)$$ and the residue is of size $e^{\pi r}$.

Thank you very much for a detailed answer. One confusion is: how did you get $r^\epsilon$ term in (*)? Don't Rankin-Selberg and Hoffstein-Lockhart give $e^{\pi r}N$?

From the spectral decomposition $$L^2(\Gamma\backslash G)=L^2_{cusp}\oplus \mathbb{C}\oplus L_{cont}^2,$$ Any function satisfying (2),(3) and (5) should be in the continuous spectrum. Therefore it can be described by the given Eisenstein series (as it has only one cusp at $\infty$. (2), (3) and (5) imply (1) and (4) with @GHfromMO's correction.

But, $\Delta(E(z,s)+c)=s(1-s)E(z,s)\neq s(1-s)(E(z,s)+c)$

To avoid confusion, in which normalization is the above result? I mean are you using $\lambda(p)=\sqrt{p}a(p)$ where $a(n)$ are the coefficients of the Maass form.