diff --git a/.github/workflows/Aqua.yml b/.github/workflows/Aqua.yml
index 73e12cb6..1c73f8f5 100644
--- a/.github/workflows/Aqua.yml
+++ b/.github/workflows/Aqua.yml
@@ -13,6 +13,13 @@ jobs:
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/qp_block_data"),
+          ])
       - name: Aqua.jl
         run: |
           PKG_SRC_PATH=`pwd`
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index be0b8658..7505fff9 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -14,6 +14,13 @@ jobs:
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/qp_block_data"),
+          ])
       - name: Install dependencies
         run: julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 7b1d59ca..630355b4 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -54,6 +54,13 @@ jobs:
             ${{ runner.os }}-test-${{ env.cache-name }}-
             ${{ runner.os }}-test-
             ${{ runner.os }}-
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MathOptInterface", rev="bl/qp_block_data"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
diff --git a/src/moi_nlp_model.jl b/src/moi_nlp_model.jl
index a4975aa3..6bfd41e5 100644
--- a/src/moi_nlp_model.jl
+++ b/src/moi_nlp_model.jl
@@ -132,10 +132,7 @@ function NLPModels.cons_nln!(nlp::MathOptNLPModel, x::AbstractVector, c::Abstrac
   offset = 0
   if nlp.quadcon.nquad > 0
     offset += nlp.quadcon.nquad
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
-    end
+    MOI.eval_constraint(nlp.quadcon.evaluator, view(c, 1:(nlp.quadcon.nquad)), x)
   end
   if nlp.nlcon.nnln > 0
     offset += nlp.nlcon.nnln
@@ -166,11 +163,8 @@ function NLPModels.cons!(nlp::MathOptNLPModel, x::AbstractVector, c::AbstractVec
   end
   if nlp.quadcon.nquad > 0
     offset += nlp.quadcon.nquad
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      c[nlp.meta.nlin + i] =
-        0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
-    end
+    ind_quad = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
+    MOI.eval_constraint(nlp.quadcon.evaluator, view(c, ind_quad), x)
   end
   if nlp.nlcon.nnln > 0
     offset += nlp.nlcon.nnln
@@ -210,14 +204,10 @@ function NLPModels.jac_nln_structure!(
 )
   offset = 0
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
-      qcon = nlp.quadcon.constraints[i]
-      ind_quad = (offset + 1):(offset + qcon.nnzg)
-      view(rows, ind_quad) .= i
-      view(cols, ind_quad) .= qcon.g
-      offset += qcon.nnzg
-    end
+    ind_quad = 1:(nlp.quadcon.nnzj)
+    view(rows, ind_quad) .= nlp.quadcon.jac_rows
+    view(cols, ind_quad) .= nlp.quadcon.jac_cols
+    offset += nlp.quadcon.nnzj
   end
   @assert offset == nlp.quadcon.nnzj
   if nlp.nlcon.nnln > 0
@@ -257,14 +247,10 @@ function NLPModels.jac_structure!(
     offset += nlp.lincon.nnzj
   end
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
-      qcon = nlp.quadcon.constraints[i]
-      ind_quad = (offset + 1):(offset + qcon.nnzg)
-      view(rows, ind_quad) .= nlp.meta.nlin .+ i
-      view(cols, ind_quad) .= qcon.g
-      offset += qcon.nnzg
-    end
+    ind_quad = (offset + 1):(offset + nlp.quadcon.nnzj)
+    view(rows, ind_quad) .= nlp.meta.nlin .+ nlp.quadcon.jac_rows
+    view(cols, ind_quad) .= nlp.quadcon.jac_cols
+    offset += nlp.quadcon.nnzj
   end
   @assert offset == nlp.lincon.nnzj + nlp.quadcon.nnzj
   if nlp.nlcon.nnln > 0
@@ -304,26 +290,8 @@ function NLPModels.jac_nln_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals:
   offset = 0
   if nlp.quadcon.nquad > 0
     ind_quad = 1:(nlp.quadcon.nnzj)
-    view(vals, ind_quad) .= 0.0
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      for (j, ind) in enumerate(qcon.b.nzind)
-        k = qcon.dg[ind]
-        vals[offset + k] += qcon.b.nzval[j]
-      end
-      for j = 1:(qcon.nnzh)
-        row = qcon.A.rows[j]
-        col = qcon.A.cols[j]
-        val = qcon.A.vals[j]
-        k1 = qcon.dg[row]
-        vals[offset + k1] += val * x[col]
-        if row != col
-          k2 = qcon.dg[col]
-          vals[offset + k2] += val * x[row]
-        end
-      end
-      offset += qcon.nnzg
-    end
+    MOI.eval_constraint_jacobian(nlp.quadcon.evaluator, view(vals, ind_quad), x)
+    offset += nlp.quadcon.nnzj
   end
   if nlp.nlcon.nnln > 0
     ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
@@ -355,26 +323,8 @@ function NLPModels.jac_coord!(nlp::MathOptNLPModel, x::AbstractVector, vals::Abs
   end
   if nlp.quadcon.nquad > 0
     ind_quad = (nlp.lincon.nnzj + 1):(nlp.lincon.nnzj + nlp.quadcon.nnzj)
-    view(vals, ind_quad) .= 0.0
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      for (j, ind) in enumerate(qcon.b.nzind)
-        k = qcon.dg[ind]
-        vals[offset + k] += qcon.b.nzval[j]
-      end
-      for j = 1:(qcon.nnzh)
-        row = qcon.A.rows[j]
-        col = qcon.A.cols[j]
-        val = qcon.A.vals[j]
-        k1 = qcon.dg[row]
-        vals[offset + k1] += val * x[col]
-        if row != col
-          k2 = qcon.dg[col]
-          vals[offset + k2] += val * x[row]
-        end
-      end
-      offset += qcon.nnzg
-    end
+    MOI.eval_constraint_jacobian(nlp.quadcon.evaluator, view(vals, ind_quad), x)
+    offset += nlp.quadcon.nnzj
   end
   if nlp.nlcon.nnln > 0
     ind_nnln = (offset + 1):(offset + nlp.nlcon.nnzj)
@@ -422,11 +372,8 @@ function NLPModels.jprod_nln!(
 )
   increment!(nlp, :neval_jprod_nln)
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      # Jv[i] += (Aᵢ * x + bᵢ)ᵀ * v
-      qcon = nlp.quadcon.constraints[i]
-      Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
-    end
+    ind_quad = 1:(nlp.quadcon.nquad)
+    MOI.eval_constraint_jacobian_product(nlp.quadcon.evaluator, view(Jv, ind_quad), x, v)
   end
   if nlp.nlcon.nnln > 0
     ind_nnln = (nlp.quadcon.nquad + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln)
@@ -434,7 +381,7 @@ function NLPModels.jprod_nln!(
   end
   if nlp.oracles.ncon > 0
     for i =
-      (nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles.ncon)
+        (nlp.quadcon.nquad + nlp.nlcon.nnln + 1):(nlp.quadcon.nquad + nlp.nlcon.nnln + nlp.oracles.ncon)
 
       Jv[i] = 0
     end
@@ -478,12 +425,8 @@ function NLPModels.jprod!(
     )
   end
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
-      qcon = nlp.quadcon.constraints[i]
-      Jv[nlp.meta.nlin + i] =
-        coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
-    end
+    ind_quad = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
+    MOI.eval_constraint_jacobian_product(nlp.quadcon.evaluator, view(Jv, ind_quad), x, v)
   end
   if nlp.nlcon.nnln > 0
     ind_nnln =
@@ -544,12 +487,16 @@ function NLPModels.jtprod_nln!(
   end
   (nlp.nlcon.nnln == 0) && (Jtv .= 0.0)
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      # Jtv += v[i] * (Aᵢ * x + bᵢ)
-      qcon = nlp.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i])
-      Jtv .+= v[i] .* qcon.b
-    end
+    # Jtv += Jᵀ * v[1:nquad], where J is the Jacobian of the quadratic block.
+    # The product of the evaluator overwrites its output, so it goes through
+    # the workspace `nlp.hv`.
+    MOI.eval_constraint_jacobian_transpose_product(
+      nlp.quadcon.evaluator,
+      nlp.hv,
+      x,
+      view(v, 1:(nlp.quadcon.nquad)),
+    )
+    Jtv .+= nlp.hv
   end
   if nlp.oracles.ncon > 0
     row_offset = nlp.quadcon.nquad + nlp.nlcon.nnln
@@ -597,11 +544,14 @@ function NLPModels.jtprod!(
     )
   end
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[nlp.meta.nlin + i])
-      Jtv .+= v[nlp.meta.nlin + i] .* qcon.b
-    end
+    ind_quad = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
+    MOI.eval_constraint_jacobian_transpose_product(
+      nlp.quadcon.evaluator,
+      nlp.hv,
+      x,
+      view(v, ind_quad),
+    )
+    Jtv .+= nlp.hv
   end
   if nlp.oracles.ncon > 0
     row_offset = nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
@@ -636,12 +586,9 @@ function NLPModels.hess_structure!(
   end
   index = nlp.obj.nnzh
   if nlp.quadcon.nquad > 0
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
-      view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols
-      index += qcon.nnzh
-    end
+    view(rows, (index + 1):(index + nlp.quadcon.nnzh)) .= nlp.quadcon.hess_rows
+    view(cols, (index + 1):(index + nlp.quadcon.nnzh)) .= nlp.quadcon.hess_cols
+    index += nlp.quadcon.nnzh
   end
   if (nlp.obj.type == "NONLINEAR") || (nlp.nlcon.nnln > 0)
     view(rows, (index + 1):(index + nlp.nlcon.nnzh)) .= nlp.nlcon.hess_rows
@@ -687,15 +634,12 @@ function NLPModels.hess_coord!(
     MOI.eval_hessian_lagrangian(nlp.eval, view(vals, ind_nnzh), x, obj_weight, λ_nnln)
   end
 
-  # 3. Quadratic constraint Hessian blocks
+  # 3. Quadratic constraint Hessian blocks. The objective of the block is
+  # zero, so the objective weight is irrelevant.
   if nlp.quadcon.nquad > 0
-    index = nlp.obj.nnzh
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      ind = (index + 1):(index + qcon.nnzh)
-      view(vals, ind) .= y[nlp.meta.nlin + i] .* qcon.A.vals
-      index += qcon.nnzh
-    end
+    ind = (nlp.obj.nnzh + 1):(nlp.obj.nnzh + nlp.quadcon.nnzh)
+    ind_quad = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
+    MOI.eval_hessian_lagrangian(nlp.quadcon.evaluator, view(vals, ind), x, 0.0, view(y, ind_quad))
   end
 
   # 4. Oracle Hessian blocks are appended at the very end
@@ -767,14 +711,11 @@ function NLPModels.jth_hess_coord!(
 
   # Quadratic constraints
   if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad
-    index = nlp.obj.nnzh
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      if j == nlp.meta.nlin + i
-        view(vals, (index + 1):(index + qcon.nnzh)) .= qcon.A.vals
-      end
-      index += qcon.nnzh
-    end
+    μ = nlp.quadcon.y_scratch
+    μ[j - nlp.meta.nlin] = 1.0
+    ind = (nlp.obj.nnzh + 1):(nlp.obj.nnzh + nlp.quadcon.nnzh)
+    MOI.eval_hessian_lagrangian(nlp.quadcon.evaluator, view(vals, ind), x, 0.0, μ)
+    μ[j - nlp.meta.nlin] = 0.0
   end
 
   # Non-oracle nonlinear constraints
@@ -847,10 +788,9 @@ function NLPModels.hprod!(
   end
   if nlp.quadcon.nquad > 0
     (nlp.obj.type == "LINEAR") && (nlp.nlcon.nnln == 0) && (hv .= 0.0)
-    for i = 1:(nlp.quadcon.nquad)
-      qcon = nlp.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nlp.meta.nlin + i])
-    end
+    ind_quad = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
+    MOI.eval_hessian_lagrangian_product(nlp.quadcon.evaluator, nlp.hv, x, v, 0.0, view(y, ind_quad))
+    hv .+= nlp.hv
   end
   if nlp.oracles.ncon > 0
     (nlp.obj.type == "LINEAR") && (nlp.meta.nnln == nlp.oracles.ncon) && (hv .= 0.0)
@@ -920,8 +860,10 @@ function NLPModels.jth_hprod!(
   @rangecheck 1 nlp.meta.ncon j
   hv .= 0.0
   if nlp.meta.nlin + 1 ≤ j ≤ nlp.meta.nlin + nlp.quadcon.nquad
-    qcon = nlp.quadcon.constraints[j - nlp.meta.nlin]
-    coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, 1.0)
+    μ = nlp.quadcon.y_scratch
+    μ[j - nlp.meta.nlin] = 1.0
+    MOI.eval_hessian_lagrangian_product(nlp.quadcon.evaluator, hv, x, v, 0.0, μ)
+    μ[j - nlp.meta.nlin] = 0.0
   elseif nlp.meta.nlin + nlp.quadcon.nquad + 1 ≤
          j ≤
          nlp.meta.nlin + nlp.quadcon.nquad + nlp.nlcon.nnln
@@ -973,9 +915,12 @@ function NLPModels.ghjvprod!(
     error("The function ghjvprod! is not supported by this MathOptNLPModel.")
   increment!(nlp, :neval_hprod)
   ghv .= 0.0
+  μ = nlp.quadcon.y_scratch
   for i = (nlp.meta.nlin + 1):(nlp.meta.nlin + nlp.quadcon.nquad)
-    qcon = nlp.quadcon.constraints[i - nlp.meta.nlin]
-    ghv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, g, v)
+    μ[i - nlp.meta.nlin] = 1.0
+    MOI.eval_hessian_lagrangian_product(nlp.quadcon.evaluator, nlp.hv, x, v, 0.0, μ)
+    μ[i - nlp.meta.nlin] = 0.0
+    ghv[i] = dot(g, nlp.hv)
   end
   for i = (nlp.meta.nlin + nlp.quadcon.nquad + 1):(nlp.meta.ncon)
     jth_hprod!(nlp, x, v, i, nlp.hv)
diff --git a/src/moi_nls_model.jl b/src/moi_nls_model.jl
index eb11cb60..7b3cef0e 100644
--- a/src/moi_nls_model.jl
+++ b/src/moi_nls_model.jl
@@ -13,6 +13,7 @@ mutable struct MathOptNLSModel <: AbstractNLSModel{Float64, Vector{Float64}}
   nlcon::NonLinearStructure
   oracles::Oracles
   λ::Vector{Float64}
+  hv::Vector{Float64}
   counters::NLSCounters
 end
 
@@ -88,6 +89,7 @@ function MathOptNLSModel(cmodel::JuMP.Model, F; hessian::Bool = true, name::Stri
     nlcon,
     oracles,
     λ,
+    zeros(nvar),
     nls_counters,
   )
 end
@@ -270,10 +272,7 @@ function NLPModels.cons_nln!(nls::MathOptNLSModel, x::AbstractVector, c::Abstrac
   offset = 0
   if nls.quadcon.nquad > 0
     offset += nls.quadcon.nquad
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      c[i] = 0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
-    end
+    MOI.eval_constraint(nls.quadcon.evaluator, view(c, 1:(nls.quadcon.nquad)), x)
   end
   if nls.nlcon.nnln > 0
     offset += nls.nlcon.nnln
@@ -303,11 +302,8 @@ function NLPModels.cons!(nls::MathOptNLSModel, x::AbstractVector, c::AbstractVec
   end
   if nls.quadcon.nquad > 0
     offset += nls.quadcon.nquad
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      c[nls.meta.nlin + i] =
-        0.5 * coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, x) + dot(qcon.b, x)
-    end
+    ind_quad = (nls.meta.nlin + 1):(nls.meta.nlin + nls.quadcon.nquad)
+    MOI.eval_constraint(nls.quadcon.evaluator, view(c, ind_quad), x)
   end
   if nls.nlcon.nnln > 0
     offset += nls.nlcon.nnln
@@ -346,14 +342,10 @@ function NLPModels.jac_nln_structure!(
 )
   offset = 0
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
-      qcon = nls.quadcon.constraints[i]
-      ind_quad = (offset + 1):(offset + qcon.nnzg)
-      view(rows, ind_quad) .= i
-      view(cols, ind_quad) .= qcon.g
-      offset += qcon.nnzg
-    end
+    ind_quad = 1:(nls.quadcon.nnzj)
+    view(rows, ind_quad) .= nls.quadcon.jac_rows
+    view(cols, ind_quad) .= nls.quadcon.jac_cols
+    offset += nls.quadcon.nnzj
   end
   if nls.nlcon.nnln > 0
     # non-oracle nonlinear constraints
@@ -390,14 +382,10 @@ function NLPModels.jac_structure!(
     offset += nls.lincon.nnzj
   end
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      # qcon.g is the sparsity pattern of the gradient of the quadratic constraint qcon
-      qcon = nls.quadcon.constraints[i]
-      ind_quad = (offset + 1):(offset + qcon.nnzg)
-      view(rows, ind_quad) .= nls.meta.nlin .+ i
-      view(cols, ind_quad) .= qcon.g
-      offset += qcon.nnzg
-    end
+    ind_quad = (offset + 1):(offset + nls.quadcon.nnzj)
+    view(rows, ind_quad) .= nls.meta.nlin .+ nls.quadcon.jac_rows
+    view(cols, ind_quad) .= nls.quadcon.jac_cols
+    offset += nls.quadcon.nnzj
   end
   if nls.nlcon.nnln > 0
     # non-oracle nonlinear constraints
@@ -434,26 +422,8 @@ function NLPModels.jac_nln_coord!(nls::MathOptNLSModel, x::AbstractVector, vals:
   offset = 0
   if nls.quadcon.nquad > 0
     ind_quad = 1:(nls.quadcon.nnzj)
-    view(vals, ind_quad) .= 0.0
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      for (j, ind) in enumerate(qcon.b.nzind)
-        k = qcon.dg[ind]
-        vals[offset + k] += qcon.b.nzval[j]
-      end
-      for j = 1:(qcon.nnzh)
-        row = qcon.A.rows[j]
-        col = qcon.A.cols[j]
-        val = qcon.A.vals[j]
-        k1 = qcon.dg[row]
-        vals[offset + k1] += val * x[col]
-        if row != col
-          k2 = qcon.dg[col]
-          vals[offset + k2] += val * x[row]
-        end
-      end
-      offset += qcon.nnzg
-    end
+    MOI.eval_constraint_jacobian(nls.quadcon.evaluator, view(vals, ind_quad), x)
+    offset += nls.quadcon.nnzj
   end
   if nls.nlcon.nnln > 0
     ind_nnln = (offset + 1):(offset + nls.nlcon.nnzj)
@@ -484,26 +454,8 @@ function NLPModels.jac_coord!(nls::MathOptNLSModel, x::AbstractVector, vals::Abs
   end
   if nls.quadcon.nquad > 0
     ind_quad = (nls.lincon.nnzj + 1):(nls.lincon.nnzj + nls.quadcon.nnzj)
-    view(vals, ind_quad) .= 0.0
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      for (j, ind) in enumerate(qcon.b.nzind)
-        k = qcon.dg[ind]
-        vals[offset + k] += qcon.b.nzval[j]
-      end
-      for j = 1:(qcon.nnzh)
-        row = qcon.A.rows[j]
-        col = qcon.A.cols[j]
-        val = qcon.A.vals[j]
-        k1 = qcon.dg[row]
-        vals[offset + k1] += val * x[col]
-        if row != col
-          k2 = qcon.dg[col]
-          vals[offset + k2] += val * x[row]
-        end
-      end
-      offset += qcon.nnzg
-    end
+    MOI.eval_constraint_jacobian(nls.quadcon.evaluator, view(vals, ind_quad), x)
+    offset += nls.quadcon.nnzj
   end
   if nls.nlcon.nnln > 0
     ind_nnln = (offset + 1):(offset + nls.nlcon.nnzj)
@@ -550,11 +502,8 @@ function NLPModels.jprod_nln!(
 )
   increment!(nls, :neval_jprod_nln)
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      # Jv[i] += (Aᵢ * x + bᵢ)ᵀ * v
-      qcon = nls.quadcon.constraints[i]
-      Jv[i] = coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
-    end
+    ind_quad = 1:(nls.quadcon.nquad)
+    MOI.eval_constraint_jacobian_product(nls.quadcon.evaluator, view(Jv, ind_quad), x, v)
   end
   if nls.nlcon.nnln > 0
     ind_nnln = (nls.quadcon.nquad + 1):(nls.quadcon.nquad + nls.nlcon.nnln)
@@ -562,7 +511,7 @@ function NLPModels.jprod_nln!(
   end
   if nls.oracles.ncon > 0
     for i =
-      (nls.quadcon.nquad + nls.nlcon.nnln + 1):(nls.quadcon.nquad + nls.nlcon.nnln + nls.oracles.ncon)
+        (nls.quadcon.nquad + nls.nlcon.nnln + 1):(nls.quadcon.nquad + nls.nlcon.nnln + nls.oracles.ncon)
 
       Jv[i] = 0
     end
@@ -606,12 +555,8 @@ function NLPModels.jprod!(
     )
   end
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      # Jv[i] = (Aᵢ * x + bᵢ)ᵀ * v
-      qcon = nls.quadcon.constraints[i]
-      Jv[nls.meta.nlin + i] =
-        coo_sym_dot(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, v) + dot(qcon.b, v)
-    end
+    ind_quad = (nls.meta.nlin + 1):(nls.meta.nlin + nls.quadcon.nquad)
+    MOI.eval_constraint_jacobian_product(nls.quadcon.evaluator, view(Jv, ind_quad), x, v)
   end
   if nls.nlcon.nnln > 0
     ind_nnln =
@@ -672,12 +617,16 @@ function NLPModels.jtprod_nln!(
   end
   (nls.nlcon.nnln == 0) && (Jtv .= 0.0)
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      # Jtv += v[i] * (Aᵢ * x + bᵢ)
-      qcon = nls.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[i])
-      Jtv .+= v[i] .* qcon.b
-    end
+    # Jtv += Jᵀ * v[1:nquad], where J is the Jacobian of the quadratic block.
+    # The product of the evaluator overwrites its output, so it goes through
+    # the workspace `nls.hv`.
+    MOI.eval_constraint_jacobian_transpose_product(
+      nls.quadcon.evaluator,
+      nls.hv,
+      x,
+      view(v, 1:(nls.quadcon.nquad)),
+    )
+    Jtv .+= nls.hv
   end
   if nls.oracles.ncon > 0
     row_offset = nls.quadcon.nquad + nls.nlcon.nnln
@@ -725,11 +674,14 @@ function NLPModels.jtprod!(
     )
   end
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, x, Jtv, v[nls.meta.nlin + i])
-      Jtv .+= v[nls.meta.nlin + i] .* qcon.b
-    end
+    ind_quad = (nls.meta.nlin + 1):(nls.meta.nlin + nls.quadcon.nquad)
+    MOI.eval_constraint_jacobian_transpose_product(
+      nls.quadcon.evaluator,
+      nls.hv,
+      x,
+      view(v, ind_quad),
+    )
+    Jtv .+= nls.hv
   end
   if nls.oracles.ncon > 0
     row_offset = nls.meta.nlin + nls.quadcon.nquad + nls.nlcon.nnln
@@ -765,12 +717,9 @@ function NLPModels.hess_structure!(
   index = nls.lls.nnzh
   if nls.quadcon.nquad > 0
     index = nls.lls.nnzh
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      view(rows, (index + 1):(index + qcon.nnzh)) .= qcon.A.rows
-      view(cols, (index + 1):(index + qcon.nnzh)) .= qcon.A.cols
-      index += qcon.nnzh
-    end
+    view(rows, (index + 1):(index + nls.quadcon.nnzh)) .= nls.quadcon.hess_rows
+    view(cols, (index + 1):(index + nls.quadcon.nnzh)) .= nls.quadcon.hess_cols
+    index += nls.quadcon.nnzh
   end
   if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0)
     view(rows, (index + 1):(index + nls.nlcon.nnzh)) .= nls.nlcon.hess_rows
@@ -803,12 +752,11 @@ function NLPModels.hess_coord!(
     view(vals, 1:(nls.lls.nnzh)) .= obj_weight .* nls.lls.hessian.vals
   end
   if nls.quadcon.nquad > 0
-    index = nls.lls.nnzh
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      view(vals, (index + 1):(index + qcon.nnzh)) .= y[nls.meta.nlin + i] .* qcon.A.vals
-      index += qcon.nnzh
-    end
+    # The objective of the block is zero, so the objective weight is
+    # irrelevant.
+    ind = (nls.lls.nnzh + 1):(nls.lls.nnzh + nls.quadcon.nnzh)
+    ind_quad = (nls.meta.nlin + 1):(nls.meta.nlin + nls.quadcon.nquad)
+    MOI.eval_hessian_lagrangian(nls.quadcon.evaluator, view(vals, ind), x, 0.0, view(y, ind_quad))
   end
   if (nls.nls_meta.nnln > 0) || (nls.nlcon.nnln > 0)
     λ = view(
@@ -904,10 +852,9 @@ function NLPModels.hprod!(
     )
   end
   if nls.quadcon.nquad > 0
-    for i = 1:(nls.quadcon.nquad)
-      qcon = nls.quadcon.constraints[i]
-      coo_sym_add_mul!(qcon.A.rows, qcon.A.cols, qcon.A.vals, v, hv, y[nls.meta.nlin + i])
-    end
+    ind_quad = (nls.meta.nlin + 1):(nls.meta.nlin + nls.quadcon.nquad)
+    MOI.eval_hessian_lagrangian_product(nls.quadcon.evaluator, nls.hv, x, v, 0.0, view(y, ind_quad))
+    hv .+= nls.hv
   end
   if nls.oracles.ncon > 0
     offset_y = nls.meta.nlin + nls.quadcon.nquad + nls.nlcon.nnln
diff --git a/src/utils.jl b/src/utils.jl
index 9013716c..68626549 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -74,21 +74,50 @@ mutable struct LinearConstraints
   nnzj::Int
 end
 
-# xᵀAx + bᵀx
-mutable struct QuadraticConstraint
-  A::COO
-  b::SparseVector{Float64}
-  g::Vector{Int}
-  dg::Dict{Int, Int}
-  nnzg::Int
-  nnzh::Int
-end
+"""
+    QuadraticConstraints
 
+The quadratic constraints, stored in a `MOI.Nonlinear.ModelWithQuad` whose
+rows are the constraints in the order they were parsed, and evaluated through
+the corresponding `MOI.Nonlinear.EvaluatorWithQuad`. The Jacobian and Hessian
+structures of the block are precomputed. `y_scratch` holds basis multiplier
+vectors for the per-constraint methods (`jth_hess_coord!`, `jth_hprod!` and
+`ghjvprod!`).
+"""
 mutable struct QuadraticConstraints
   nquad::Int
-  constraints::Vector{QuadraticConstraint}
+  evaluator::MOI.Nonlinear.EvaluatorWithQuad
+  jac_rows::Vector{Int}
+  jac_cols::Vector{Int}
   nnzj::Int
+  hess_rows::Vector{Int}
+  hess_cols::Vector{Int}
   nnzh::Int
+  y_scratch::Vector{Float64}
+end
+
+function QuadraticConstraints(quad_model::MOI.Nonlinear.ModelWithQuad)
+  nquad = length(quad_model)
+  evaluator = MOI.Nonlinear.Evaluator(quad_model, MOI.Nonlinear.SparseReverseMode())
+  MOI.initialize(evaluator, [:Grad, :Jac, :JacVec, :Hess, :HessVec])
+  jac_structure = MOI.jacobian_structure(evaluator)
+  jac_rows = [r for (r, _) in jac_structure]
+  jac_cols = [c for (_, c) in jac_structure]
+  hess_structure = MOI.hessian_lagrangian_structure(evaluator)
+  # NLPModels expects the lower triangle
+  hess_rows = [max(r, c) for (r, c) in hess_structure]
+  hess_cols = [min(r, c) for (r, c) in hess_structure]
+  return QuadraticConstraints(
+    nquad,
+    evaluator,
+    jac_rows,
+    jac_cols,
+    length(jac_structure),
+    hess_rows,
+    hess_cols,
+    length(hess_structure),
+    zeros(nquad),
+  )
 end
 
 """
@@ -319,127 +348,26 @@ end
 Parse a `ScalarQuadraticFunction` fun with its associated set.
 `qcons`, `quad_lcon`, `quad_ucon` are updated.
 """
-function parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
-  _index(v::MOI.VariableIndex) = index_map[v].value
-
-  b = spzeros(Float64, nvar)
-  rows = Int[]
-  cols = Int[]
-  vals = Float64[]
-
-  # Parse a ScalarAffineTerm{Float64}(coefficient, variable_index)
-  for term in fun.affine_terms
-    b[_index(term.variable)] = term.coefficient
-  end
-
-  # Parse a ScalarQuadraticTerm{Float64}(coefficient, variable_index_1, variable_index_2)
-  for term in fun.quadratic_terms
-    i = _index(term.variable_1)
-    j = _index(term.variable_2)
-    if i ≥ j
-      push!(rows, i)
-      push!(cols, j)
-    else
-      push!(rows, j)
-      push!(cols, i)
-    end
-    push!(vals, term.coefficient)
-  end
-
-  if typeof(set) in (MOI.Interval{Float64}, MOI.GreaterThan{Float64})
-    push!(quad_lcon, -fun.constant + set.lower)
-  elseif typeof(set) == MOI.EqualTo{Float64}
-    push!(quad_lcon, -fun.constant + set.value)
-  else
-    push!(quad_lcon, -Inf)
-  end
-
-  if typeof(set) in (MOI.Interval{Float64}, MOI.LessThan{Float64})
-    push!(quad_ucon, -fun.constant + set.upper)
-  elseif typeof(set) == MOI.EqualTo{Float64}
-    push!(quad_ucon, -fun.constant + set.value)
-  else
-    push!(quad_ucon, Inf)
-  end
-
-  A = COO(rows, cols, vals)
-  g = unique(vcat(rows, cols, b.nzind))  # sparsity pattern of Ax + b
-  nnzg = length(g)
-  # dg is a dictionary where:
-  # - The key `r` specifies a row index in the vector Ax + b.
-  # - The value `dg[r]` is a position in the vector (of length nnzg)
-  # where the non-zero entries of the Jacobian for row `r` are stored.
-  dg = Dict{Int, Int}(g[p] => p for p = 1:nnzg)
-  nnzh = length(vals)
-  qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh)
-  push!(qcons, qcon)
+function parser_SQF(fun, set, quad_model, index_map)
+  f = MOI.Utilities.map_indices(index_map, fun)
+  # The constant is moved into the set, as MOI requires for
+  # scalar-function-in-set constraints.
+  g = SQF(f.quadratic_terms, f.affine_terms, 0.0)
+  MOI.add_constraint(quad_model, g, MOI.Utilities.shift_constant(set, -f.constant))
+  return
 end
 
-"""
-    parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
-
-Parse a `VectorQuadraticFunction` fun with its associated set.
-`qcons`, `quad_lcon`, `quad_ucon` are updated.
-"""
-function parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
-  _index(v::MOI.VariableIndex) = index_map[v].value
-
-  ncon = length(fun.constants)
-  for k = 1:ncon
-    b = spzeros(Float64, nvar)
-    rows = Int[]
-    cols = Int[]
-    vals = Float64[]
-
-    # Parse a VectorAffineTerm{Float64}(output_index, scalar_term)
-    for affine_term in fun.affine_terms
-      if affine_term.output_index == k
-        b[_index(affine_term.scalar_term.variable)] = affine_term.scalar_term.coefficient
-      end
-    end
-
-    # Parse a VectorQuadraticTerm{Float64}(output_index, scalar_term)
-    for quadratic_term in fun.quadratic_terms
-      if quadratic_term.output_index == k
-        i = _index(quadratic_term.scalar_term.variable_1)
-        j = _index(quadratic_term.scalar_term.variable_2)
-        if i ≥ j
-          push!(rows, i)
-          push!(cols, j)
-        else
-          push!(rows, j)
-          push!(cols, i)
-        end
-        push!(vals, quadratic_term.scalar_term.coefficient)
-      end
-    end
+_scalar_set(::MOI.Nonnegatives) = MOI.GreaterThan(0.0)
+_scalar_set(::MOI.Nonpositives) = MOI.LessThan(0.0)
+_scalar_set(::MOI.Zeros) = MOI.EqualTo(0.0)
 
-    constant = fun.constants[k]
-
-    if typeof(set) in (MOI.Nonnegatives, MOI.Zeros)
-      append!(quad_lcon, constant)
-    else
-      append!(quad_lcon, -Inf)
-    end
-
-    if typeof(set) in (MOI.Nonpositives, MOI.Zeros)
-      append!(quad_ucon, -constant)
-    else
-      append!(quad_ucon, Inf)
-    end
-
-    A = COO(rows, cols, vals)
-    g = unique(vcat(rows, cols, b.nzind))  # sparsity pattern of Ax + b
-    nnzg = length(g)
-    # dg is a dictionary where:
-    # - The key `r` specifies a row index in the vector Ax + b.
-    # - The value `dg[r]` is a position in the vector (of length nnzg)
-    # where the non-zero entries of the Jacobian for row `r` are stored.
-    dg = Dict{Int, Int}(g[p] => p for p = 1:nnzg)
-    nnzh = length(vals)
-    qcon = QuadraticConstraint(A, b, g, dg, nnzg, nnzh)
-    push!(qcons, qcon)
+function parser_VQF(fun, set, quad_model, index_map)
+  f = MOI.Utilities.map_indices(index_map, fun)
+  for fi in MOI.Utilities.scalarize(f)
+    g = SQF(fi.quadratic_terms, fi.affine_terms, 0.0)
+    MOI.add_constraint(quad_model, g, MOI.Utilities.shift_constant(_scalar_set(set), -fi.constant))
   end
+  return
 end
 
 """
@@ -459,9 +387,10 @@ function parser_MOI(moimodel, index_map, nvar)
 
   # Variables associated to quadratic constraints
   nquad = 0
-  qcons = QuadraticConstraint[]
-  quad_lcon = Float64[]
-  quad_ucon = Float64[]
+  quad_model = MOI.Nonlinear.ModelWithQuad(MOI.Nonlinear.Model())
+  for _ = 1:nvar
+    MOI.add_variable(quad_model)
+  end
 
   contypes = MOI.get(moimodel, MOI.ListOfConstraintTypesPresent())
   for (F, S) in contypes
@@ -494,11 +423,11 @@ function parser_MOI(moimodel, index_map, nvar)
         nlin += set.dimension
       end
       if typeof(fun) <: SQF
-        parser_SQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+        parser_SQF(fun, set, quad_model, index_map)
         nquad += 1
       end
       if typeof(fun) <: VQF
-        parser_VQF(fun, set, nvar, qcons, quad_lcon, quad_ucon, index_map)
+        parser_VQF(fun, set, quad_model, index_map)
         nquad += set.dimension
       end
     end
@@ -506,13 +435,10 @@ function parser_MOI(moimodel, index_map, nvar)
   coo = COO(linrows, lincols, linvals)
   lin_nnzj = length(linvals)
   lincon = LinearConstraints(coo, lin_nnzj)
-  quad_nnzj = 0
-  quad_nnzh = 0
-  for i = 1:nquad
-    quad_nnzj += qcons[i].nnzg
-    quad_nnzh += qcons[i].nnzh
-  end
-  quadcon = QuadraticConstraints(nquad, qcons, quad_nnzj, quad_nnzh)
+  quadcon = QuadraticConstraints(quad_model)
+  quad_bounds = MOI.NLPBlockData(quadcon.evaluator).constraint_bounds
+  quad_lcon = [b.lower for b in quad_bounds]
+  quad_ucon = [b.upper for b in quad_bounds]
 
   return nlin, lincon, lin_lcon, lin_ucon, quadcon, quad_lcon, quad_ucon
 end