diff --git a/examples/JuMP_utils.jl b/examples/JuMP_utils.jl
index ab8cd21f0..f4cd7be70 100644
--- a/examples/JuMP_utils.jl
+++ b/examples/JuMP_utils.jl
@@ -184,8 +184,8 @@ function solve_check(
             Hypatia.rescale_affine(cone, z_k)
         end
         if Hypatia.needs_permute(cone)
-            Hypatia.permute_affine(cone, s_k)
-            Hypatia.permute_affine(cone, z_k)
+            s_k = Hypatia.permute_affine(cone, s_k)
+            z_k = Hypatia.permute_affine(cone, z_k)
         end
         append!(s, s_k)
         append!(z, z_k)
diff --git a/src/MathOptInterface/cones.jl b/src/MathOptInterface/cones.jl
index 6ee7293df..902202f99 100644
--- a/src/MathOptInterface/cones.jl
+++ b/src/MathOptInterface/cones.jl
@@ -83,8 +83,8 @@ end
 needs_permute(cone::SpecNucCone) = needs_untransform(cone)
 
 function permute_affine(cone::SpecNucCone, vals::AbstractVector{T}) where {T}
-    @views vals[2:end] = reshape(vals[2:end], cone.row_dim, cone.column_dim)'
-    return vals
+    w_vals = reshape(vals[2:end], cone.row_dim, cone.column_dim)'
+    return vcat(vals[1], vec(w_vals))
 end
 
 function permute_affine(cone::SpecNucCone, func::VAF{T}) where {T}
diff --git a/src/MathOptInterface/wrapper.jl b/src/MathOptInterface/wrapper.jl
index 948d536e3..0e8cfd8b8 100644
--- a/src/MathOptInterface/wrapper.jl
+++ b/src/MathOptInterface/wrapper.jl
@@ -212,6 +212,62 @@ end
 
 MOI.optimize!(opt::Optimizer) = Solvers.solve(opt.solver)
 
+function MOI.modify(
+    opt::Optimizer{T},
+    ::MOI.ObjectiveFunction{SAF{T}},
+    chg::MOI.ScalarConstantChange{T},
+    ) where {T}
+    obj_offset = chg.new_constant
+    if opt.obj_sense == MOI.MAX_SENSE
+        obj_offset = -obj_offset
+    end
+    Solvers.modify_obj_offset(opt.solver, obj_offset)
+    return
+end
+
+function MOI.modify(
+    opt::Optimizer{T},
+    ::MOI.ObjectiveFunction{SAF{T}},
+    chg::MOI.ScalarCoefficientChange{T},
+    ) where {T}
+    new_c = chg.new_coefficient
+    if opt.obj_sense == MOI.MAX_SENSE
+        new_c = -new_c
+    end
+    Solvers.modify_c(opt.solver, [chg.variable.value], [new_c])
+    return
+end
+
+function MOI.modify(
+    opt::Optimizer{T},
+    ci::MOI.ConstraintIndex{VAF{T}, MOI.Zeros},
+    chg::MOI.VectorConstantChange{T},
+    ) where {T}
+    idxs = opt.zeros_idxs[ci.value]
+    Solvers.modify_b(opt.solver, idxs, chg.new_constant)
+    return
+end
+
+function MOI.modify(
+    opt::Optimizer{T},
+    ci::MOI.ConstraintIndex{VAF{T}, <:SupportedCone{T}},
+    chg::MOI.VectorConstantChange{T},
+    ) where {T}
+    i = ci.value
+    idxs = opt.moi_cone_idxs[i]
+    set = opt.moi_cones[i]
+    new_h = chg.new_constant
+    if needs_permute(set)
+        @assert !needs_rescale(set)
+        new_h = permute_affine(set, new_h)
+    end
+    if needs_rescale(set)
+        rescale_affine(set, new_h)
+    end
+    Solvers.modify_h(opt.solver, idxs, new_h)
+    return
+end
+
 function MOI.get(opt::Optimizer, ::MOI.TerminationStatus)
     status = opt.solver.status
     if status in (Solvers.NotLoaded, Solvers.Loaded)
diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl
index 7966512a2..7a81c93ff 100644
--- a/src/Solvers/Solvers.jl
+++ b/src/Solvers/Solvers.jl
@@ -34,6 +34,7 @@ include("point.jl")
 @enum Status begin
     NotLoaded
     Loaded
+    Modified
     SolveCalled
     PrimalInconsistent
     DualInconsistent
@@ -110,8 +111,13 @@ mutable struct Solver{T <: Real}
     model::Models.Model{T}
     x_keep_idxs::AbstractVector{Int}
     y_keep_idxs::AbstractVector{Int}
+    Ap_fact::Factorization{T}
+    Ap_rank::Int
     Ap_R::UpperTriangular{T, <:AbstractMatrix{T}}
     Ap_Q::Union{UniformScaling, AbstractMatrix{T}}
+    AG_fact::Factorization{T}
+    AG_rank::Int
+    AG_R::UpperTriangular{T, <:AbstractMatrix{T}}
     reduce_cQ1
     reduce_Rpib0
     reduce_GQ1
@@ -197,7 +203,7 @@ mutable struct Solver{T <: Real}
         if reduce
             @assert preprocess # cannot use reduction without preprocessing # TODO only need primal eq preprocessing
         end
-        # @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point
+        @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point
         @assert near_factor >= 1 # factor to relax tolerances by if fail to converge should be at least 1
 
         if isnothing(default_tol_power)
@@ -223,7 +229,7 @@ mutable struct Solver{T <: Real}
             tol_infeas = default_tol_tight
         end
         if isnothing(tol_illposed)
-            tol_illposed = default_tol_tight / 100
+            tol_illposed = default_tol_tight
         end
         @assert min(tol_rel_opt, tol_abs_opt, tol_feas,
             tol_infeas, tol_illposed, tol_slow) >= 0
@@ -254,13 +260,90 @@ mutable struct Solver{T <: Real}
     end
 end
 
+function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real}
+    solver.orig_model = model
+    solver.status = Loaded
+    return
+end
+
+function modify_obj_offset(solver::Solver{T}, offset_new::T) where {T <: Real}
+    solver.orig_model.obj_offset = offset_new
+    solver.status = Modified
+    return
+end
+
+function modify_c(solver::Solver{T}, c_new::AbstractVector{T}) where {T <: Real}
+    @assert length(c_new) == solver.orig_model.n
+    solver.orig_model.c = c_new
+    solver.status = Modified
+    return
+end
+
+function modify_c(
+    solver::Solver{T},
+    idxs::AbstractVector{Int},
+    c_new::Vector{T},
+    ) where {T <: Real}
+    @assert length(c_new) == length(idxs) <= solver.orig_model.n
+    solver.orig_model.c[idxs] .= c_new
+    solver.status = Modified
+    return
+end
+
+function modify_b(solver::Solver{T}, b_new::AbstractVector{T}) where {T <: Real}
+    @assert length(b_new) == solver.orig_model.p
+    solver.orig_model.b = b_new
+    solver.status = Modified
+    return
+end
+
+function modify_b(
+    solver::Solver{T},
+    idxs::AbstractVector{Int},
+    b_new::Vector{T},
+    ) where {T <: Real}
+    @assert length(b_new) == length(idxs) <= solver.orig_model.p
+    solver.orig_model.b[idxs] .= b_new
+    solver.status = Modified
+    return
+end
+
+function modify_h(solver::Solver{T}, h_new::AbstractVector{T}) where {T <: Real}
+    @assert length(h_new) == solver.orig_model.q
+    solver.orig_model.h = h_new
+    solver.status = Modified
+    return
+end
+
+function modify_h(
+    solver::Solver{T},
+    idxs::AbstractVector{Int},
+    h_new::Vector{T},
+    ) where {T <: Real}
+    @assert length(h_new) == length(idxs) <= solver.orig_model.q
+    solver.orig_model.h[idxs] .= h_new
+    solver.status = Modified
+    return
+end
+
 function solve(solver::Solver{T}) where {T <: Real}
-    (solver.status == Loaded) || return solver
+    init_status = solver.status
+    if !in(init_status, (Loaded, Modified))
+        @warn("solve called when solver status is $init_status")
+        return
+    end
+
     setup_solver(solver)
-    setup_point(solver)
+    if init_status == Modified
+        setup_modified(solver)
+    else
+        setup_loaded(solver)
+    end
 
-    if solver.status == SolveCalled
+    if !in(solver.status, (PrimalInconsistent, DualInconsistent))
+        initialize_point(solver)
         setup_stepping(solver)
+        solver.verbose && print_header(solver.stepper, solver)
 
         while true
             step_and_check(solver) && break
@@ -281,7 +364,132 @@ function solve(solver::Solver{T}) where {T <: Real}
 
     free_memory(solver.syssolver)
     flush(stdout)
-    return solver
+    return
+end
+
+function setup_loaded(solver::Solver{T}) where {T <: Real}
+    orig_model = solver.orig_model
+    if solver.use_dense_model
+        Models.densify!(orig_model)
+    end
+    solver.result = Point(orig_model)
+
+    # copy original model to solver.model, which may be modified
+    model = solver.model = Models.Model{T}(
+        orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h,
+        orig_model.cones, obj_offset = orig_model.obj_offset)
+
+    solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver)
+
+    setup_point(solver)
+    if solver.status != SolveCalled
+        return
+    end
+
+    load(solver.stepper, solver)
+    solver.time_loadsys = @elapsed load(solver.syssolver, solver)
+
+    solver.x_residual = zero(model.c)
+    solver.y_residual = zero(model.b)
+    solver.z_residual = zero(model.h)
+    solver.tau_residual = 0
+    return
+end
+
+function setup_modified(solver::Solver{T}) where {T <: Real}
+    orig_model = solver.orig_model
+    model = solver.model
+
+    model.obj_offset = orig_model.obj_offset
+    model.c = copy(orig_model.c)
+    model.b = copy(orig_model.b)
+    model.h = copy(orig_model.h)
+    if solver.used_rescaling
+        # rescale using scalings computed during first load
+        model.c ./= solver.c_scale
+        model.b ./= solver.b_scale
+        model.h ./= solver.h_scale
+    end
+
+    update_initial_point(solver)
+
+    # TODO easier to update these system solvers (just need an update to setup_point_sub):
+    @assert solver.syssolver isa Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}}
+    set_point_sub_rhs(solver.syssolver, solver.model)
+    return
+end
+
+# preprocess and find initial point
+function setup_point(solver::Solver{T}) where {T <: Real}
+    (init_z, init_s) = initialize_cone_point(solver.orig_model)
+
+    if solver.reduce
+        solver.time_inity = @elapsed handle_primal_eq(solver)
+        solver.status == PrimalInconsistent && return
+        init_y = update_primal_eq(solver, init_z)
+
+        solver.time_initx = @elapsed handle_dual_eq(solver)
+        solver.status == DualInconsistent && return
+        init_x = update_dual_eq(solver, init_s)
+    else
+        solver.time_initx = @elapsed handle_dual_eq(solver)
+        solver.status == DualInconsistent && return
+        init_x = update_dual_eq(solver, init_s)
+
+        solver.time_inity = @elapsed handle_primal_eq(solver)
+        solver.status == PrimalInconsistent && return
+        init_y = update_primal_eq(solver, init_z)
+    end
+
+    point = solver.point = Point(solver.model)
+    point.x .= init_x
+    point.y .= init_y
+    point.z .= init_z
+    point.s .= init_s
+    return
+end
+
+function update_initial_point(solver::Solver{T}) where {T <: Real}
+    (init_z, init_s) = initialize_cone_point(solver.orig_model)
+    point = solver.point
+    if solver.reduce
+        point.y .= update_primal_eq(solver, init_z)
+        point.x .= update_dual_eq(solver, init_s)
+    else
+        point.x .= update_dual_eq(solver, init_s)
+        point.y .= update_primal_eq(solver, init_z)
+    end
+    point.z .= init_z
+    point.s .= init_s
+    return
+end
+
+function initialize_point(solver::Solver{T}) where {T <: Real}
+    point = solver.point
+    point.tau[] = one(T)
+    point.kap[] = one(T)
+
+    calc_mu(solver)
+    if isnan(solver.mu) || abs(one(T) - solver.mu) > sqrt(eps(T))
+        @warn("initial mu is $(solver.mu) but should be 1 (this could " *
+            "indicate a problem with cone barrier oracles)")
+    end
+
+    model = solver.model
+    Cones.load_point.(model.cones, point.primal_views)
+    Cones.load_dual_point.(model.cones, point.dual_views)
+    return
+end
+
+function setup_stepping(solver::Solver{T}) where {T <: Real}
+    model = solver.model
+    solver.x_conv_tol = inv(1 + norm(model.c, Inf))
+    solver.y_conv_tol = inv(1 + norm(model.b, Inf))
+    solver.z_conv_tol = inv(1 + norm(model.h, Inf))
+    solver.prev_is_slow = false
+    solver.prev2_is_slow = false
+    solver.worst_dir_res = 0
+    return
 end
 
 function setup_solver(solver::Solver{T}) where {T <: Real}
@@ -318,75 +526,6 @@ function setup_solver(solver::Solver{T}) where {T <: Real}
     solver.y_feas = NaN
     solver.z_feas = NaN
     solver.tau_feas = NaN
-
-    orig_model = solver.orig_model
-    if solver.use_dense_model
-        Models.densify!(orig_model)
-    end
-    solver.result = Point(orig_model)
-
-    # copy original model to solver.model, which may be modified
-    solver.model = Models.Model{T}(
-        orig_model.c, orig_model.A, orig_model.b, orig_model.G, orig_model.h,
-        orig_model.cones, obj_offset = orig_model.obj_offset)
-    return
-end
-
-function setup_point(solver::Solver{T}) where {T <: Real}
-    # preprocess and find initial point
-    (init_z, init_s) = initialize_cone_point(solver.orig_model)
-    solver.time_rescale = @elapsed solver.used_rescaling = rescale_data(solver)
-
-    if solver.reduce
-        # TODO don't find point / unnecessary stuff before reduce
-        solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, true)
-        solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s)
-    else
-        solver.time_initx = @elapsed init_x = find_initial_x(solver, init_s)
-        solver.time_inity = @elapsed init_y = find_initial_y(solver, init_z, false)
-    end
-
-    if solver.status == SolveCalled
-        model = solver.model
-        point = solver.point = Point(model)
-        point.x .= init_x
-        point.y .= init_y
-        point.z .= init_z
-        point.s .= init_s
-        point.tau[] = one(T)
-        point.kap[] = one(T)
-
-        calc_mu(solver)
-        if isnan(solver.mu) || abs(one(T) - solver.mu) > sqrt(eps(T))
-            @warn("initial mu is $(solver.mu) but should be 1 (this could " *
-                "indicate a problem with cone barrier oracles)")
-        end
-        Cones.load_point.(model.cones, point.primal_views)
-        Cones.load_dual_point.(model.cones, point.dual_views)
-    end
-    return
-end
-
-function setup_stepping(solver::Solver{T}) where {T <: Real}
-    model = solver.model
-    # setup iteration helpers
-    solver.x_residual = zero(model.c)
-    solver.y_residual = zero(model.b)
-    solver.z_residual = zero(model.h)
-    solver.tau_residual = 0
-
-    solver.x_conv_tol = inv(1 + norm(model.c, Inf))
-    solver.y_conv_tol = inv(1 + norm(model.b, Inf))
-    solver.z_conv_tol = inv(1 + norm(model.h, Inf))
-    solver.prev_is_slow = false
-    solver.prev2_is_slow = false
-    solver.worst_dir_res = 0
-
-    load(solver.stepper, solver)
-    solver.time_loadsys = @elapsed load(solver.syssolver, solver)
-
-    solver.verbose && print_header(solver.stepper, solver)
-    flush(stdout)
     return
 end
 
@@ -396,7 +535,6 @@ function step_and_check(solver::Solver{T}) where {T <: Real}
 
     if solver.verbose
         print_iteration(stepper, solver)
-        flush(stdout)
     end
 
     check_converged(solver) && return true
@@ -564,8 +702,7 @@ function check_converged(
     end
 
     # TODO experiment with ill-posedness check
-    max_illp = max(solver.mu, tau / min(one(T), solver.point.kap[]))
-    if max_illp <= near_factor * solver.tol_illposed
+    if max(tau, solver.point.kap[]) <= near_factor * solver.tol_illposed
         solver.verbose && println("ill-posedness detected; terminating")
         solver.status = (check_near ? NearIllPosed : IllPosed)
         return true
@@ -594,6 +731,8 @@ function initialize_cone_point(model::Models.Model{T}) where {T <: Real}
     return (init_z, init_s)
 end
 
+get_model(solver::Solver) = solver.orig_model
+
 get_status(solver::Solver) = solver.status
 get_solve_time(solver::Solver) = solver.solve_time
 get_num_iters(solver::Solver) = solver.num_iters
@@ -610,12 +749,6 @@ get_tau(solver::Solver) = solver.point.tau[]
 get_kappa(solver::Solver) = solver.point.kap[]
 get_mu(solver::Solver) = solver.mu
 
-function load(solver::Solver{T}, model::Models.Model{T}) where {T <: Real}
-    solver.orig_model = model
-    solver.status = Loaded
-    return solver
-end
-
 include("process.jl")
 
 include("search.jl")
@@ -642,6 +775,7 @@ function print_header(stepper::Stepper, solver::Solver)
 
     print_header_more(stepper, solver)
     println()
+    flush(stdout)
     return
 end
 print_header_more(stepper::Stepper, solver::Solver) = nothing
@@ -662,6 +796,7 @@ function print_iteration(stepper::Stepper, solver::Solver)
         print_iteration_more(stepper, solver)
     end
     println()
+    flush(stdout)
     return
 end
 print_iteration_more(stepper::Stepper, solver::Solver) = nothing
diff --git a/src/Solvers/process.jl b/src/Solvers/process.jl
index 8be11cd97..f28baf8e2 100644
--- a/src/Solvers/process.jl
+++ b/src/Solvers/process.jl
@@ -14,7 +14,9 @@ function rescale_data(solver::Solver{T}) where {T <: Real}
     solver.rescale || return false
     model = solver.model
     (c, A, b, G, h) = (model.c, model.A, model.b, model.G, model.h)
-    (isa(A, MatrixyAG) && isa(G, MatrixyAG)) || return false
+    if !isa(A, MatrixyAG) || !isa(G, MatrixyAG)
+        return false
+    end
 
     minval = sqrt(eps(T))
     maxabsmin(v::AbstractVecOrMat) = mapreduce(abs, max, v; init = minval)
@@ -59,100 +61,81 @@ function rescale_data(solver::Solver{T}) where {T <: Real}
     return true
 end
 
-# optionally preprocess dual equalities and solve for x as least squares
-# solution to Ax = b, Gx = h - s
-function find_initial_x(
-    solver::Solver{T},
-    init_s::Vector{T},
-    ) where {T <: Real}
-    if solver.status != SolveCalled
-        return zeros(T, 0)
-    end
+function handle_dual_eq(solver::Solver{T}) where {T <: Real}
     model = solver.model
     n = model.n
-    if iszero(n) # x is empty (no primal variables)
-        solver.x_keep_idxs = Int[]
-        return zeros(T, 0)
-    end
-    p = model.p
-    q = model.q
-    A = model.A
-    G = model.G
     solver.x_keep_idxs = 1:n
+    iszero(n) && return
 
-    rhs = vcat(model.b, model.h - init_s)
-
-    # indirect method
+    A = model.A
+    G = model.G
     if solver.init_use_indirect || !isa(A, MatrixyAG) || !isa(G, MatrixyAG)
-        # TODO pick lsqr or lsmr
-        if iszero(p)
-            AG = G
-        else
-            linmap(M::UniformScaling) = LinearMaps.LinearMap(M, n)
-            linmap(M) = LinearMaps.LinearMap(M)
-            AG = vcat(linmap(A), linmap(G))
-        end
-        init_x = IterativeSolvers.lsqr(AG, rhs)
-        return init_x
+        return
     end
 
     # direct method
-    if iszero(p)
+    AG = if iszero(model.p)
         # A is empty
-        if issparse(G)
-            AG = G
+        if G isa UniformScaling
+            sparse(G, n, n)
+        elseif issparse(G)
+            G
         elseif G isa Matrix{T}
-            AG = copy(G)
+            copy(G)
         else
-            AG = Matrix(G)
+            Matrix(G)
         end
     else
-        AG = vcat(A, G)
+        vcat(A, G)
     end
-    if issparse(AG)
+    AG_fact = solver.AG_fact = if issparse(AG)
         if !(T <: Float64)
             @warn("using dense factorization of [A; G] in preprocessing and " *
                 "initial point finding because sparse factorization for number " *
-                "type $T is not supported by SuiteSparse packages")
-            AG_fact = qr!(Matrix(AG), ColumnNorm())
+                "type $T is not supported by SuiteSparse packages", maxlog = 1)
+            qr!(Matrix(AG), ColumnNorm())
         else
-            AG_fact = qr(AG, tol = solver.init_tol_qr)
+            qr(AG, tol = solver.init_tol_qr)
         end
     else
-        AG_fact = qr!(AG, ColumnNorm())
+        qr!(AG, ColumnNorm())
     end
-    AG_rank = get_rank_est(AG_fact, solver.init_tol_qr)
+    AG_rank = solver.AG_rank = get_rank_est(AG_fact, solver.init_tol_qr)
 
     if !solver.preprocess || (AG_rank == n)
         if AG_rank < n
             @warn("some dual equalities appear to be dependent " *
             "(possibly inconsistent); try using preprocess = true")
         end
-        init_x = AG_fact \ rhs
-        return init_x
+        return
     end
 
-    # preprocess dual equalities
+    # preprocess
+    AG_R = solver.AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank])
     col_piv = (AG_fact isa QRPivoted{T, Matrix{T}}) ? AG_fact.p : AG_fact.pcol
     x_keep_idxs = col_piv[1:AG_rank]
-    AG_R = UpperTriangular(AG_fact.R[1:AG_rank, 1:AG_rank])
 
-    c_sub = model.c[x_keep_idxs]
     # yz_sub = AG_fact.Q * vcat((AG_R' \ c_sub), zeros(p + q - AG_rank))
-    yz_sub = zeros(T, p + q)
+    p = model.p
+    yz_sub = zeros(T, p + model.q)
     yz_sub1 = view(yz_sub, 1:AG_rank)
-    copyto!(yz_sub1, c_sub)
+    @views copyto!(yz_sub1, model.c[x_keep_idxs])
     ldiv!(AG_R', yz_sub1)
     lmul!(AG_fact.Q, yz_sub)
     if !(AG_fact isa QRPivoted{T, Matrix{T}})
         yz_sub = yz_sub[AG_fact.rpivinv]
     end
-    @views residual = norm(A' * yz_sub[1:p] + G' *
-        yz_sub[(p + 1):end] - model.c, Inf)
-    if residual > solver.init_tol_qr
+
+    # residual = A' * yz_sub[1:p] + G' * yz_sub[(p + 1):end] - model.c
+    residual = copy(model.c)
+    @views mul!(residual, G', yz_sub[(p + 1):end], true, -1)
+    @views mul!(residual, A', yz_sub[1:p], true, true)
+    res_norm = norm(residual, Inf)
+
+    if res_norm > solver.init_tol_qr
         if solver.verbose
             println("some dual equality constraints are inconsistent " *
-            "(residual $residual, tolerance $(solver.init_tol_qr))")
+            "(residual norm $res_norm, tolerance $(solver.init_tol_qr))")
         end
         solver.status = DualInconsistent
         return zeros(T, 0)
@@ -161,102 +144,111 @@ function find_initial_x(
         println("$(n - AG_rank) of $n dual equality constraints are dependent")
     end
 
-    # modify solver.model to remove/reorder some primal variables x
-    model.c = c_sub
+    # modify model to remove/reorder some primal variables x
     model.A = A[:, x_keep_idxs]
     model.G = G[:, x_keep_idxs]
     model.n = AG_rank
     solver.x_keep_idxs = x_keep_idxs
+    return
+end
+
+# update data for b, c, h from dual equality preprocessing/reduction
+# and return initial x as least squares solution to Ax = b, Gx = h - s
+function update_dual_eq(
+    solver::Solver{T},
+    init_s::Vector{T},
+    ) where {T <: Real}
+    model = solver.model
+
+    model.c = model.c[solver.x_keep_idxs]
+    iszero(model.n) && return zeros(T, 0)
+
+    rhs_x = vcat(model.b, model.h - init_s)
+
+    if solver.init_use_indirect ||
+        !isa(model.A, MatrixyAG) || !isa(model.G, MatrixyAG)
+        # use indirect method TODO pick lsqr or lsmr
+        if iszero(model.p)
+            AG = model.G
+        else
+            linmap(M::UniformScaling) = LinearMaps.LinearMap(M, model.n)
+            linmap(M) = LinearMaps.LinearMap(M)
+            AG = vcat(linmap(model.A), linmap(model.G))
+        end
+        return IterativeSolvers.lsqr(AG, rhs_x)
+    end
+
+    AG_fact = solver.AG_fact
+    if !solver.preprocess || (solver.AG_rank == size(AG_fact, 2))
+        return AG_fact \ rhs_x
+    end
 
     # init_x = AG_R \ ((AG_fact.Q' * vcat(b, h - point.s))[1:AG_rank])
-    temp = vcat(model.b, model.h - init_s)
-    lmul!(AG_fact.Q', temp)
-    init_x = temp[1:model.n]
-    ldiv!(AG_R, init_x)
+    lmul!(AG_fact.Q', rhs_x)
+    init_x = rhs_x[1:model.n]
+    ldiv!(solver.AG_R, init_x)
 
     return init_x
 end
 
-# optionally preprocess primal equalities and solve for y as least squares
-# solution to A'y = -c - G'z
-function find_initial_y(
-    solver::Solver{T},
-    init_z::Vector{T},
-    reduce::Bool,
-    ) where {T <: Real}
-    if solver.status != SolveCalled
-        return zeros(T, 0)
-    end
+function handle_primal_eq(solver::Solver{T}) where {T <: Real}
     model = solver.model
     p = model.p
-    if iszero(p) # y is empty (no primal variables)
+
+    if iszero(p)
+        # y is empty (no primal variables)
         solver.y_keep_idxs = Int[]
         solver.Ap_R = UpperTriangular(zeros(T, 0, 0))
         solver.Ap_Q = I
         return zeros(T, 0)
     end
-    n = model.n
-    q = model.q
-    A = model.A
     solver.y_keep_idxs = 1:p
 
-    if !reduce || !isa(A, MatrixyAG)
-        # rhs = -c - G' * point.z
-        rhs = copy(model.c)
-        mul!(rhs, model.G', init_z, -1, -1)
-
-        # indirect method
-        if solver.init_use_indirect
-            # TODO pick lsqr or lsmr
-            init_y = IterativeSolvers.lsqr(A', rhs)
-            return init_y
-        end
-    end
+    solver.init_use_indirect && return
 
     # factorize A'
-    if issparse(A)
+    A = model.A
+    Ap_fact = solver.Ap_fact = if issparse(A)
         if !(T <: Float64)
             @warn("using dense factorization of A' in preprocessing and initial " *
                 "point finding because sparse factorization for number type $T " *
-                "is not supported by SuiteSparse packages")
-            Ap_fact = qr!(Matrix(A'), ColumnNorm())
+                "is not supported by SuiteSparse packages", maxlog = 1)
+            qr!(Matrix(A'), ColumnNorm())
         else
-            Ap_fact = qr(sparse(A'), tol = solver.init_tol_qr)
+            qr(sparse(A'), tol = solver.init_tol_qr)
         end
     else
-        Ap_fact = qr!(Matrix(A'), ColumnNorm())
+        qr!(Matrix(A'), ColumnNorm())
     end
-    Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr)
+    Ap_rank = solver.Ap_rank = get_rank_est(Ap_fact, solver.init_tol_qr)
 
-    if !reduce && !solver.preprocess
+    if !solver.preprocess
         if Ap_rank < p
             @warn("some primal equalities appear to be dependent " *
             "(possibly inconsistent); try using preprocess = true")
         end
-        init_y = Ap_fact \ rhs
-        return init_y
+        return
     end
 
-    # preprocess dual equalities
+    # preprocess
+    Ap_Q = Ap_fact.Q
     Ap_R = UpperTriangular(Ap_fact.R[1:Ap_rank, 1:Ap_rank])
     col_piv = (Ap_fact isa QRPivoted{T, Matrix{T}}) ? Ap_fact.p : Ap_fact.pcol
     y_keep_idxs = col_piv[1:Ap_rank]
-    Ap_Q = Ap_fact.Q
 
-    b_sub = model.b[y_keep_idxs]
     if Ap_rank < p
         # some dependent primal equalities, so check if they are consistent
         # x_sub = Ap_Q * vcat((Ap_R' \ b_sub), zeros(n - Ap_rank))
-        x_sub = zeros(T, n)
-        x_sub1 = view(x_sub, 1:Ap_rank)
-        copyto!(x_sub1, b_sub)
+        x_sub = zeros(T, model.n)
+        @views x_sub1 = x_sub[1:Ap_rank]
+        @views copyto!(x_sub1, model.b[y_keep_idxs])
         ldiv!(Ap_R', x_sub1)
         lmul!(Ap_Q, x_sub)
-
         if !(Ap_fact isa QRPivoted{T, Matrix{T}})
             x_sub = x_sub[Ap_fact.rpivinv]
         end
-        residual = norm(A * x_sub - model.b, Inf)
+        residual = norm(model.A * x_sub - model.b, Inf)
+
         if residual > solver.init_tol_qr
             if solver.verbose
                 println("some primal equality constraints are inconsistent " *
@@ -266,102 +258,129 @@ function find_initial_y(
             return zeros(T, 0)
         end
         if solver.verbose
+            p = model.p
             println("$(p - Ap_rank) of $p primal equality constraints " *
                 "are dependent")
         end
     end
 
-    if reduce && isa(model.G, MatrixyAG)
-        # remove all primal equalities by making A and b empty with
-        # n = n0 - p0 and p = 0
-        # TODO improve efficiency
-        # TODO avoid calculating GQ1 explicitly if possible
-        # recover original-space solution using:
-        # x0 = Q * [(R' \ b0), x]
-        # y0 = R \ (-cQ1' - GQ1' * z0)
-        if !(Ap_fact isa QRPivoted{T, Matrix{T}})
+    if !(solver.reduce && isa(model.G, MatrixyAG))
+        # not reducing
+        # modify model to remove/reorder some dual variables y
+        if Ap_fact isa QRPivoted{T, Matrix{T}}
+            model.A = model.A[y_keep_idxs, :]
+        else
             row_piv = Ap_fact.prow
-            model.c = model.c[row_piv]
+            model.A = model.A[y_keep_idxs, row_piv]
             model.G = model.G[:, row_piv]
-            solver.reduce_row_piv_inv = Ap_fact.rpivinv
-        else
-            solver.reduce_row_piv_inv = Int[]
+            solver.x_keep_idxs = solver.x_keep_idxs[row_piv]
         end
 
-        Q1_idxs = 1:Ap_rank
-        Q2_idxs = (Ap_rank + 1):n
-
-        # [cQ1 cQ2] = c0' * Q
-        cQ = model.c' * Ap_Q
-        cQ1 = solver.reduce_cQ1 = cQ[Q1_idxs]
-        cQ2 = cQ[Q2_idxs]
-        # c = cQ2
-        model.c = cQ2
-        model.n = length(model.c)
-        # offset = offset0 + cQ1 * (R' \ b0)
-        Rpib0 = solver.reduce_Rpib0 = ldiv!(Ap_R', b_sub)
-        # solver.Rpib0 = Rpib0 # TODO
-        model.obj_offset += dot(cQ1, Rpib0)
-
-        # [GQ1 GQ2] = G0 * Q
-        # very inefficient method used for sparse G * QRSparseQ
-        # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818
-        if model.G isa UniformScaling
-            side = size(Ap_Q, 1)
-            G_mul = Matrix{T}(model.G, side, side)
-        elseif !isa(model.G, Matrix)
-            G_mul = Matrix{T}(model.G)
-        else
-            G_mul = model.G
-        end
-        GQ = G_mul * Ap_Q
-        GQ1 = solver.reduce_GQ1 = GQ[:, Q1_idxs]
-        GQ2 = GQ[:, Q2_idxs]
-        # h = h0 - GQ1 * (R' \ b0)
-        mul!(model.h, GQ1, Rpib0, -1, true)
-
-        # G = GQ2
-        model.G = GQ2
-
-        # A and b empty
-        model.p = 0
-        model.A = zeros(T, 0, model.n)
-        model.b = zeros(T, 0)
-        solver.reduce_Ap_R = Ap_R
-        solver.reduce_Ap_Q = Ap_Q
-
-        solver.reduce_y_keep_idxs = y_keep_idxs
-        solver.Ap_R = UpperTriangular(zeros(T, 0, 0))
-        solver.Ap_Q = I
-
-        return zeros(T, 0)
+        model.p = Ap_rank
+        solver.y_keep_idxs = y_keep_idxs
+        solver.Ap_R = Ap_R
+        solver.Ap_Q = Ap_Q
+        return
     end
 
-    # init_y = Ap_R \ ((Ap_fact.Q' * (-c - G' * point.z))[1:Ap_rank])
-    temp = copy(model.c)
-    mul!(temp, model.G', init_z, true, true)
-    lmul!(Ap_fact.Q', temp)
-    init_y = temp[1:Ap_rank]
-    init_y .*= -1
-    ldiv!(Ap_R, init_y)
-
-    # modify solver.model to remove/reorder some dual variables y
-    if !(Ap_fact isa QRPivoted{T, Matrix{T}})
+    # reduce
+    if Ap_fact isa QRPivoted{T, Matrix{T}}
+        solver.reduce_row_piv_inv = Int[]
+    else
         row_piv = Ap_fact.prow
-        model.A = A[y_keep_idxs, row_piv]
-        model.c = model.c[row_piv]
         model.G = model.G[:, row_piv]
-        solver.x_keep_idxs = solver.x_keep_idxs[row_piv]
+        solver.reduce_row_piv_inv = Ap_fact.rpivinv
+    end
+
+    model.n -= Ap_rank
+    Ap_Q = Ap_fact.Q
+
+    # [GQ1 GQ2] = G0 * Q
+    # very inefficient method used for sparse G * QRSparseQ
+    # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818
+    G_mul = if model.G isa UniformScaling
+        side = size(Ap_Q, 1)
+        Matrix{T}(model.G, side, side)
+    elseif !isa(model.G, Matrix)
+        Matrix{T}(model.G)
     else
-        model.A = A[y_keep_idxs, :]
+        model.G
+    end
+    GQ = G_mul * Ap_Q
+    solver.reduce_GQ1 = GQ[:, 1:Ap_rank]
+    # G = GQ2
+    model.G = GQ[:, (Ap_rank + 1):end]
+
+    # A and b empty
+    model.p = 0
+    model.A = zeros(T, 0, model.n)
+    solver.reduce_Ap_R = Ap_R
+    solver.reduce_Ap_Q = Ap_Q
+
+    solver.reduce_y_keep_idxs = y_keep_idxs
+    solver.Ap_R = UpperTriangular(zeros(T, 0, 0))
+    solver.Ap_Q = I
+    return
+end
+
+# update data for b, c, h from primal equality preprocessing/reduction
+# and return initial y as least squares solution to A'y = -c - G'z
+function update_primal_eq(
+    solver::Solver{T},
+    init_z::Vector{T},
+    ) where {T <: Real}
+    iszero(solver.orig_model.p) && return zeros(T, 0)
+    model = solver.model
+
+    if solver.init_use_indirect
+        # use indirect method TODO pick lsqr or lsmr
+        rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1)
+        return IterativeSolvers.lsqr(model.A', rhs_y)
+    end
+
+    Ap_fact = solver.Ap_fact
+    if !solver.preprocess
+        rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1)
+        return Ap_fact \ rhs_y
+    end
+
+    # preprocess
+    Ap_rank = solver.Ap_rank
+
+    if !(solver.reduce && isa(model.G, MatrixyAG))
+        # not reducing
+        model.b = model.b[solver.y_keep_idxs]
+        if !(Ap_fact isa QRPivoted{T, Matrix{T}})
+            model.c = model.c[Ap_fact.prow]
+        end
+
+        # init_y = Ap_R \ ((Ap_Q' * (-c - G' * z))[1:Ap_rank])
+        rhs_y = mul!(copy(model.c), model.G', init_z, -1, -1)
+        lmul!(solver.Ap_Q', rhs_y)
+        init_y = rhs_y[1:Ap_rank]
+        ldiv!(solver.Ap_R, init_y)
+        return init_y
     end
-    model.b = b_sub
-    model.p = Ap_rank
-    solver.y_keep_idxs = y_keep_idxs
-    solver.Ap_R = Ap_R
-    solver.Ap_Q = Ap_Q
 
-    return init_y
+    # reduce
+    # [cQ1 cQ2] = c0' * Q
+    # c = cQ2
+    cQ = model.c' * solver.reduce_Ap_Q
+    cQ1 = solver.reduce_cQ1 = cQ[1:Ap_rank]
+    model.c = cQ[(Ap_rank + 1):end]
+
+    Rpib0 = solver.reduce_Rpib0 = model.b[solver.reduce_y_keep_idxs]
+    ldiv!(solver.reduce_Ap_R', Rpib0)
+    # offset = offset0 + cQ1 * (R' \ b0)
+    model.obj_offset += dot(cQ1, Rpib0)
+
+    model.b = zeros(T, 0)
+
+    GQ1 = solver.reduce_GQ1
+    # h = h0 - GQ1 * (R' \ b0)
+    model.h = mul!(copy(model.h), GQ1, Rpib0, -1, true)
+
+    return zeros(T, 0)
 end
 
 # (pivoted) QR factorizations are usually rank-revealing but may be unreliable
@@ -385,65 +404,46 @@ end
 function postprocess(solver::Solver{T}) where {T <: Real}
     point = solver.point
     result = solver.result
-    if solver.status in infeas_statuses
-        tau = one(T)
-    else
-        tau = point.tau[]
-        if tau <= 0
-            result.vec .= NaN
-            return
-        end
+    tau = point.tau[]
+    if tau <= 0 || any(isnan, point.vec)
+        result.vec .= NaN
+        return
     end
 
-    # finalize s,z
-    @. result.s = point.s / tau
-    @. result.z = point.z / tau
+    copyto!(result.s, point.s)
+    copyto!(result.z, point.z)
+    x = copy(point.x)
+    y = copy(point.y)
+    if !in(solver.status, infeas_statuses)
+        # rescale non-infeasible certificates by 1/tau
+        result.s ./= tau
+        result.z ./= tau
+        x ./= tau
+        y ./= tau
+    end
 
     # finalize x
-    if solver.preprocess && !iszero(solver.orig_model.n) && !any(isnan, point.x)
+    if solver.preprocess && !iszero(solver.orig_model.n)
         # unpreprocess solver's solution
         if solver.reduce && !iszero(solver.orig_model.p)
-            # unreduce solver's solution
-            # x = Q * [(R' \ b0), x]
-            xa = zeros(T, solver.orig_model.n - length(solver.reduce_Rpib0))
-            @. @views xa[solver.x_keep_idxs] = point.x / tau
-            if solver.status in infeas_statuses
-                Rpib0 = zeros(T, length(solver.reduce_Rpib0))
-            else
-                Rpib0 = solver.reduce_Rpib0
-            end
-            xb = vcat(Rpib0, xa)
-            lmul!(solver.reduce_Ap_Q, xb)
-            if isempty(solver.reduce_row_piv_inv)
-                result.x .= xb
-            else
-                @. @views result.x = xb[solver.reduce_row_piv_inv]
-            end
+            unreduce_x(solver, x)
         else
-            @. @views result.x[solver.x_keep_idxs] = point.x / tau
+            @views copyto!(result.x[solver.x_keep_idxs], x)
         end
     else
-        @. result.x = point.x / tau
+        copyto!(result.x, x)
     end
 
     # finalize y
-    if solver.preprocess && !iszero(solver.orig_model.p) && !any(isnan, point.y)
+    if solver.preprocess && !iszero(solver.orig_model.p)
         # unpreprocess solver's solution
         if solver.reduce
-            # unreduce solver's solution
-            # y = R \ (-cQ1' - GQ1' * z)
-            ya = solver.reduce_GQ1' * result.z
-            if !(solver.status in infeas_statuses)
-                ya .+= solver.reduce_cQ1
-            end
-            @views ldiv!(solver.reduce_Ap_R,
-                ya[1:length(solver.reduce_y_keep_idxs)])
-            @. @views result.y[solver.reduce_y_keep_idxs] = -ya
+            unreduce_y(solver, y)
         else
-            @. @views result.y[solver.y_keep_idxs] = point.y / tau
+            @views copyto!(result.y[solver.y_keep_idxs], y)
         end
     else
-        @. result.y = point.y / tau
+        copyto!(result.y, y)
     end
 
     if solver.used_rescaling
@@ -456,3 +456,42 @@ function postprocess(solver::Solver{T}) where {T <: Real}
 
     return
 end
+
+function unreduce_x(
+    solver::Solver{T},
+    x::Vector{T},
+    ) where {T <: Real}
+    # x = Q * [(R' \ b0), x]
+    Rpib0 = if solver.status in infeas_statuses
+        zero(solver.reduce_Rpib0)
+    else
+        solver.reduce_Rpib0
+    end
+    xa = zeros(T, solver.orig_model.n - length(Rpib0))
+    @views copyto!(xa[solver.x_keep_idxs], x)
+    xb = vcat(Rpib0, xa)
+    lmul!(solver.reduce_Ap_Q, xb)
+    result_x = solver.result.x
+    if isempty(solver.reduce_row_piv_inv)
+        copyto!(result_x, xb)
+    else
+        @views copyto!(result_x, xb[solver.reduce_row_piv_inv])
+    end
+    return
+end
+
+function unreduce_y(
+    solver::Solver{T},
+    y::Vector{T},
+    ) where {T <: Real}
+    # y = R \ (-cQ1' - GQ1' * z)
+    ya = solver.reduce_GQ1' * solver.result.z
+    if !in(solver.status, infeas_statuses)
+        ya .+= solver.reduce_cQ1
+    end
+    @views ya_sub = ya[1:length(solver.reduce_y_keep_idxs)]
+    ldiv!(solver.reduce_Ap_R, ya_sub)
+    @views y_sub = solver.result.y[solver.reduce_y_keep_idxs]
+    @. y_sub = -ya
+    return
+end
diff --git a/src/Solvers/systemsolvers/common.jl b/src/Solvers/systemsolvers/common.jl
index 7aac2cd62..2f8c39e4e 100644
--- a/src/Solvers/systemsolvers/common.jl
+++ b/src/Solvers/systemsolvers/common.jl
@@ -186,24 +186,29 @@ function setup_point_sub(
     model::Models.Model{T},
     ) where {T <: Real}
     (n, p, q) = (model.n, model.p, model.q)
-    dim_sub = n + p + q
-    z_start = model.n + model.p
-
     sol_sub = syssolver.sol_sub = Point{T}()
     rhs_sub = syssolver.rhs_sub = Point{T}()
     rhs_const = syssolver.rhs_const = Point{T}()
     sol_const = syssolver.sol_const = Point{T}()
     for point_sub in (sol_sub, rhs_sub, rhs_const, sol_const)
-        point_sub.vec = zeros(T, dim_sub)
+        point_sub.vec = zeros(T, n + p + q)
         @views point_sub.x = point_sub.vec[1:n]
         @views point_sub.y = point_sub.vec[n .+ (1:p)]
         @views point_sub.z = point_sub.vec[n + p .+ (1:q)]
         point_sub.z_views = [view(point_sub.z, idxs) for idxs in model.cone_idxs]
     end
+    set_point_sub_rhs(syssolver, model)
+    return
+end
+
+function set_point_sub_rhs(
+    syssolver::Union{QRCholSystemSolver{T}, SymIndefSystemSolver{T}},
+    model::Models.Model{T},
+    ) where {T <: Real}
+    rhs_const = syssolver.rhs_const
     @. rhs_const.x = -model.c
     @. rhs_const.y = model.b
     @. rhs_const.z = model.h
-
     return
 end
 
diff --git a/src/Solvers/systemsolvers/qrchol.jl b/src/Solvers/systemsolvers/qrchol.jl
index c13661d1e..c88c99148 100644
--- a/src/Solvers/systemsolvers/qrchol.jl
+++ b/src/Solvers/systemsolvers/qrchol.jl
@@ -151,7 +151,11 @@ function load(
 
     # very inefficient method used for sparse G * QRSparseQ
     # see https://github.com/JuliaLang/julia/issues/31124#issuecomment-501540818
-    GQ = model.G * solver.Ap_Q
+    GQ = if solver.Ap_Q == I
+        model.G
+    else
+        model.G * solver.Ap_Q
+    end
 
     syssolver.GQ2 = GQ[:, (p + 1):end]
     syssolver.HGQ2 = zeros(T, q, nmp)
diff --git a/test/nativeinstances.jl b/test/nativeinstances.jl
index 2e9e01060..7a3326d76 100644
--- a/test/nativeinstances.jl
+++ b/test/nativeinstances.jl
@@ -37,14 +37,17 @@ function build_solve_check(
     G,
     h::Vector{T},
     cones::Vector{Cone{T}},
-    tol::Real,
-    ;
+    tol::Real;
     obj_offset::T = zero(T),
     solver::Solvers.Solver{T} = Solvers.Solver{T}(),
+    already_loaded::Bool = false,
     ) where {T <: Real}
-    model = Hypatia.Models.Model{T}(c, A, b, G, h, cones, obj_offset = obj_offset)
-
-    Solvers.load(solver, model)
+    if already_loaded
+        model = Solvers.get_model(solver)
+    else
+        model = Hypatia.Models.Model{T}(c, A, b, G, h, cones, obj_offset = obj_offset)
+        Solvers.load(solver, model)
+    end
     Solvers.solve(solver)
 
     status = Solvers.get_status(solver)
@@ -2754,3 +2757,83 @@ function indirect5(T; options...)
     r = build_solve_check(c, A, b, G, h, cones, tol; options...)
     @test r.status == Solvers.PrimalInfeasible
 end
+
+function modify1(T; options...)
+    tol = test_tol(T)
+    c = T[1, 0]
+    A = T[1 1]
+    b = [-T(2)]
+    G = SparseMatrixCSC(-one(T) * I, 2, 2)
+    h = zeros(T, 2)
+    cones = Cone{T}[Cones.Nonnegative{T}(2)]
+
+    r = build_solve_check(c, A, b, G, h, cones, tol; options...)
+    @test r.status == Solvers.PrimalInfeasible
+    solver = r.solver
+
+    b = [-T(1)]
+    Solvers.modify_b(solver, b)
+    @test Solvers.get_status(solver) == Solvers.Modified
+    r = build_solve_check(c, A, b, G, h, cones, tol;
+        solver = solver, already_loaded = true, options...)
+    @test r.status == Solvers.PrimalInfeasible
+    @test iszero(solver.time_rescale)
+    @test iszero(solver.time_loadsys)
+
+    c = T[2, 1]
+    b = [T(2)]
+    obj_offset = T(1)
+    Solvers.modify_c(solver, c)
+    Solvers.modify_b(solver, b)
+    Solvers.modify_obj_offset(solver, obj_offset)
+    @test Solvers.get_status(solver) == Solvers.Modified
+    r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = obj_offset,
+        solver = solver, already_loaded = true, options...)
+    @test r.status == Solvers.Optimal
+    @test r.primal_obj ≈ 3 atol=tol rtol=tol
+    @test r.x ≈ [0, 2] atol=tol rtol=tol
+    @test iszero(solver.time_rescale)
+    @test iszero(solver.time_loadsys)
+
+    c = T[3, 1]
+    h = T[1, -2]
+    obj_offset = T(-1)
+    Solvers.modify_c(solver, c)
+    Solvers.modify_h(solver, h)
+    Solvers.modify_obj_offset(solver, obj_offset)
+    @test Solvers.get_status(solver) == Solvers.Modified
+    r = build_solve_check(c, A, b, G, h, cones, tol; obj_offset = obj_offset,
+        solver = solver, already_loaded = true, options...)
+    @test r.status == Solvers.Optimal
+    @test r.primal_obj ≈ -1 atol=tol rtol=tol
+    @test r.x ≈ [-1, 3] atol=tol rtol=tol
+    @test iszero(solver.time_rescale)
+    @test iszero(solver.time_loadsys)
+end
+
+function modify2(T; options...)
+    tol = test_tol(T)
+    c = T[-1, -1, 0]
+    A = zeros(T, 0, 3)
+    b = T[]
+    G = repeat(SparseMatrixCSC(-one(T) * I, 3, 3), 2, 1)
+    h = zeros(T, 6)
+    cones = Cone{T}[Cones.EpiNormInf{T, T}(3),
+        Cones.EpiNormInf{T, T}(3, use_dual = true)]
+
+    r = build_solve_check(c, A, b, G, h, cones, tol; options...)
+    @test r.status == Solvers.DualInfeasible
+    solver = r.solver
+
+    c = zeros(T, 3)
+    Solvers.modify_c(solver, c)
+    @test Solvers.get_status(solver) == Solvers.Modified
+    r = build_solve_check(c, A, b, G, h, cones, tol;
+        solver = solver, already_loaded = true, options...)
+    @test r.status == Solvers.Optimal
+    @test r.primal_obj ≈ 0 atol=tol rtol=tol
+    @test r.dual_obj ≈ 0 atol=tol rtol=tol
+    @test norm(r.z) ≈ 0 atol=tol rtol=tol
+    @test iszero(solver.time_rescale)
+    @test iszero(solver.time_loadsys)
+end
\ No newline at end of file
diff --git a/test/nativesets.jl b/test/nativesets.jl
index 3ecd15131..9de6c5a6e 100644
--- a/test/nativesets.jl
+++ b/test/nativesets.jl
@@ -177,3 +177,8 @@ inst_indirect = [
     "indirect4",
     "indirect5",
     ]
+
+inst_modify = [
+    "modify1",
+    "modify2",
+    ]
diff --git a/test/runmoitests.jl b/test/runmoitests.jl
index 8c1f5c58d..77b26e21c 100644
--- a/test/runmoitests.jl
+++ b/test/runmoitests.jl
@@ -29,7 +29,7 @@ end
 test_T = [
     (Float64, 2 * sqrt(sqrt(eps())), 4, String[]),
     # TODO add test_linear after MOI 0.10.7 is tagged:
-    (BigFloat, 2 * eps(BigFloat)^0.2, 1, String["test_conic"]),
+    (BigFloat, 2 * eps(BigFloat)^0.2, 1, String["test_conic", "test_modification"]),
 ]
 
 @testset "MOI.Test tests: $T" for (T, tol_test, tol_relax, includes) in test_T
diff --git a/test/runnativetests.jl b/test/runnativetests.jl
index 2eb2853ce..36b18b71f 100644
--- a/test/runnativetests.jl
+++ b/test/runnativetests.jl
@@ -61,8 +61,8 @@ end
 @testset "no preprocess tests" begin
     println("\nstarting no preprocess tests")
     for inst_name in inst_cones_few, T in diff_reals
-        options = (; default_options..., preprocess = false, reduce = false,
-            syssolver = Solvers.SymIndefDenseSystemSolver{T}())
+        options = (; default_options..., rescale = false, preprocess = false,
+            reduce = false, syssolver = Solvers.SymIndefDenseSystemSolver{T}())
         test_instance_solver(inst_name, T, options)
     end
 end
@@ -137,5 +137,25 @@ end
     end
 end
 
+@testset "model modification tests" begin
+    println("\nstarting model modification tests")
+    syssolvers = [
+        (Solvers.SymIndefDenseSystemSolver, diff_reals),
+        (Solvers.SymIndefSparseSystemSolver, [Float64,]),
+        (Solvers.QRCholDenseSystemSolver, diff_reals),
+        ]
+    for inst_name in inst_modify, (syssolver, real_types) in syssolvers,
+        T in real_types, rescale in (false, true), preprocess in (false, true),
+        reduce in (false, true)
+        if !preprocess
+            reduce && continue
+            syssolver == Solvers.QRCholDenseSystemSolver && continue
+        end
+        options = (; default_options..., syssolver = syssolver{T}(),
+            rescale = rescale, preprocess = preprocess, reduce = reduce)
+        test_instance_solver(inst_name, T, options)
+    end
+end
+
 end
 ;