We propose and analyze several multilevel algorithms for the fast simulation of possibly nonstationary Gaussian random fields (GRFs) indexed, for example, by the closure of a bounded domain 𝒟⊂ℝn or, more generally, by a compact metric space 𝒳 such as a compact n-manifold ℳ. A colored GRF 𝒵, admissible for our algorithms, solves the stochastic fractional-order equation 𝒜β𝒵=𝒲 for some β>n/4, where 𝒜 is a linear, local, second-order elliptic self-adjoint differential operator in divergence form and 𝒲 is white noise on 𝒳. We thus consider GRFs on 𝒳 with covariance operators of the form 𝒞=𝒜−2β. The proposed algorithms numerically approximate samples of 𝒵 on nested sequences {𝒯ℓ}ℓ≥0 of regular, simplicial partitions 𝒯ℓ of 𝒟 and ℳ, respectively. Work and memory to compute one approximate realization of the GRF 𝒵 on the triangulation 𝒯ℓ of 𝒳 with consistency 𝒪(N−ρℓ), for some consistency order ρ>0, scale essentially linearly in Nℓ=#(𝒯ℓ), independent of the possibly low regularity of the GRF. The algorithms are based on a sinc quadrature for an integral representation of (the application of) the negative fractional-order elliptic “coloring” operator 𝒜−β to white noise 𝒲. For the proposed numerical approximation, we prove bounds of the computational cost and the consistency error in various norms.