In a single application, it is rather likely that you only will need a small subset of all the possible extremes of the generalized hypergeometric function. It is a very general function, after all. Having an idea about the range of $z$ and of the parameters $a_i, b_i$ would allow giving more specific advice.
In general, the standard method, assuming $p \le q + 1$, is of course to use the defining power series when $|z|$ is small. If $p < q + 1$, it is best to switch to an asymptotic expansion when $|z|$ is large, either because the Taylor series converges too slowly and/or because it becomes too inaccurate due to catastrophic cancellation. The best cutoff between these algorithms depends on the parameters and the accuracy requirements.
For ${}_1F_2$ the asymptotic series is given by http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F2/06/02/03/ It looks rather horrible, but if your $a_1, b_1, b_2$ are fixed, you can compute numerical values for the coefficients in advance. General formulas are found in the DLMF: http://dlmf.nist.gov/16.11 (Note that some care is required to select the correct branch cuts.)
If there is a range where neither the Taylor series nor the asymptotic series works well enough, "exponentially-improved expansions" might be useful. Another possibility worth mentioning is that you could just plug the hypergeometric differential equation into a general-purpose ODE solver. This should work quite well especially if you only need 4-5 digits. This can be used to do analytic continuation from a small $z$ (where the power series works well) to a larger one, or in reverse from a value obtained via an asymptotic series (you might need to do a bit more work to get all the derivatives needed as initial values).
If you need functions with $p = q + 1$ on the whole complex plane then $1/z$ transformation formulas can be used to map the exterior of the unit disk to the interior. Some convergence acceleration algorithm or other method, like numerical integration of the ODE, must be used close to the unit circle. If $p > q + 1$ the radius is convergence is zero, so if the function you want to evaluate is given by such a divergent series, you might need to apply a Borel transform (numerically or symbolically) to reduce it to a convergent series.
For a complete implementation, there are other issues to consider as well (for example, dealing with parameters which are extremely large or very close to negative integers). For sufficiently bad parameters, it will be very difficult to obtain accurate values with double precision, no matter what you do, so arbitrary-precision arithmetic might be needed.
I should note that I have written a nearly complete numerical implementation of the generalized hypergeometric function for the mpmath library (it is currently missing asymptotic series for functions higher than ${}_2F_3$), which might be useful to study or run tests against (assuming it is not fast enough already for your purposes).
