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 HN:ΣN→R on the Hamming hypercube ΣN={±1}N, which are defined by the property that {HN(σ)}σ∈ΣN is a centered Gaussian process with covariance E{HN(σ1)HN(σ2)} depending only on the scalar product ⟨σ1,σ2⟩. The asymptotic value of the optimum maxσ∈ΣNHN(σ) 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 σ 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 HN, and characterize the typical energy value it achieves. When the p-spin model HN satisfies a certain no-overlap gap assumption, for any ε>0, the algorithm outputs σ∈ΣN such that HN(σ)≥(1−ε)maxσ′HN(σ′), with high probability. The number of iterations is bounded in N and depends uniquely on ε. 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.