Mean-field spin glasses are families of random energy functions (Hamiltonians) on high-dimensional product spaces. In this paper we consider the case of Ising mixed $p$-spin models, namely Hamiltonians $H_N:Sigma_Nto {mathbb R}$ on the Hamming hypercube $Sigma_N = {pm 1}^N$, which are defined by the property that ${H_N({boldsymbol sigma})}_{{boldsymbol sigma}in Sigma_N}$ is a centered Gaussian process with covariance ${mathbb E}{H_N({boldsymbol sigma}_1)H_N({boldsymbol sigma}_2)}$ depending only on the scalar product $langle {boldsymbol sigma}_1,{boldsymbol sigma}_2rangle$. The asymptotic value of the optimum $max_{{boldsymbol sigma}in Sigma_N}H_N({boldsymbol sigma})$ was characterized in terms of a variational principle known as the Parisi formula, first proved by Talagrand and, in a more general setting, by Panchenko. The structure of superlevel sets is extremely rich and has been studied by a number of authors. Here we ask whether a near optimal configuration ${boldsymbol sigma}$ can be computed in polynomial time. We develop a message passing algorithm whose complexity per-iteration is of the same order as the complexity of evaluating the gradient of $H_N$, and characterize the typical energy value it achieves. When the $p$-spin model $H_N$ satisfies a certain no-overlap gap assumption, for any $varepsilon>0$, the algorithm outputs ${boldsymbol sigma}inSigma_N$ such that $H_N({boldsymbol sigma})ge (1-varepsilon)max_{{boldsymbol sigma}} H_N({boldsymbol sigma})$, with high probability. The number of iterations is bounded in $N$ and depends uniquely on $varepsilon$. More generally, regardless of whether the no-overlap gap assumption holds, the energy achieved is given by an extended variational principle, which generalizes the Parisi formula.