Struct Solver 

pub struct Solver { /* private fields */ }

Session-style solver: holds an IpoptApplication, its TNLP, and the converged factor between calls.

Implementations

impl Solver

pub fn new(app: IpoptApplication, tnlp: Rc<RefCell<dyn TNLP>>) -> Self

Build a new session. The app should already have its options configured and initialize() called.

Examples found in repository

examples/sensitivity_session.rs (line 166)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 178)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 167)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

pub fn app(&self) -> &IpoptApplication

Borrow the underlying IpoptApplication (e.g. to read its options table after a solve). Mutation between solve calls is supported via Self::app_mut.

pub fn app_mut(&mut self) -> &mut IpoptApplication

Mutable borrow of the underlying IpoptApplication. Useful for reconfiguring options before a follow-up solve(). Note that changing options that affect the KKT linear system between calls will invalidate the cached factor; the next solve() rebuilds it.

pub fn solve(&mut self) -> ApplicationReturnStatus

Run the IPM to convergence. On a successful solve the ConvergedState (including the KKT backsolver) is stashed inside the Solver and accessible via Self::converged.

Each call to solve() overwrites the previous converged state; the previously held factor is dropped.

Examples found in repository

examples/sensitivity_session.rs (line 167)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 180)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 168)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

pub fn converged(&self) -> Option<Ref<'_, ConvergedState>>

Borrow the converged state, if a successful solve has been run. Returns None if no solve has run or if the most recent solve failed before reaching convergence.

Examples found in repository

examples/sensitivity_session.rs (line 169)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 186)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

pub fn kkt_dim(&self) -> Option<usize>

Total dimension of the compound KKT vector (sum of block_dims). Returns None if no converged factor is held.

Examples found in repository

examples/sensitivity_session.rs (line 188)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

pub fn block_dims(&self) -> Option<[usize; 8]>

Block dimensions of the compound KKT vector in (x, s, y_c, y_d, z_l, z_u, v_l, v_u) order. Returns None if no converged factor is held.

pub fn kkt_solve( &self, rhs: &[Number], lhs: &mut [Number], ) -> Result<(), SolverError>

Solve K · lhs = rhs against the converged KKT factor. Both slices must have length kkt_dim(); the layout is the flat x || s || y_c || y_d || z_l || z_u || v_l || v_u packing.

Examples found in repository

examples/sensitivity_session.rs (line 191)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

pub fn kkt_solve_many( &self, rhs_flat: &[Number], lhs_flat: &mut [Number], n_rhs: usize, ) -> Result<(), SolverError>

Batched-RHS back-solve. rhs_flat and lhs_flat are row-major (n_rhs, kkt_dim) buffers; each row is solved against the same converged factor. Equivalent in result to looping Self::kkt_solve but reuses one IteratesVector for the RHS and one for the result across all n_rhs calls — see crate::algorithm_backsolver::PdSensBacksolver::solve_many.

pub fn parametric_step( &self, pin_constraint_indices: &[Index], deltas: &[Number], ) -> Result<Vec<Number>, SolverError>

First-order parametric step Δx ≈ ∂x*/∂p · Δp for a set of pinned equality constraints. pin_constraint_indices are 0-based indices into the user’s g(x); deltas is the perturbation Δp (same length).

Returns the n_x-long primal step. For the full KKT-space step, use Self::kkt_solve directly.

Examples found in repository

examples/sensitivity_session.rs (line 176)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 197)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 174)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

pub fn compute_reduced_hessian( &self, pin_constraint_indices: &[Index], obj_scal: Number, ) -> Result<Vec<Number>, SolverError>

Reduced Hessian H_R = obj_scal · B K⁻¹ Bᵀ over the pinned equality-constraint rows, where B selects the pin_constraint_indices rows of the y_c block. Returns the n²-long column-major dense matrix (n = pin_constraint_indices.len()).

Equivalent to crate::SensSolve::with_reduced_hessian but usable post-hoc on a held Solver.

Examples found in repository

examples/sensitivity_session.rs (line 183)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}