A novel very efficient algorithm based on the Geometrical Optics (GO) is presented for the analysis of graded index (GRIN) lens antennas, namely dielectric inhomogeneous lenses with a three-dimensional arbitrary varying refractive index. A family of curved ray-paths are traced starting from a set of points defined on the lens input interface, which is illuminated by a feed antenna, up to a corresponding set of points on the output interface; i.e., the lens radiating aperture. The ray-tracing is numerically performed in combination with the field transportation along the ray by exploiting an additional wavefront-curvature transport equation, thus providing a single independent vector Ordinary Differential Equation (ODE) for each ray. This scheme allows a complete parallelization of the algorithm by assigning any ray to a different thread. Crossing the lens input/output interfaces at the two end-points of the curved ray is accomplished by augmenting the refraction Snell's law and the Fresnel transmission coefficients with a novel compact closed form dyadic formula for the transmitted wavefront curvature. The ODE output provides the field aperture distribution on the output lens interface, permitting the far field calculation as a radiation integral. To this end, an ad-hoc quadrature rule is exploited, which requires an aperture sampling rate corresponding to the slowly varying part of the integrand, while the rapid phase variation is accounted for analytically; thus resulting in an extremely efficient and frequency independent scheme. The effectiveness and accuracy of the algorithm is shown by resorting to a couple of well-known analytical benchmarks and to two more general examples with real-life antennas: a Spaceborne Weather Radar lens antenna and a feed antenna with a zero-focal lens.