Function jacobian_batch 

pub fn jacobian_batch<F>(
    points: &[f64],
    n: usize,
    settings: DiffSettings,
    backend: ComputeBackend,
    f: F,
) -> JacobianBatch
where
    F: Fn(usize, &[f64], &mut Vec<f64>) + Sync,

Assemble N independent Jacobians, one per lane, on the chosen backend.

This is the batched form: the parallel axis is the lane, not the column, so it is the right entry point when there are many small Jacobians (a per-cell chemistry Jacobian over a mesh, an ensemble of ODE systems). Use jacobian when there is one large Jacobian instead.

The flat point layout

points is a flat, row-major buffer of lanes * n values: lane i’s point is points[i * n .. (i + 1) * n]. A &[Vec<f64>] would be the obvious alternative and is rejected on purpose — it costs one allocation and one pointer chase per lane, and it is not the layout a GPU buffer would ever take. points.len() must be an exact multiple of n.

Arguments

• points — flat lanes * n buffer as above, in the caller’s own units.
• n — the dimension of each point, dimensionless. Must be non-zero.
• settings — scheme and step-size policy, shared by every lane.
• backend — requested backend; see jacobian_batch_backend_for. Each lane’s columns are computed serially, since the lane axis is already the parallel one.
• f — f(i, x, out) must fill out with lane i’s n function components at x. It must be a pure deterministic function of its arguments.

Returns

A JacobianBatch with one solution per lane, in lane order. An empty points slice, or n == 0, returns an empty batch and calls f zero times. A points.len() that is not a multiple of n returns an empty batch as well — it is a caller bug, not a numerical failure, and there is no sensible lane count to report per-lane statuses against.

Cost

lanes * DiffScheme::evaluations_per_jacobian(n) calls to f.

Example

use outram_foam_basic_lib::compute::ComputeBackend;
use outram_foam_basic_lib::math::differentiate::{jacobian_batch, DiffSettings};

// 3 lanes of the 2-D rotation-like system f = [-k_i * y, k_i * x],
// whose Jacobian is [[0, -k_i], [k_i, 0]].
let k = [1.0_f64, 2.5, 7.0];
let points: Vec<f64> = vec![0.3, -0.7, 0.3, -0.7, 0.3, -0.7];

let batch = jacobian_batch(
    &points,
    2,
    DiffSettings::central(),
    ComputeBackend::CpuMulti,
    |i, v: &[f64], out: &mut Vec<f64>| {
        out.push(-k[i] * v[1]);
        out.push(k[i] * v[0]);
    },
);

let mats = batch.matrices().expect("linear system, exact everywhere");
for (i, m) in mats.iter().enumerate() {
    assert!((m.get(0, 1) + k[i]).abs() < 1e-9);
    assert!((m.get(1, 0) - k[i]).abs() < 1e-9);
}