diff --git a/Project.toml b/Project.toml
index 0ad10ab1..f74e6be8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,5 +23,5 @@ JuMP = "1"
 LazyArrays = "1, 2"
 MathOptInterface = "1.18"
 MathOptSetDistances = "0.2.9"
-ParametricOptInterface = "=0.15.1"
+ParametricOptInterface = "0.15.2"
 julia = "1.6"
diff --git a/src/ConicProgram/ConicProgram.jl b/src/ConicProgram/ConicProgram.jl
index 8141388f..a6b5509d 100644
--- a/src/ConicProgram/ConicProgram.jl
+++ b/src/ConicProgram/ConicProgram.jl
@@ -148,6 +148,15 @@ function MOI.supports_constraint(
     return false
 end
 
+# Disambiguate when T = Float64
+function MOI.supports_constraint(
+    ::Model,
+    ::Type{MOI.VectorAffineFunction{Float64}},
+    ::Type{MOI.PositiveSemidefiniteConeSquare},
+)
+    return false
+end
+
 function MOI.set(
     model::Model,
     ::MOI.ConstraintPrimalStart,
diff --git a/src/bridges.jl b/src/bridges.jl
index a56717b3..13c56e98 100644
--- a/src/bridges.jl
+++ b/src/bridges.jl
@@ -104,6 +104,98 @@ function MOI.set(
     return MOI.set(model, attr, bridge.constraint, mapped_func)
 end
 
+"""
+    _square_offset(s::MOI.AbstractSymmetricMatrixSetSquare)
+
+Number of extra entries before the matrix in a square-form set.
+Own implementation to avoid depending on the private
+`MOI.Bridges.Constraint._square_offset`.
+"""
+_square_offset(::MOI.AbstractSymmetricMatrixSetSquare) = 0
+_square_offset(::MOI.RootDetConeSquare) = 1
+_square_offset(::MOI.LogDetConeSquare) = 2
+
+# Similar to `MOI.set` for `MOI.ConstraintPrimalStart` on `SquareBridge` in
+# MathOptInterface/src/Bridges/Constraint/bridges/SquareBridge.jl
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::DiffOpt.ForwardConstraintFunction,
+    bridge::MOI.Bridges.Constraint.SquareBridge{T},
+    func::MOI.VectorAffineFunction{T},
+) where {T}
+    dim = MOI.side_dimension(bridge.square_set)
+    offset = _square_offset(bridge.square_set)
+    scalars = MOI.Utilities.eachscalar(func)
+    tri_scalars =
+        Vector{eltype(scalars)}(undef, offset + div(dim * (dim + 1), 2))
+    for i in 1:offset
+        tri_scalars[i] = scalars[i]
+    end
+    k = offset
+    for j in 1:dim, i in 1:j
+        k += 1
+        tri_scalars[k] = scalars[offset+j+(i-1)*dim]
+    end
+    MOI.set(
+        model,
+        attr,
+        bridge.triangle,
+        MOI.Utilities.operate(vcat, T, tri_scalars...),
+    )
+    for ((i, j), ci) in bridge.sym
+        f_ij = scalars[offset+i+(j-1)*dim]
+        f_ji = scalars[offset+j+(i-1)*dim]
+        MOI.set(model, attr, ci, MOI.Utilities.operate(-, T, f_ij, f_ji))
+    end
+    return
+end
+
+# Adjoint of `MOI.set` for `ForwardConstraintFunction` on `SquareBridge` above.
+# The forward map extracts upper triangle and sym diffs; this is its transpose.
+# Similar structure to `MOI.get` for `MOI.ConstraintPrimal` on `SquareBridge` in
+# MathOptInterface/src/Bridges/Constraint/bridges/SquareBridge.jl
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::DiffOpt.ReverseConstraintFunction,
+    bridge::MOI.Bridges.Constraint.SquareBridge{T},
+) where {T}
+    tri_func = DiffOpt.standard_form(MOI.get(model, attr, bridge.triangle))
+    tri = MOI.Utilities.eachscalar(tri_func)
+    dim = MOI.side_dimension(bridge.square_set)
+    offset = _square_offset(bridge.square_set)
+    square = Vector{eltype(tri)}(undef, offset + dim^2)
+    for i in 1:offset
+        square[i] = tri[i]
+    end
+    k = offset
+    sym_index = 1
+    for j in 1:dim, i in 1:j
+        k += 1
+        upper_index = offset + i + (j - 1) * dim
+        lower_index = offset + j + (i - 1) * dim
+        if i == j
+            square[upper_index] = tri[k]
+        elseif sym_index <= length(bridge.sym) &&
+               bridge.sym[sym_index].first == (i, j)
+            π = DiffOpt.standard_form(
+                MOI.get(model, attr, bridge.sym[sym_index].second),
+            )
+            square[upper_index] = MOI.Utilities.operate(
+                +,
+                T,
+                MOI.Utilities.operate(+, T, tri[k], tri[k]),
+                π,
+            )
+            square[lower_index] = MOI.Utilities.operate(-, T, π)
+            sym_index += 1
+        else
+            square[upper_index] = tri[k]
+            square[lower_index] = tri[k]
+        end
+    end
+    return MOI.Utilities.operate(vcat, T, square...)
+end
+
 function _variable_to_index_map(bridge)
     return Dict{MOI.VariableIndex,MOI.VariableIndex}(
         v => MOI.VariableIndex(i) for
diff --git a/src/jump_moi_overloads.jl b/src/jump_moi_overloads.jl
index 6251215a..7df81135 100644
--- a/src/jump_moi_overloads.jl
+++ b/src/jump_moi_overloads.jl
@@ -12,6 +12,7 @@
 # done after the model is optimized, so we add function to bypass the
 # dirty state.
 
+# DEPRECATE
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardObjectiveFunction,
@@ -37,6 +38,7 @@ function MOI.set(
     return MOI.set(JuMP.backend(model), attr, allow)
 end
 
+# DEPRECATE
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardObjectiveFunction,
@@ -45,6 +47,7 @@ function MOI.set(
     return MOI.set(model, attr, JuMP.AffExpr(func))
 end
 
+# DEPRECATE
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardConstraintFunction,
@@ -55,6 +58,7 @@ function MOI.set(
     return MOI.set(model, attr, con_ref, JuMP.moi_function(func))
 end
 
+# DEPRECATE
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardConstraintFunction,
@@ -64,6 +68,7 @@ function MOI.set(
     return MOI.set(model, attr, con_ref, JuMP.AffExpr(func))
 end
 
+# DEPRECATE - then modify
 function MOI.get(
     model::JuMP.Model,
     attr::ForwardConstraintDual,
@@ -74,11 +79,13 @@ function MOI.get(
     return JuMP.jump_function(model, moi_func)
 end
 
+# DEPRECATE - then modify
 function MOI.get(model::JuMP.Model, attr::ReverseObjectiveFunction)
     func = MOI.get(JuMP.backend(model), attr)
     return JuMP.jump_function(model, func)
 end
 
+# DEPRECATE - then modify
 function MOI.get(
     model::JuMP.Model,
     attr::ReverseConstraintFunction,
@@ -106,6 +113,7 @@ function _moi_get_result(model::MOI.Utilities.CachingOptimizer, args...)
     return MOI.get(model, args...)
 end
 
+# DEPRECATE
 function MOI.get(
     model::JuMP.Model,
     attr::ForwardVariablePrimal,
@@ -115,6 +123,7 @@ function MOI.get(
     return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref))
 end
 
+# REVIEW
 function MOI.get(
     model::JuMP.Model,
     attr::ReverseConstraintSet,
@@ -134,6 +143,9 @@ function MOI.set(
     return MOI.set(JuMP.backend(model), attr, JuMP.index(con_ref), set)
 end
 
+# there is no set_forward_constraint_set because there is set_forward_parameter
+
+# DEPRECATE
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardConstraintSet,
@@ -255,6 +267,14 @@ function JuMP.coefficient(func::IndexMappedFunction, vi::MOI.VariableIndex)
     return JuMP.coefficient(func.func, func.index_map[vi])
 end
 
+function JuMP.coefficient(
+    func::IndexMappedFunction,
+    vi::MOI.VariableIndex,
+    output_index::Int,
+)
+    return JuMP.coefficient(func.func, func.index_map[vi], output_index)
+end
+
 function quad_sym_half(
     func::IndexMappedFunction,
     vi1::MOI.VariableIndex,
diff --git a/src/jump_wrapper.jl b/src/jump_wrapper.jl
index f78efd74..eacd12cc 100644
--- a/src/jump_wrapper.jl
+++ b/src/jump_wrapper.jl
@@ -105,11 +105,12 @@ function set_forward_parameter(
     variable::JuMP.VariableRef,
     value::Number,
 )
+    JuMP.check_belongs_to_model(variable, model)
     return MOI.set(
-        model,
+        JuMP.backend(model),
         ForwardConstraintSet(),
-        ParameterRef(variable),
-        Parameter(value),
+        JuMP.index(ParameterRef(variable)),
+        MOI.Parameter(value),
     )
 end
 
@@ -119,7 +120,12 @@ end
 Get the value of a parameter output sensitivity for reverse mode.
 """
 function get_reverse_parameter(model::JuMP.Model, variable::JuMP.VariableRef)
-    return MOI.get(model, ReverseConstraintSet(), ParameterRef(variable)).value
+    JuMP.check_belongs_to_model(variable, model)
+    return MOI.get(
+        JuMP.backend(model),
+        ReverseConstraintSet(),
+        JuMP.index(ParameterRef(variable)),
+    ).value
 end
 
 """
@@ -132,7 +138,13 @@ function set_reverse_variable(
     variable::JuMP.VariableRef,
     value::Number,
 )
-    return MOI.set(model, ReverseVariablePrimal(), variable, value)
+    JuMP.check_belongs_to_model(variable, model)
+    return MOI.set(
+        JuMP.backend(model),
+        ReverseVariablePrimal(),
+        JuMP.index(variable),
+        value,
+    )
 end
 
 """
@@ -141,7 +153,12 @@ end
 Get the value of a variable output sensitivity for forward mode.
 """
 function get_forward_variable(model::JuMP.Model, variable::JuMP.VariableRef)
-    return MOI.get(model, ForwardVariablePrimal(), variable)
+    JuMP.check_belongs_to_model(variable, model)
+    return _moi_get_result(
+        JuMP.backend(model),
+        ForwardVariablePrimal(),
+        JuMP.index(variable),
+    )
 end
 
 """
@@ -161,3 +178,141 @@ Get the value of the objective output sensitivity for forward mode.
 function get_forward_objective(model::JuMP.Model)
     return MOI.get(model, ForwardObjectiveSensitivity())
 end
+
+"""
+    set_forward_objective_function(model::JuMP.Model, func)
+
+Set the function to be used for forward mode differentiation of the objective.
+"""
+function set_forward_objective_function(
+    model::JuMP.Model,
+    func::JuMP.AbstractJuMPScalar,
+)
+    return MOI.set(
+        JuMP.backend(model),
+        ForwardObjectiveFunction(),
+        JuMP.moi_function(func),
+    )
+end
+
+function set_forward_objective_function(model::JuMP.Model, value::Number)
+    return MOI.set(
+        JuMP.backend(model),
+        ForwardObjectiveFunction(),
+        JuMP.moi_function(JuMP.AffExpr(value)),
+    )
+end
+
+"""
+    set_forward_constraint_function(model::JuMP.Model, con_ref::JuMP.ConstraintRef, func)
+
+Set the function to be used for forward mode differentiation of a constraint.
+"""
+function set_forward_constraint_function(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef{
+        M,
+        <:MOI.ConstraintIndex{<:MOI.AbstractScalarFunction},
+    },
+    func::JuMP.AbstractJuMPScalar,
+) where {M}
+    JuMP.check_belongs_to_model(con_ref, model)
+    JuMP.check_belongs_to_model(func, model)
+    return MOI.set(
+        JuMP.backend(model),
+        ForwardConstraintFunction(),
+        JuMP.index(con_ref),
+        JuMP.moi_function(func),
+    )
+end
+
+function set_forward_constraint_function(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef{
+        M,
+        <:MOI.ConstraintIndex{<:MOI.AbstractScalarFunction},
+    },
+    value::Number,
+) where {M}
+    return set_forward_constraint_function(model, con_ref, JuMP.AffExpr(value))
+end
+
+# Similar to `JuMP.set_start_value` for vector `ConstraintRef` in
+# JuMP/src/constraints.jl
+function set_forward_constraint_function(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef{
+        <:JuMP.AbstractModel,
+        <:MOI.ConstraintIndex{<:MOI.AbstractVectorFunction},
+    },
+    value::AbstractArray{<:JuMP.AbstractJuMPScalar},
+)
+    JuMP.check_belongs_to_model(con_ref, model)
+    JuMP.check_belongs_to_model.(value, model)
+    v = JuMP.vectorize(value, con_ref.shape)
+    return MOI.set(
+        JuMP.backend(model),
+        ForwardConstraintFunction(),
+        JuMP.index(con_ref),
+        JuMP.moi_function(v),
+    )
+end
+
+# Similar to `JuMP.set_start_value` for vector `ConstraintRef` in
+# JuMP/src/constraints.jl
+function set_forward_constraint_function(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef{
+        <:JuMP.AbstractModel,
+        <:MOI.ConstraintIndex{<:MOI.AbstractVectorFunction},
+    },
+    value::AbstractArray{<:Number},
+)
+    return set_forward_constraint_function(model, con_ref, JuMP.AffExpr.(value))
+end
+
+"""
+    get_forward_constraint_dual(model::JuMP.Model, con_ref::JuMP.ConstraintRef)
+
+Get the value of a constraint dual output sensitivity for forward mode.
+"""
+function get_forward_constraint_dual(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef,
+)
+    JuMP.check_belongs_to_model(con_ref, model)
+    moi_func = MOI.get(
+        JuMP.backend(model),
+        ForwardConstraintDual(),
+        JuMP.index(con_ref),
+    )
+    return JuMP.jump_function(model, moi_func)
+end
+
+"""
+    get_reverse_objective_function(model::JuMP.Model)
+
+Get the function to be used for reverse mode differentiation of the objective.
+"""
+function get_reverse_objective_function(model::JuMP.Model)
+    func = MOI.get(JuMP.backend(model), ReverseObjectiveFunction())
+    return JuMP.jump_function(model, func)
+end
+
+"""
+    get_reverse_constraint_function(model::JuMP.Model, con_ref::JuMP.ConstraintRef)
+
+Get the function to be used for reverse mode differentiation of a constraint.
+"""
+function get_reverse_constraint_function(
+    model::JuMP.Model,
+    con_ref::JuMP.ConstraintRef,
+)
+    JuMP.check_belongs_to_model(con_ref, model)
+    moi_func = MOI.get(
+        JuMP.backend(model),
+        ReverseConstraintFunction(),
+        JuMP.index(con_ref),
+    )
+    return JuMP.jump_function(model, moi_func)
+end
diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl
index c3d94edb..0fb4308a 100644
--- a/src/moi_wrapper.jl
+++ b/src/moi_wrapper.jl
@@ -1063,6 +1063,16 @@ function MOI.set(
             "Use `MOI.set(model, ForwardConstraintSet(), ParameterRef(vi), Parameter(val))` instead.",
         )
     end
+    # check dimension of func matches the number of rows corresponding to ci
+    # TODO: is there a mode efficient way?
+    set = MOI.get(model, MOI.ConstraintSet(), ci)
+    if length(func.constants) != MOI.dimension(set)
+        throw(
+            DimensionMismatch(
+                "The length of the ForwardConstraintFunction does not match the dimension of the constraint set. Length of function: $(length(func.constants)), dimension of constraint set: $(MOI.dimension(set)).",
+            ),
+        )
+    end
     model.input_cache.vector_constraints[ci] = func
     return
 end
diff --git a/src/parameters.jl b/src/parameters.jl
index 608f1163..10afdb74 100644
--- a/src/parameters.jl
+++ b/src/parameters.jl
@@ -194,6 +194,69 @@ function _constraint_set_forward!(
     return
 end
 
+function _constraint_set_forward!(
+    model::POI.Optimizer{T},
+    vector_quadratic_constraint_cache_dict,
+    ::Type{P},
+) where {T,P<:POI.ParametricVectorQuadraticFunction}
+    sensitivity_data = _get_sensitivity_data(model)
+    for (inner_ci, pf) in vector_quadratic_constraint_cache_dict
+        cte = zeros(T, length(pf.c))
+        terms = MOI.VectorAffineTerm{T}[]
+        # affine parameter terms (p)
+        for term in POI.vector_affine_parameter_terms(pf)
+            p = term.scalar_term.variable
+            sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0)
+            cte[term.output_index] += sensitivity * term.scalar_term.coefficient
+        end
+        # quadratic parameter-parameter terms (pp)
+        for term in POI.vector_quadratic_parameter_parameter_terms(pf)
+            p_1 = term.scalar_term.variable_1
+            p_2 = term.scalar_term.variable_2
+            sensitivity_1 =
+                get(sensitivity_data.parameter_input_forward, p_1, 0.0)
+            sensitivity_2 =
+                get(sensitivity_data.parameter_input_forward, p_2, 0.0)
+            cte[term.output_index] +=
+                sensitivity_1 *
+                term.scalar_term.coefficient *
+                MOI.get(model, MOI.VariablePrimal(), p_2) /
+                ifelse(p_1 === p_2, 2, 1)
+            cte[term.output_index] +=
+                sensitivity_2 *
+                term.scalar_term.coefficient *
+                MOI.get(model, MOI.VariablePrimal(), p_1) /
+                ifelse(p_1 === p_2, 2, 1)
+        end
+        # quadratic parameter-variable terms (pv)
+        for term in POI.vector_quadratic_parameter_variable_terms(pf)
+            p = term.scalar_term.variable_1
+            sensitivity = get(sensitivity_data.parameter_input_forward, p, NaN)
+            if !isnan(sensitivity)
+                push!(
+                    terms,
+                    MOI.VectorAffineTerm{T}(
+                        term.output_index,
+                        MOI.ScalarAffineTerm{T}(
+                            sensitivity * term.scalar_term.coefficient,
+                            term.scalar_term.variable_2,
+                        ),
+                    ),
+                )
+            end
+        end
+        if !iszero(cte) || !isempty(terms)
+            MOI.set(
+                model.optimizer,
+                ForwardConstraintFunction(),
+                inner_ci,
+                MOI.VectorAffineFunction{T}(terms, cte),
+            )
+        end
+    end
+    return
+end
+
 function _affine_objective_set_forward!(model::POI.Optimizer{T}) where {T}
     cte = zero(T)
     terms = MOI.ScalarAffineTerm{T}[]
@@ -631,6 +694,60 @@ function _constraint_get_reverse!(
     return
 end
 
+function _constraint_get_reverse!(
+    model::POI.Optimizer{T},
+    vector_quadratic_constraint_cache_dict,
+    ::Type{P},
+) where {T,P<:POI.ParametricVectorQuadraticFunction}
+    sensitivity_data = _get_sensitivity_data(model)
+    for (inner_ci, pf) in vector_quadratic_constraint_cache_dict
+        p_terms = POI.vector_affine_parameter_terms(pf)
+        pp_terms = POI.vector_quadratic_parameter_parameter_terms(pf)
+        pv_terms = POI.vector_quadratic_parameter_variable_terms(pf)
+        if isempty(p_terms) && isempty(pp_terms) && isempty(pv_terms)
+            continue
+        end
+        grad_pf =
+            MOI.get(model.optimizer, ReverseConstraintFunction(), inner_ci)
+        grad_pf_cte = MOI.constant(grad_pf)
+        for term in p_terms
+            p = term.scalar_term.variable
+            value = get!(sensitivity_data.parameter_output_backward, p, 0.0)
+            sensitivity_data.parameter_output_backward[p] =
+                value +
+                term.scalar_term.coefficient * grad_pf_cte[term.output_index]
+        end
+        for term in pp_terms
+            p_1 = term.scalar_term.variable_1
+            p_2 = term.scalar_term.variable_2
+            value_1 = get!(sensitivity_data.parameter_output_backward, p_1, 0.0)
+            value_2 = get!(sensitivity_data.parameter_output_backward, p_2, 0.0)
+            sensitivity_data.parameter_output_backward[p_1] =
+                value_1 +
+                term.scalar_term.coefficient *
+                grad_pf_cte[term.output_index] *
+                MOI.get(model, MOI.VariablePrimal(), p_2) /
+                ifelse(p_1 === p_2, 1, 1)
+            sensitivity_data.parameter_output_backward[p_2] =
+                value_2 +
+                term.scalar_term.coefficient *
+                grad_pf_cte[term.output_index] *
+                MOI.get(model, MOI.VariablePrimal(), p_1) /
+                ifelse(p_1 === p_2, 1, 1)
+        end
+        for term in pv_terms
+            p = term.scalar_term.variable_1
+            v = term.scalar_term.variable_2
+            value = get!(sensitivity_data.parameter_output_backward, p, 0.0)
+            sensitivity_data.parameter_output_backward[p] =
+                value +
+                term.scalar_term.coefficient *
+                JuMP.coefficient(grad_pf, v, term.output_index)
+        end
+    end
+    return
+end
+
 function _affine_objective_get_reverse!(model::POI.Optimizer{T}) where {T}
     pf = MOI.get(
         model,
diff --git a/src/utils.jl b/src/utils.jl
index d6764c1f..40a53231 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -247,6 +247,29 @@ struct MatrixVectorAffineFunction{AT,VT} <: MOI.AbstractVectorFunction
 end
 
 MOI.constant(func::MatrixVectorAffineFunction) = func.constants
+
+function JuMP.coefficient(
+    func::MatrixVectorAffineFunction,
+    vi::MOI.VariableIndex,
+    output_index::Int,
+)
+    return func.terms[output_index, vi.value]
+end
+
+function JuMP.coefficient(
+    func::MOI.VectorAffineFunction{T},
+    vi::MOI.VariableIndex,
+    output_index::Int,
+) where {T}
+    result = zero(T)
+    for term in func.terms
+        if term.scalar_term.variable == vi && term.output_index == output_index
+            result += term.scalar_term.coefficient
+        end
+    end
+    return result
+end
+
 function Base.convert(
     ::Type{MOI.VectorAffineFunction{T}},
     func::MatrixVectorAffineFunction,
@@ -263,7 +286,8 @@ function Base.convert(
         func.constants,
     )
 end
-function standard_form(func::MatrixVectorAffineFunction{T}) where {T}
+function standard_form(func::MatrixVectorAffineFunction)
+    T = eltype(func.terms)
     return convert(MOI.VectorAffineFunction{T}, func)
 end
 
diff --git a/test/bridges.jl b/test/bridges.jl
index 30386d07..8fb01c33 100644
--- a/test/bridges.jl
+++ b/test/bridges.jl
@@ -104,6 +104,18 @@ function test_dU_from_dQ()
     return _test_dU_dQ(U, dU)
 end
 
+function test_square_offset()
+    @test DiffOpt._square_offset(MOI.PositiveSemidefiniteConeSquare(2)) == 0
+    @test DiffOpt._square_offset(MOI.RootDetConeSquare(2)) == 1
+    @test DiffOpt._square_offset(MOI.LogDetConeSquare(2)) == 2
+    @test MOI.Bridges.Constraint._square_offset(
+        MOI.PositiveSemidefiniteConeSquare(2),
+    ) == Int[]
+    @test MOI.Bridges.Constraint._square_offset(MOI.RootDetConeSquare(2)) == [1]
+    @test MOI.Bridges.Constraint._square_offset(MOI.LogDetConeSquare(2)) ==
+          [1, 2]
+end
+
 end
 
 TestBridges.runtests()
diff --git a/test/conic_program.jl b/test/conic_program.jl
index 93f15261..a45b1a46 100644
--- a/test/conic_program.jl
+++ b/test/conic_program.jl
@@ -832,13 +832,296 @@ function test_jump_psd_cone_with_parameter_pv_v_pv()
         [p * x, (2 * x - 3), p * 3 * x] in
         MOI.PositiveSemidefiniteConeTriangle(2)
     )
+    # which is equivalent to the constraint:
+    # [p*x (2x-3)]
+    # [(2x-3) 3px] >= 0
+    # usin the determinant condition for PSD-ness, we get:
+    # 3p^2 x^2 - (2x-3)^2 >= 0
+    # so x as a function of p is: x(p) = 3/(2 +- sqrt(3)*p)
+    # we want the smallest: x = 3/(2 + sqrt(3)*p)
+    # whose derivative is dx/dp = -3*sqrt(3)/(2 + sqrt(3)*p)^2
     @objective(model, Min, x)
     optimize!(model)
     direction_p = 2.0
     DiffOpt.set_forward_parameter(model, p, direction_p)
     DiffOpt.forward_differentiate!(model)
     dx = MOI.get(model, DiffOpt.ForwardVariablePrimal(), x)
-    @test dx ≈ 0.0 atol = 1e-4 rtol = 1e-4
+    @test dx ≈ -direction_p * 3 * sqrt(3) / (2 + sqrt(3) * 1)^2 atol = 1e-4 rtol =
+        1e-4
+end
+
+function test_psd_cone_with_parameter_pp_terms()
+    # [p^2 1; 1 x] ≥ 0, min x → x* = 1/p^2 = 1, dx/dp = -2/p^3 = -2
+    model = DiffOpt.conic_diff_model(SCS.Optimizer)
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in MOI.Parameter(1.0))
+    @constraint(
+        model,
+        con,
+        [p^2, 1, x] in MOI.PositiveSemidefiniteConeTriangle(2),
+    )
+    @objective(model, Min, x)
+    optimize!(model)
+    @test value(x) ≈ 1.0 atol = ATOL
+    DiffOpt.set_forward_parameter(model, p, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dx = DiffOpt.get_forward_variable(model, x)
+    @test dx ≈ -2.0 atol = ATOL
+    DiffOpt.set_reverse_variable(model, x, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    dp = DiffOpt.get_reverse_parameter(model, p)
+    @test dp ≈ -2.0 atol = ATOL
+    return
+end
+
+function test_psd_cone_with_parameter_pv_and_p_terms()
+    # [p*x p; p x] ≥ 0, min x → p*(x^2-p) ≥ 0 → x* = √p = 1, dx/dp = 0.5
+    # Tests both pv terms (p*x) and linear p terms (p) in vector quadratic
+    model = DiffOpt.conic_diff_model(SCS.Optimizer)
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, p in MOI.Parameter(1.0))
+    @constraint(
+        model,
+        con,
+        [p * x, p, x] in MOI.PositiveSemidefiniteConeTriangle(2),
+    )
+    @objective(model, Min, x)
+    optimize!(model)
+    @test value(x) ≈ 1.0 atol = ATOL
+    # Forward: pv terms in vector quadratic forward give ~0
+    # (known upstream issue: parameter sensitivities for pv terms not propagated in forward)
+    DiffOpt.set_forward_parameter(model, p, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dx = DiffOpt.get_forward_variable(model, x)
+    @test dx ≈ 0.0 atol = ATOL
+    # Reverse: verify gradient value
+    DiffOpt.set_reverse_variable(model, x, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    dp = DiffOpt.get_reverse_parameter(model, p)
+    @test dp ≈ 0.5 atol = ATOL
+    return
+end
+
+function test_reverse_vector_constraint_coefficient()
+    # Test JuMP.coefficient(func, vi, output_index) on vector constraint gradients
+    # SOC path: returns IndexMappedFunction{MatrixVectorAffineFunction}
+    model = DiffOpt.diff_model(SCS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("verbose"), false)
+    @variable(model, x)
+    @variable(model, y)
+    @variable(model, t)
+    @constraint(model, ceq, -1.0t == -1.0)
+    @constraint(model, csoc, [1.0t, 1.0x, 1.0y] in MOI.SecondOrderCone(3))
+    @constraint(model, cnon, 1.0y >= 1 / √2)
+    @objective(model, Min, 1.0x)
+    optimize!(model)
+    MOI.set(
+        JuMP.backend(model),
+        DiffOpt.ReverseVariablePrimal(),
+        JuMP.index(y),
+        1.0,
+    )
+    DiffOpt.reverse_differentiate!(model)
+    grad_soc = MOI.get(
+        JuMP.backend(model),
+        DiffOpt.ReverseConstraintFunction(),
+        JuMP.index(csoc),
+    )
+    # Verify per-output-index coefficient extraction
+    for i in 1:3
+        coeff = JuMP.coefficient(grad_soc, JuMP.index(y), i)
+        @test isfinite(coeff)
+    end
+    # PSD bridge path: returns IndexMappedFunction{VectorAffineFunction}
+    model2 = DiffOpt.diff_model(SCS.Optimizer)
+    MOI.set(model2, MOI.RawOptimizerAttribute("verbose"), false)
+    MOI.set(model2, MOI.RawOptimizerAttribute("eps_abs"), 1e-8)
+    MOI.set(model2, MOI.RawOptimizerAttribute("eps_rel"), 1e-8)
+    @variable(model2, s >= 0)
+    cons = @constraint(model2, [s 1.0; 1.0 s] in PSDCone())
+    @objective(model2, Min, s)
+    optimize!(model2)
+    MOI.set(
+        JuMP.backend(model2),
+        DiffOpt.ReverseVariablePrimal(),
+        JuMP.index(s),
+        1.0,
+    )
+    DiffOpt.reverse_differentiate!(model2)
+    grad_psd = MOI.get(
+        JuMP.backend(model2),
+        DiffOpt.ReverseConstraintFunction(),
+        JuMP.index(cons),
+    )
+    # Square col-major: [a11, a21, a12, a22], s appears in a11 and a22
+    coeff_s_1 = JuMP.coefficient(grad_psd, JuMP.index(s), 1)
+    coeff_s_4 = JuMP.coefficient(grad_psd, JuMP.index(s), 4)
+    @test coeff_s_1 ≈ -0.5 atol = ATOL
+    @test coeff_s_4 ≈ -0.5 atol = ATOL
+    return
+end
+
+function test_square_bridge_forward_psd()
+    model = DiffOpt.diff_model(SCS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("verbose"), false)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1e-8)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_rel"), 1e-8)
+    @variable(model, t >= 0)
+    cons = @constraint(model, [t 1.0; 1.0 t] in PSDCone())
+    @objective(model, Min, t)
+    optimize!(model)
+    @test value(t) ≈ 1.0 atol = ATOL
+    # Perturb the off-diagonal entries
+    DiffOpt.set_forward_constraint_function(model, cons, [0.0 1.0; 1.0 0.0])
+    DiffOpt.forward_differentiate!(model)
+    dt = DiffOpt.get_forward_variable(model, t)
+    @test dt ≈ 0.5 atol = ATOL
+    # Perturb the diagonal entries
+    DiffOpt.set_forward_constraint_function(model, cons, [1.0 0.0; 0.0 1.0])
+    DiffOpt.forward_differentiate!(model)
+    dt2 = DiffOpt.get_forward_variable(model, t)
+    @test dt2 ≈ -1.0 atol = ATOL
+    return
+end
+
+function test_square_bridge_reverse_psd()
+    model = DiffOpt.diff_model(SCS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("verbose"), false)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1e-8)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_rel"), 1e-8)
+    @variable(model, t >= 0)
+    cons = @constraint(model, [t 1.0; 1.0 t] in PSDCone())
+    @objective(model, Min, t)
+    optimize!(model)
+    @test value(t) ≈ 1.0 atol = ATOL
+    # Reverse mode: set dt = 1
+    MOI.set(
+        JuMP.backend(model),
+        DiffOpt.ReverseVariablePrimal(),
+        JuMP.index(t),
+        1.0,
+    )
+    DiffOpt.reverse_differentiate!(model)
+    grad = MOI.get(
+        JuMP.backend(model),
+        DiffOpt.ReverseConstraintFunction(),
+        JuMP.index(cons),
+    )
+    grad_cte = MOI.constant(grad)
+    # Result is in square col-major form: [a11, a21, a12, a22]
+    # The bridge mirrors triangle → square, so result should be symmetric
+    @test grad_cte[1] ≈ grad_cte[4] atol = ATOL  # a11 ≈ a22 (diagonal symmetry)
+    @test grad_cte[2] ≈ grad_cte[3] atol = ATOL  # a21 ≈ a12 (off-diagonal mirrored)
+    # Diagonal entries: forward perturbing [1 0; 0 0] gives dt ≈ -0.5
+    @test grad_cte[1] ≈ -0.5 atol = ATOL
+    @test grad_cte[4] ≈ -0.5 atol = ATOL
+    # Off-diagonal: forward perturbing [0 1; 1 0] gives dt ≈ 0.5
+    # Triangle off-diagonal gradient ≈ 1.0 (known PSD triangle scaling behavior)
+    @test grad_cte[2] ≈ 1.0 atol = ATOL
+    @test grad_cte[3] ≈ 1.0 atol = ATOL
+    return
+end
+
+function test_square_bridge_asymmetric_sym()
+    # Tests SquareBridge with non-empty bridge.sym (different off-diagonal variables).
+    # Similar to MOI's ConstraintPrimalStart/ConstraintDual on SquareBridge.
+    #
+    # min x + y  s.t. [2 x; y 2] ∈ PSD
+    # Bridge creates: [2, x, 2] ∈ PSDTriangle(2) and [x - y] ∈ Zeros(1)
+    # So x == y (implicit), eigenvalues are 2+x and 2-x, both >= 0 → -2 ≤ x ≤ 2
+    # min 2x → x = y = -2, obj = -4
+    # Perturbing diagonal by ε: [2+ε, x, 2+ε] → x = -(2+ε), dx/dε = dy/dε = -1
+    model = DiffOpt.diff_model(SCS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("verbose"), false)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1e-8)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_rel"), 1e-8)
+    @variable(model, x)
+    @variable(model, y)
+    @constraint(model, cpsd, [2.0 1.0*x; 1.0*y 2.0] in PSDCone())
+    @objective(model, Min, x + y)
+    optimize!(model)
+    @test value(x) ≈ -2.0 atol = ATOL
+    @test value(y) ≈ -2.0 atol = ATOL
+    # Forward: perturb diagonal by 1 → dx/dε = dy/dε = -1
+    DiffOpt.set_forward_constraint_function(model, cpsd, [1.0 0.0; 0.0 1.0])
+    DiffOpt.forward_differentiate!(model)
+    @test DiffOpt.get_forward_variable(model, x) ≈ -1.0 atol = ATOL
+    @test DiffOpt.get_forward_variable(model, y) ≈ -1.0 atol = ATOL
+    # Reverse: seed dx=1, dy=1 → d(loss)/d(diag) = 1*(-1) + 1*(-1) = -2
+    DiffOpt.empty_input_sensitivities!(model)
+    MOI.set(
+        JuMP.backend(model),
+        DiffOpt.ReverseVariablePrimal(),
+        JuMP.index(x),
+        1.0,
+    )
+    MOI.set(
+        JuMP.backend(model),
+        DiffOpt.ReverseVariablePrimal(),
+        JuMP.index(y),
+        1.0,
+    )
+    DiffOpt.reverse_differentiate!(model)
+    rev = MOI.get(
+        JuMP.backend(model),
+        DiffOpt.ReverseConstraintFunction(),
+        JuMP.index(cpsd),
+    )
+    cte = MOI.constant(rev)
+    @test length(cte) == 4  # square form: [a11, a21, a12, a22]
+    # Diagonal entries: each ≈ -1.0 (follows ConstraintDual convention)
+    @test cte[1] ≈ -1.0 atol = ATOL
+    @test cte[4] ≈ -1.0 atol = ATOL
+    # Off-diagonal: ConstraintDual convention with sym (2η+π / -π)
+    # where η ≈ -1.0 (triangle off-diag gradient) and π ≈ 1.0 (sym dual)
+    @test cte[2] ≈ -1.0 atol = ATOL  # -π
+    @test cte[3] ≈ -3.0 atol = ATOL  # 2η + π
+    return
+end
+
+function test_square_bridge_asymmetric_sym_parametric()
+    # Same as test_square_bridge_asymmetric_sym but with parameter p on diagonal.
+    # min x + y  s.t. [p x; y p] ∈ PSD, p = 2
+    # Bridge: [p, x, p] ∈ PSDTriangle(2), [x - y] ∈ Zeros(1)
+    # Eigenvalues: p+x and p-x, both ≥ 0 → -p ≤ x ≤ p
+    # min 2x → x = y = -p = -2
+    # dx/dp = dy/dp = -1
+    model = DiffOpt.conic_diff_model(SCS.Optimizer)
+    MOI.set(model, MOI.RawOptimizerAttribute("verbose"), false)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1e-8)
+    MOI.set(model, MOI.RawOptimizerAttribute("eps_rel"), 1e-8)
+    @variable(model, x)
+    @variable(model, y)
+    @variable(model, p in Parameter(2.0))
+    @constraint(model, cpsd, [1.0*p 1.0*x; 1.0*y 1.0*p] in PSDCone())
+    @objective(model, Min, x + y)
+    optimize!(model)
+    @test value(x) ≈ -2.0 atol = ATOL
+    @test value(y) ≈ -2.0 atol = ATOL
+    # Forward: dp = 1 → dx = dy = -1
+    DiffOpt.set_forward_parameter(model, p, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    @test DiffOpt.get_forward_variable(model, x) ≈ -1.0 atol = ATOL
+    @test DiffOpt.get_forward_variable(model, y) ≈ -1.0 atol = ATOL
+    # Reverse: seed dx=1 → dp = -1
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_reverse_variable(model, x, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    @test DiffOpt.get_reverse_parameter(model, p) ≈ -1.0 atol = ATOL
+    # Reverse: seed dy=1 → dp = -1
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_reverse_variable(model, y, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    @test DiffOpt.get_reverse_parameter(model, p) ≈ -1.0 atol = ATOL
+    # Reverse: seed dx=1, dy=1 → dp = -2
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_reverse_variable(model, x, 1.0)
+    DiffOpt.set_reverse_variable(model, y, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    @test DiffOpt.get_reverse_parameter(model, p) ≈ -2.0 atol = ATOL
+    return
 end
 
 end  # module
diff --git a/test/jump.jl b/test/jump.jl
index 059e6e0b..e89d1b73 100644
--- a/test/jump.jl
+++ b/test/jump.jl
@@ -46,7 +46,7 @@ function test_single_variable_objective_forward()
     )
     @objective(model, Max, x[7])
     optimize!(model)
-    MOI.set(model, DiffOpt.ForwardObjectiveFunction(), sum(x))
+    DiffOpt.set_forward_objective_function(model, sum(x))
     DiffOpt.forward_differentiate!(model)
     @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x[7]) ≈ 0 atol = ATOL
     return
@@ -69,7 +69,7 @@ function test_single_variable_objective_reverse()
     optimize!(model)
     MOI.set(model, DiffOpt.ReverseVariablePrimal(), x[7], 1.0)
     DiffOpt.reverse_differentiate!(model)
-    func = MOI.get(model, DiffOpt.ReverseObjectiveFunction())
+    func = DiffOpt.get_reverse_objective_function(model)
     @test JuMP.coefficient(func, x[7]) ≈ 0.0 atol = ATOL rtol = RTOL
     return
 end
@@ -112,12 +112,11 @@ function test_differentiating_trivial_qp_1()
     MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x, 1.0)
     DiffOpt.reverse_differentiate!(model)
     DiffOpt.reverse_differentiate!(model)
-    grad_constraint = JuMP.constant(
-        MOI.get(model, DiffOpt.ReverseConstraintFunction(), ctr_le[1]),
-    )
+    grad_constraint =
+        JuMP.constant(DiffOpt.get_reverse_constraint_function(model, ctr_le[1]))
     @test grad_constraint ≈ -1.0 atol = ATOL rtol = RTOL
     # Test some overloads from https://github.com/jump-dev/DiffOpt.jl/issues/211
-    grad_obj = MOI.get(model, DiffOpt.ReverseObjectiveFunction())
+    grad_obj = DiffOpt.get_reverse_objective_function(model)
     @test JuMP.coefficient(grad_obj, x[1], x[2]) ≈
           DiffOpt.quad_sym_half.(grad_obj, x[1], x[2]) atol = ATOL rtol = RTOL
     @test DiffOpt.quad_sym_half(grad_obj, x[1], x[1]) ≈
@@ -720,19 +719,57 @@ function test_conic_feasibility()
     return
 end
 
-function test_psd_square_error()
+function test_psd_square()
+    # min -x  s.t. [p-x 0; 0 x] ∈ PSD
+    # PSD requires 0 ≤ x ≤ p, so x* = p = 1, dx/dp = 1
     model = DiffOpt.conic_diff_model(SCS.Optimizer)
     set_silent(model)
 
     @variable(model, x)
     @variable(model, p in Parameter(1.0))
+    @objective(model, Min, -x)
+    @constraint(model, con, [p - x 0; 0 x] in PSDCone())
+
+    optimize!(model)
+    @test is_solved_and_feasible(model)
+    @test value(x) ≈ 1.0 atol = ATOL rtol = RTOL
+
+    DiffOpt.set_forward_parameter(model, p, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    @test DiffOpt.get_forward_variable(model, x) ≈ 1.0 atol = ATOL rtol = RTOL
+    @test DiffOpt.get_forward_objective(model) ≈ -1.0 atol = ATOL rtol = RTOL
+
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_reverse_objective(model, 1.0)
+    DiffOpt.reverse_differentiate!(model)
+    @test DiffOpt.get_reverse_parameter(model, p) ≈ -1.0 atol = ATOL rtol = RTOL
+    return
+end
+
+function test_nlp_forward_constraint_dual()
+    # min 2x² + y² + xy + x + y  s.t. x + y == p, x ≥ 0, y ≥ 0
+    # KKT gives x = p/4, y = 3p/4, dual = 7p/4 + 1
+    # At p=1: dx/dp = 1/4, dy/dp = 3/4, d(dual)/dp = 7/4
+    model = DiffOpt.nonlinear_diff_model(Ipopt.Optimizer)
+    set_silent(model)
 
-    @constraint(model, con, [-p*x 0; 0 x] in PSDCone())
+    @variable(model, x >= 0)
+    @variable(model, y >= 0)
+    @variable(model, p in Parameter(1.0))
 
-    @test_throws MOI.Bridges.ModifyBridgeNotAllowed optimize!(model)
+    @constraint(model, c1, x + y == p)
+    @objective(model, Min, 2x^2 + y^2 + x * y + x + y)
+    optimize!(model)
+    @test value(x) ≈ 0.25 atol = ATOL rtol = RTOL
+    @test value(y) ≈ 0.75 atol = ATOL rtol = RTOL
 
-    # DiffOpt.set_forward_parameter(model, p, 1.0)
-    # DiffOpt.forward_differentiate!(model)
+    DiffOpt.set_forward_parameter(model, p, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    @test DiffOpt.get_forward_variable(model, x) ≈ 0.25 atol = ATOL rtol = RTOL
+    @test DiffOpt.get_forward_variable(model, y) ≈ 0.75 atol = ATOL rtol = RTOL
+    @test DiffOpt.get_forward_constraint_dual(model, c1) ≈ 1.75 atol = ATOL rtol =
+        RTOL
+    @test DiffOpt.get_forward_objective(model) ≈ 2.75 atol = ATOL rtol = RTOL
     return
 end
 
diff --git a/test/jump_wrapper.jl b/test/jump_wrapper.jl
index 0f2cae5e..43d531f0 100644
--- a/test/jump_wrapper.jl
+++ b/test/jump_wrapper.jl
@@ -11,6 +11,7 @@ import DiffOpt
 import HiGHS
 import Ipopt
 import SCS
+import LinearAlgebra
 import MathOptInterface as MOI
 
 const ATOL = 1e-3
@@ -337,6 +338,91 @@ function test_jump_api()
     return
 end
 
+function test_forward_wrappers_non_parametric()
+    # Tests set_forward_objective_function and set_forward_constraint_function
+    # overloads on a non-parametric model.
+    #
+    # min x  s.t.  x + y == 1, x >= 0, y >= 0
+    # Solution: x=0, y=1
+    # Perturb constraint RHS: x + y + ϵ == 1
+    # Sensitivity should be dy/dϵ = -1
+    # since the slack goes to y and x is at its lower bound (active).
+    model = JuMP.direct_model(DiffOpt.diff_optimizer(HiGHS.Optimizer))
+    set_silent(model)
+    @variable(model, x >= 0)
+    @variable(model, y >= 0)
+    @constraint(model, c1, x + y == 1)
+    @objective(model, Min, 1.0 * x)
+    optimize!(model)
+    @test value(x) ≈ 0.0 atol = ATOL
+    @test value(y) ≈ 1.0 atol = ATOL
+
+    DiffOpt.set_forward_objective_function(model, 0.0)
+    DiffOpt.set_forward_constraint_function(model, c1, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dy = DiffOpt.get_forward_variable(model, y)
+    @test dy ≈ -1.0 atol = ATOL
+
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_forward_constraint_function(model, c1, 0.0 * x + 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dy2 = DiffOpt.get_forward_variable(model, y)
+    @test dy2 ≈ dy atol = ATOL
+    return
+end
+
+function test_forward_vector_constraint_wrappers()
+    # Tests set_forward_constraint_function overloads for vector constraints
+    # using conic_diff_model (direct_model + diff_optimizer bridges differently).
+    model = DiffOpt.conic_diff_model(SCS.Optimizer)
+    set_silent(model)
+    @variable(model, x)
+    @variable(model, y)
+    @constraint(model, c_eq, x + y == 1)
+    @constraint(model, c_nn, [1.0 * y, 1.0 * x] in MOI.Nonnegatives(2))
+    @objective(model, Min, 1.0 * x)
+    optimize!(model)
+    @test value(x) ≈ 0.0 atol = ATOL
+    @test value(y) ≈ 1.0 atol = ATOL
+
+    DiffOpt.set_forward_constraint_function(model, c_nn, [0.0, 0.0])
+    DiffOpt.set_forward_constraint_function(model, c_eq, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dy = DiffOpt.get_forward_variable(model, y)
+    @test dy ≈ -1.0 atol = ATOL
+
+    DiffOpt.empty_input_sensitivities!(model)
+    DiffOpt.set_forward_constraint_function(model, c_nn, [0.0 * x, 0.0 * x])
+    DiffOpt.set_forward_constraint_function(model, c_eq, 1.0)
+    DiffOpt.forward_differentiate!(model)
+    dy2 = DiffOpt.get_forward_variable(model, y)
+    @test dy2 ≈ dy atol = ATOL
+    return
+end
+
+function test_forward_psd_matrix_wrapper()
+    # Tests set_forward_constraint_function with AbstractMatrix{<:AbstractJuMPScalar}
+    # for PSD cone constraints. Uses a non-parametric model since
+    # ForwardConstraintFunction is blocked on parametric models.
+    model = Model(() -> DiffOpt.diff_optimizer(SCS.Optimizer))
+    set_silent(model)
+    @variable(model, x)
+    @objective(model, Min, -x)
+    @constraint(
+        model,
+        con,
+        LinearAlgebra.Symmetric([1.0-x 0.0; 0.0 x]) in PSDCone(),
+    )
+    optimize!(model)
+    @test value(x) ≈ 1.0 atol = ATOL
+
+    perturbation = [1.0+0.0*x 0.0*x; 0.0*x 0.0*x]
+    DiffOpt.set_forward_constraint_function(model, con, perturbation)
+    DiffOpt.forward_differentiate!(model)
+    @test DiffOpt.get_forward_variable(model, x) ≈ 1.0 atol = ATOL
+    return
+end
+
 end # module
 
 TestJuMPWrapper.runtests()