A Riemannian gradient descent algorithm and a truncated variant will be presented for solving systems of phaseless equations |Ax|^2 = y. The algorithms are developed by exploiting the inherent low rank structure of the problem based on the embedded manifold of rank-1 positive semidefinite matrices. Theoretical recovery guarantee has been established for the truncated variant, showing that the algorithm is able to achieve successful recovery when the number of equations is proportional to the number of unknowns. In addition, we will present a loss function without spurious local minima when the sampling complexity is optimal.