This question is not to be confused with the similarly titled question here.

In the above lined question, I gave a complete answer, but noticed that things are apparently not so simple in the elliptic case. I am wondering if Theorem 3.1 in Gilbarg and Trudinger is true without modification on a Riemannian manifold.

In the parabolic case we can use $e^{-t}$ to get a strict subsolution and this carries over nicely to manifolds. But the $e^{\gamma x_1}$ in the elliptic case does not work so nicely. Indeed, the elliptic proof seems to work with so few hypotheses because $\nabla^2x_1=0$. In general, we have $$Le^{\gamma f}=(\gamma a^{ij}\nabla_{ij}f+\gamma^2a^{ij}\nabla_if\nabla_if+\gamma b^i\nabla_if)e^{\gamma f}.$$
So the question is which $f$ to choose. A general manifold definitely does not have a nonzero function with $\nabla_{ij}f=0$ globally, and it would be surprising if such a function existed for a given bounded domain. So we have to estimate this term.

When $\gamma\ge 0$, we can estimate $\gamma a^{ij}\nabla_{ij}f\ge C \gamma |\nabla^2 f|\lambda$. A natural candidate for $f$ is $r$, the distance to some point in the domain. Then we have $$Le^{\gamma r}\ge\lambda(\gamma^2-C\gamma)e^{\gamma r}$$ under the hypotheses of GT Theorem 3.1. This formula holds at smooth points, but the cut locus causes issues. The issue is that when $\gamma\ge 0$, Calabi's trick (Schoen--Yau, Lectures on Diff Geo, p. 21) does not work. If we choose $\gamma<0$ then this does work, but $\gamma a^{ij}\nabla_{ij}f\ge C \gamma |\nabla^2 f|\lambda$ is no longer true. One needs to replace $\lambda$ with $\Lambda$, which is $\le C\lambda$ for uniform ellipticity. On the other hand, strict ellipticity is also sufficient because then one can ignore the $\lambda$ terms by choosing $\gamma$ even larger.

So how does one prove the analogue of Theorem 3.1 in Gilbarg--Trudinger without making stronger assumptions on the ellipticity of $a$?

Note: $|\nabla^2 r|$ can be estimated via the Hessian comparison theorem, noting that the sectional curvature is in a compact interval over a bounded domain.