In phase retrieval, the goal is to recover a signal $\mathbf{x}\in\mathbb{C}^N$ from the magnitudes of linear measurements $\mathbf{Ax}\in\mathbb{C}^M$. While recent theory has established that $M\approx 4N$ intensity measurements are necessary and sufficient to recover generic $\mathbf{x}$, there is great interest in reducing the number of measurements through the exploitation of sparse $\mathbf{x}$, which is known as compressive phase retrieval. In this work, we detail a novel, probabilistic approach to compressive phase retrieval based on the generalized approximate message passing (GAMP) algorithm. We then present a numerical study of the proposed PR-GAMP algorithm, demonstrating its excellent phase-transition behavior, robustness to noise, and runtime. Our experiments suggest that approximately $M\geq 2K\log_2(N/K)$ intensity measurements suffice to recover $K$-sparse Bernoulli-Gaussian signals for $\mathbf{A}$ with i.i.d Gaussian entries and $K\ll N$. Meanwhile, when recovering a 6k-sparse 65k-pixel grayscale image from 32k randomly masked and blurred Fourier intensity measurements at 30~dB measurement SNR, PR-GAMP achieved an output SNR of no less than 28~dB in all of 100 random trials, with a median runtime of only 7.3 seconds. Compared to the recently proposed CPRL, sparse-Fienup, and GESPAR algorithms, our experiments suggest that PR-GAMP has a superior phase transition and orders-of-magnitude faster runtimes as the sparsity and problem dimensions increase.