gpuGEM wraps NVIDIA cuOpt PDLP solver with validated default settings derived from systematic benchmarking on whole-body metabolic models, including the Harvey model and personalised microbiome whole-body models. It provides a minimal, COBRA-compatible API so existing modelling workflows can switch to GPU solving with one line of code.

Why gpuGEM?

GPU first-order LP solvers (PDLP) offer substantial speedups on genome-scale FBA, but getting correct results requires non-obvious configuration. Out-of-the-box cuOpt settings produce large constraint violations or silently return wrong solutions on whole-body models. gpuGEM ships with settings that have been validated for correctness:

Installation

pip install gpugem

This installs the full tool: numpy, scipy, NVIDIA cuOpt (the GPU solver) and COBRApy. gpuGEM is a GPU solver, so solving needs an NVIDIA GPU with CUDA 12 and Python 3.11 or later. cuOpt is a large download (about 1.4 GB, mostly libcudf-cu12 and libcuopt-cu12). The cuopt-cu12 entry on PyPI is a stub that fetches the real wheels from NVIDIA's index. Behind a proxy that blocks pypi.nvidia.com it can fail with a confusing build error; if that happens, name the index explicitly:

pip install gpugem --extra-index-url https://pypi.nvidia.com

import gpugem works without a usable GPU; the GPU is checked when you call solve, which raises a DependencyError explaining what is missing.

Optional extras:

pip install "gpugem[gurobi]"       # the Gurobi CPU baseline (needs a licence)
pip install "gpugem[application]"  # the candidate-gene application, see application/
pip install "gpugem[dev]"          # tests and linting

Dependencies

Troubleshooting

A missing or unusable dependency raises gpugem.DependencyError (benchmark and example scripts print the same message; benchmark scripts exit with status 3). The message always says what is wrong, why gpuGEM needs it, and how to fix it. The original exception is kept as __cause__. err.kind is one of:

try:
    result = gpugem.solve_cobra(model)
except gpugem.DependencyError as err:
    print(err.kind, err.dependency, err.remedy)

Quick start

From a COBRApy model:

import cobra
import gpugem

model = cobra.io.read_sbml_model("e_coli_core.xml")
result = gpugem.solve_cobra(model)

print(result)                       # FBAResult(status='Optimal', objective=0.8739, time=0.42s)
print(result.objective)
print(result.fluxes)                # numpy array, same order as model.reactions

From a COBRA .mat file (MATLAB format):

import gpugem

arrays = gpugem.loaders.from_mat("Harvey_1_03c_reduced.mat")
result = gpugem.solve(**arrays)

From raw matrices (scipy sparse):

import gpugem

result = gpugem.solve(
    S=S,           # stoichiometric matrix (sparse)
    b=b,           # RHS equalities
    lb=lb, ub=ub,  # variable bounds
    c=c,           # objective coefficients
    maximize=True,
)

Whole-body models with coupling constraints:

result = gpugem.solve(S=S, b=b, lb=lb, ub=ub, c=c,
                      C=C, d_lb=d_lb, d_ub=d_ub,
                      maximize=True)

Default settings

gpuGEM automatically selects solver settings based on model size:

Small models (≤ 100K reactions) — e.g. Harvey whole-body model:

pdlp_precision             = 1      # mixed FP32/FP64 — double precision diverges on these models
absolute_primal_tolerance  = 1e-8   # tight tolerance: 0 stoich violations at 1e-6 max residual
relative_primal_tolerance  = 1e-8
absolute_dual_tolerance    = 1e-8
relative_dual_tolerance    = 1e-8
per_constraint_residual    = 1      # max-norm convergence check → better accuracy

Large models (> 100K reactions) — e.g. personalised microbiome whole-body models:

presolve                   = 1      # PaPILO — required; default PSLP falsely reports Infeasible
per_constraint_residual    = 1      # max-norm convergence → 44% fewer violations vs default
absolute_primal_tolerance  = 1e-4   # cuOpt's internal default tolerance
relative_primal_tolerance  = 1e-4
absolute_dual_tolerance    = 1e-4
relative_dual_tolerance    = 1e-4

Any cuOpt parameter can be overridden:

# Use Stable1 PDLP mode for maximum accuracy (slower)
result = gpugem.solve(**arrays, pdlp_solver_mode=0, time_limit=300)

# Disable presolve explicitly
result = gpugem.solve(**arrays, presolve=0)

Time limit for large models

time_limit is the wall-clock budget in seconds and defaults to 60. That is enough for genome-scale models such as iML1515 (about 4 s) and the Harvey whole-body model (a few seconds), but not for the million-variable personalised microbiome models, which took roughly 25 to 880 s on an NVIDIA RTX A4500. For those, pass a larger limit, for example time_limit=900. If the limit is reached, result.status is "TimeLimit" and the solution has not converged (its residual can be far above the 1e-4 criterion), so check result.status and result.feasibility["stoich_max_residual"] before using the fluxes.

Stateful solver (repeated solves)

solver = gpugem.FBASolver(S=S, b=b, lb=lb, ub=ub, c=c, maximize=True)

for patient_c in patient_objectives:
    result = solver.solve(c=patient_c)
    print(result.objective)

Result object

result.status       # 'Optimal', 'TimeLimit', 'Infeasible', ...
result.objective    # float
result.fluxes       # np.ndarray (n_vars,)
result.wall_time_s  # float
result.solver_settings  # dict of applied cuOpt parameters
result.feasibility  # dict: stoich_max_residual, stoich_rows_violated_1e6, ...

Known limitations

• Lifted, extreme-scale (S84/S85-class) models are not reliably solvable by every Gurobi entry point. benchmarks/run_matlab_python_lifting_comparison.py found that Gurobi solves the lifted S84/S85 systems to Optimal via COBRA Toolbox's solveCobraLP (MATLAB) but reports Infeasible for the structurally identical system (confirmed: identical auxiliary-variable and lifted-row counts) via a direct gurobipy call (Python). Checked and ruled out as causes: a Gurobi solver-settings mismatch (Method, TimeLimit, FeasibilityTol/OptimalityTol all confirmed identical), a Gurobi version mismatch (pip install gurobipy initially resolved a different major version, 13.0.3, than MATLAB's licensed 11.0.3; re-run with the versions pinned to match — Infeasible persisted), and a borderline numerical-tolerance flip (persisted after relaxing FeasibilityTol 100x). This reproduces, under a second independent solver pathway, the same S84/S85 scale boundary already found for lifted cuOpt solves (see "Version x lifting runtime comparison" above) — lifting's mathematical correctness is not in question (structural fidelity is exact), but its practical solvability at this scale is entry-point-dependent, for a reason not yet isolated. See specs/015-matlab-python-lifting-comparison/.

• GPU/driver-absent detection is a heuristic, checked on one kind of machine. gpuGEM reports DependencyError(kind="no_gpu") only when cuOpt raises and there is positive evidence of a GPU problem (nvidia-smi fails or lists no GPU, or it is absent and no /dev/nvidia* device exists). Observed with cuOpt 26.8.0 on a machine with no NVIDIA driver and no nvidia-smi: import cuopt succeeds and the first solve raises RuntimeError: CUDA error ... cudaErrorInsufficientDriver, which gpuGEM reports as no_gpu. A machine with a working nvidia-smi but a broken CUDA setup has not been tested, and its failure may surface as cuOpt's own exception. The Gurobi licence error numbers (10009, 10010) used for kind="license" are unverified. See specs/016-missing-dependency-errors/.

• The dependency inventory is static. tests/test_dependency_inventory.py finds imports and subprocess / shutil.which calls with literal names; a module or executable whose name is built at run time is not seen (the test prints those it cannot analyse).

Running the tests

pip install -e ".[dev]"
pytest -q -rs

No GPU is required to run the suite. On a CPU-only machine the tests that need cuOpt or a Gurobi licence are reported as skipped rather than failed, so a clean run looks like:

206 passed, 10 skipped

If you see failures rather than skips, that is a genuine problem and we would like to hear about it. The 10 skips become real tests on a machine with an NVIDIA GPU (CUDA 12) and pip install gpugem.

CI runs this suite on Python 3.11 and 3.12 on every push. Because GitHub's hosted runners have no GPU, a green badge means "everything checkable without a GPU is correct"; it does not assert that the GPU solving path was exercised. That path is covered by the benchmark runs reported in the paper.

Support and maintenance

gpuGEM is maintained by the Digital Metabolic Twin Centre at the University of Galway. We commit to keeping the software available and to responding to issues for at least two years following publication. Please report problems through the issue tracker.

Every released version is archived on Zenodo with a permanent DOI, so the exact code behind a published result stays retrievable even if this repository moves.

Citation

If you use gpuGEM in your research, please cite:

gpuGEM: validated GPU solving of genome-scale metabolic LPs
Digital Metabolic Twin Centre, 2026. https://github.com/Digital-Metabolic-Twin-Centre/gpuGEM

License

MIT © 2026 Digital Metabolic Twin Centre