diff --git a/.github/workflows/freethreading_tests.yml b/.github/workflows/freethreading_tests.yml
new file mode 100644
index 00000000..8d3126e7
--- /dev/null
+++ b/.github/workflows/freethreading_tests.yml
@@ -0,0 +1,128 @@
+name: Free-threading tests
+
+on:
+  pull_request:
+  push:
+    branches:
+      - master
+
+jobs:
+  test_freethreading:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.13t", "3.14t"]
+
+    steps:
+      - uses: actions/checkout@v6
+        with:
+          submodules: recursive
+
+      - name: Install uv
+        uses: astral-sh/setup-uv@v6
+
+      - name: Set up free-threaded Python ${{ matrix.python-version }}
+        run: uv python install ${{ matrix.python-version }}
+
+      - name: Install OpenBLAS
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y libopenblas-dev
+
+      - name: Create venv and install dependencies
+        run: |
+          uv venv --python ${{ matrix.python-version }} .venv
+          source .venv/bin/activate
+          uv pip install numpy scipy pytest pytest-run-parallel
+
+      - name: Build and install scs
+        run: |
+          source .venv/bin/activate
+          uv pip install -v .
+
+      - name: Run tests
+        run: |
+          source .venv/bin/activate
+          python -m pytest test/ -v
+
+      - name: Run tests in parallel threads (race detection)
+        run: |
+          source .venv/bin/activate
+          python -m pytest test/ -v --parallel-threads=4 --iterations=3
+
+  test_sanitizers:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - sanitizer: tsan
+            cpython_configure_flag: "--with-thread-sanitizer"
+            cflags: "-fsanitize=thread -g"
+          - sanitizer: asan
+            cpython_configure_flag: "--with-address-sanitizer"
+            cflags: "-fsanitize=address -g"
+
+    steps:
+      - uses: actions/checkout@v6
+        with:
+          submodules: recursive
+
+      - name: Install build dependencies
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y clang libopenblas-dev libssl-dev zlib1g-dev \
+            libbz2-dev libreadline-dev libsqlite3-dev libncurses5-dev \
+            libncursesw5-dev xz-utils libffi-dev liblzma-dev
+
+      - name: Cache CPython ${{ matrix.sanitizer }} build
+        id: cache-cpython
+        uses: actions/cache@v4
+        with:
+          path: cpython-${{ matrix.sanitizer }}
+          key: cpython-${{ matrix.sanitizer }}-3.14-${{ runner.os }}-v1
+
+      - name: Build CPython 3.14t with ${{ matrix.sanitizer }}
+        if: steps.cache-cpython.outputs.cache-hit != 'true'
+        run: |
+          git clone --depth 1 https://github.com/python/cpython.git -b 3.14 cpython-src
+          cd cpython-src
+          CC=clang CXX=clang++ ./configure --disable-gil ${{ matrix.cpython_configure_flag }} \
+            --prefix $GITHUB_WORKSPACE/cpython-${{ matrix.sanitizer }}
+          make -j$(nproc)
+          make install
+
+      - name: Create venv and install dependencies
+        run: |
+          $GITHUB_WORKSPACE/cpython-${{ matrix.sanitizer }}/bin/python3.14t -m venv .san-venv
+          source .san-venv/bin/activate
+          pip install numpy scipy pytest meson-python meson ninja
+
+      - name: Set sanitizer runtime options
+        run: |
+          if [ "${{ matrix.sanitizer }}" = "tsan" ]; then
+            echo "TSAN_OPTIONS=halt_on_error=1 allocator_may_return_null=1 suppressions=${{ github.workspace }}/test/tsan-suppressions.txt" >> $GITHUB_ENV
+            echo "OPENBLAS_NUM_THREADS=1" >> $GITHUB_ENV
+          elif [ "${{ matrix.sanitizer }}" = "asan" ]; then
+            echo "ASAN_OPTIONS=halt_on_error=1 allocator_may_return_null=1" >> $GITHUB_ENV
+            echo "LSAN_OPTIONS=suppressions=${{ github.workspace }}/test/lsan-suppressions.txt" >> $GITHUB_ENV
+          fi
+
+      - name: Build and install scs
+        run: |
+          source .san-venv/bin/activate
+          pip install -v . --no-build-isolation
+        env:
+          CC: clang
+          CFLAGS: ${{ matrix.cflags }}
+
+      - name: Run tests under ${{ matrix.sanitizer }}
+        run: |
+          source .san-venv/bin/activate
+          python -m pytest test/ -v -s
+
+      - name: Run threading stress tests under ${{ matrix.sanitizer }}
+        run: |
+          source .san-venv/bin/activate
+          python -m pytest test/test_free_threading.py test/test_thread_safety.py -v -s
diff --git a/meson.build b/meson.build
index c08f6d09..ae80eb57 100644
--- a/meson.build
+++ b/meson.build
@@ -132,6 +132,7 @@ if get_option('native_arch')
 endif
 
 _deps = [blas_deps]
+_deps += dependency('threads')
 if get_option('use_openmp')
   _deps += dependency('openmp')
 endif
diff --git a/pyproject.toml b/pyproject.toml
index 0a803ae0..dd719a0f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,6 +23,7 @@ classifiers = [
     'Programming Language :: Python :: 3.13',
     'Programming Language :: Python :: 3.14',
     'Programming Language :: Python :: 3 :: Only',
+    'Programming Language :: Python :: Free Threading :: 3 - Stable',
     'Programming Language :: Python :: Implementation :: CPython',
     'Operating System :: Microsoft :: Windows',
     'Operating System :: POSIX',
@@ -96,3 +97,11 @@ before-all = [
 testpaths = [
     "test",
 ]
+faulthandler_timeout = 600
+faulthandler_exit_on_timeout = true
+markers = [
+    "thread_unsafe: mark test as unsafe to run in parallel threads (pytest-run-parallel)",
+]
+# pytest-run-parallel: functions that are not thread-safe
+# (capsys/capfd are detected automatically by pytest-run-parallel)
+thread_unsafe_functions = ["numpy.random.seed"]
diff --git a/scs/py/__init__.py b/scs/py/__init__.py
index 98cb47e6..d5387640 100644
--- a/scs/py/__init__.py
+++ b/scs/py/__init__.py
@@ -95,6 +95,11 @@ def __init__(self, data, cone, **settings):
     @param data     Dictionary containing keys `P`, `A`, `b`, `c`.
     @param cone     Dictionary containing cone information.
     @param settings Settings as kwargs, see docs.
+
+    Thread safety: construction is assumed to be thread-local. Calling
+    `__init__` on a live SCS instance from another thread (i.e. while
+    `solve` or `update` may be running on it) is undefined behavior.
+    Use a fresh `SCS(...)` instance instead.
     """
     self._settings = settings
     if not data or not cone:
diff --git a/scs/scsmodule.h b/scs/scsmodule.h
index 3673b5c5..6b9ae36d 100644
--- a/scs/scsmodule.h
+++ b/scs/scsmodule.h
@@ -64,6 +64,10 @@ static PyObject *moduleinit(void) {
     return NULL;
   }
 
+#ifdef Py_GIL_DISABLED
+  PyUnstable_Module_SetGIL(m, Py_MOD_GIL_NOT_USED);
+#endif
+
   /* Initialize SCS_Type */
   SCS_Type.tp_new = PyType_GenericNew;
   if (PyType_Ready(&SCS_Type) < 0)
diff --git a/scs/scsobject.h b/scs/scsobject.h
index d0ac98c7..a687eed9 100644
--- a/scs/scsobject.h
+++ b/scs/scsobject.h
@@ -5,7 +5,8 @@
 
 /* SCS Object type */
 typedef struct {
-  PyObject_HEAD ScsWork *work; /* Workspace */
+  PyObject_HEAD
+  ScsWork *work; /* Workspace */
   ScsSolution *sol;            /* Solution, keep around for warm-starts */
   scs_int m, n;
   PyThread_type_lock lock;     /* Per-instance lock protecting work/sol */
@@ -88,11 +89,17 @@ static int get_pos_int_param(char *key, scs_int *v, scs_int defVal,
                              PyObject *opts) {
   *v = defVal;
   if (opts) {
-    PyObject *obj = PyDict_GetItemString(opts, key);
+    PyObject *obj = NULL;
+    int rc = PyDict_GetItemStringRef(opts, key, &obj);
+    if (rc < 0) {
+      return printErr(key);
+    }
     if (obj) {
       if (parse_pos_scs_int(obj, v) < 0) {
+        Py_DECREF(obj);
         return printErr(key);
       }
+      Py_DECREF(obj);
     }
   }
   return 0;
@@ -116,26 +123,37 @@ static int get_cone_arr_dim(char *key, scs_int **varr, scs_int *vsize,
   /* get cone['key'] */
   scs_int i, n = 0;
   scs_int *q = NULL;
-  PyObject *obj = PyDict_GetItemString(cone, key);
+  PyObject *obj = NULL;
+  int rc = PyDict_GetItemStringRef(cone, key, &obj);
+  if (rc < 0) {
+    return printErr(key);
+  }
   if (obj) {
     if (PyList_Check(obj)) {
       n = (scs_int)PyList_Size(obj);
       q = (scs_int *)scs_calloc(n, sizeof(scs_int));
       for (i = 0; i < n; ++i) {
-        PyObject *qi = PyList_GetItem(obj, i);
-        if (parse_pos_scs_int(qi, &(q[i])) < 0) {
+        PyObject *qi = PyList_GetItemRef(obj, i);
+        if (!qi || parse_pos_scs_int(qi, &(q[i])) < 0) {
+          Py_XDECREF(qi);
+          scs_free(q);
+          Py_DECREF(obj);
           return printErr(key);
         }
+        Py_DECREF(qi);
       }
     } else if (PyInt_Check(obj) || PyLong_Check(obj)) {
       n = 1;
       q = (scs_int *)scs_malloc(sizeof(scs_int));
       if (parse_pos_scs_int(obj, q) < 0) {
+        scs_free(q);
+        Py_DECREF(obj);
         return printErr(key);
       }
     } else if (PyArray_Check(obj)) {
       PyArrayObject *pobj = (PyArrayObject *)obj;
       if (!PyArray_ISINTEGER(pobj) || PyArray_NDIM(pobj) != 1) {
+        Py_DECREF(obj);
         return printErr(key);
       }
       n = (scs_int)PyArray_Size((PyObject *)obj);
@@ -144,12 +162,16 @@ static int get_cone_arr_dim(char *key, scs_int **varr, scs_int *vsize,
       memcpy(q, (scs_int *)PyArray_DATA(px0), n * sizeof(scs_int));
       Py_DECREF(px0);
     } else {
+      Py_DECREF(obj);
       return printErr(key);
     }
     if (PyErr_Occurred()) {
       /* potentially could have been triggered before */
+      scs_free(q);
+      Py_DECREF(obj);
       return printErr(key);
     }
+    Py_DECREF(obj);
   }
   *vsize = n;
   *varr = q;
@@ -161,14 +183,24 @@ static int get_cone_float_arr(char *key, scs_float **varr, scs_int *vsize,
   /* get cone['key'] */
   scs_int i, n = 0;
   scs_float *q = NULL;
-  PyObject *obj = PyDict_GetItemString(cone, key);
+  PyObject *obj = NULL;
+  int rc = PyDict_GetItemStringRef(cone, key, &obj);
+  if (rc < 0) {
+    return printErr(key);
+  }
   if (obj) {
     if (PyList_Check(obj)) {
       n = (scs_int)PyList_Size(obj);
       q = (scs_float *)scs_calloc(n, sizeof(scs_float));
       for (i = 0; i < n; ++i) {
-        PyObject *qi = PyList_GetItem(obj, i);
+        PyObject *qi = PyList_GetItemRef(obj, i);
+        if (!qi) {
+          scs_free(q);
+          Py_DECREF(obj);
+          return printErr(key);
+        }
         q[i] = (scs_float)PyFloat_AsDouble(qi);
+        Py_DECREF(qi);
       }
     } else if (PyInt_Check(obj) || PyLong_Check(obj) || PyFloat_Check(obj)) {
       n = 1;
@@ -177,6 +209,7 @@ static int get_cone_float_arr(char *key, scs_float **varr, scs_int *vsize,
     } else if (PyArray_Check(obj)) {
       PyArrayObject *pobj = (PyArrayObject *)obj;
       if (!PyArray_ISFLOAT(pobj) || PyArray_NDIM(pobj) != 1) {
+        Py_DECREF(obj);
         return printErr(key);
       }
       n = (scs_int)PyArray_Size((PyObject *)obj);
@@ -185,12 +218,16 @@ static int get_cone_float_arr(char *key, scs_float **varr, scs_int *vsize,
       memcpy(q, (scs_float *)PyArray_DATA(px0), n * sizeof(scs_float));
       Py_DECREF(px0);
     } else {
+      Py_DECREF(obj);
       return printErr(key);
     }
     if (PyErr_Occurred()) {
       /* potentially could have been triggered before */
+      scs_free(q);
+      Py_DECREF(obj);
       return printErr(key);
     }
+    Py_DECREF(obj);
   }
   *vsize = n;
   *varr = q;
@@ -674,10 +711,6 @@ static PyObject *SCS_solve(SCS *self, PyObject *args) {
   npy_intp veclen[1];
   int scs_float_type = scs_get_float_type();
 
-  if (!self->work) {
-    return none_with_error("Workspace not initialized!");
-  }
-
   PyArrayObject *warm_x, *warm_y, *warm_s;
   PyObject *warm_start;
 
@@ -697,10 +730,21 @@ static PyObject *SCS_solve(SCS *self, PyObject *args) {
   /* Acquire per-instance lock. Release the GIL first to avoid deadlock:
    * another thread may hold this lock inside scs_solve (with GIL released),
    * so we must not hold the GIL while waiting for the lock. */
+  int lock_ok;
   Py_BEGIN_ALLOW_THREADS;
-  PyThread_acquire_lock(self->lock, WAIT_LOCK);
+  lock_ok = (PyThread_acquire_lock(self->lock, WAIT_LOCK) == PY_LOCK_ACQUIRED);
   Py_END_ALLOW_THREADS;
 
+  if (!lock_ok) {
+    return none_with_error("Failed to acquire instance lock");
+  }
+
+  /* Check workspace under lock to avoid TOCTOU race with SCS_finish */
+  if (!self->work) {
+    PyThread_release_lock(self->lock);
+    return none_with_error("Workspace not initialized!");
+  }
+
   if (_warm_start) {
     /* If any of these of missing, we use the values in sol */
     if (!Py_IsNone((PyObject *)warm_x)) {
@@ -734,7 +778,9 @@ static PyObject *SCS_solve(SCS *self, PyObject *args) {
   Py_END_ALLOW_THREADS;
 
   /* Copy results out of sol while still holding the lock, because another
-   * thread's solve could overwrite sol as soon as we release. */
+   * thread's solve could overwrite sol as soon as we release.
+   * Note: unlike SCS_update, we release the lock after Py_END_ALLOW_THREADS
+   * because we need to read from sol (shared state) under lock protection. */
   veclen[0] = self->n;
   _x = scs_malloc(self->n * sizeof(scs_float));
   memcpy(_x, sol->x, self->n * sizeof(scs_float));
@@ -828,11 +874,6 @@ PyObject *SCS_update(SCS *self, PyObject *args) {
   PyArrayObject *b_new, *c_new;
   scs_float *b = NULL, *c = NULL;
 
-  /* Check that the workspace is already initialized */
-  if (!self->work) {
-    return none_with_error("Workspace not initialized!");
-  }
-
   /* b, c can be None, so don't check is PyArray_Type */
   if (!PyArg_ParseTuple(args, "OO", &b_new, &c_new)) {
     return none_with_error("Error parsing inputs");
@@ -863,8 +904,32 @@ PyObject *SCS_update(SCS *self, PyObject *args) {
   }
 
   /* Acquire per-instance lock (release GIL first to avoid deadlock) */
+  int lock_ok;
+  Py_BEGIN_ALLOW_THREADS;
+  lock_ok = (PyThread_acquire_lock(self->lock, WAIT_LOCK) == PY_LOCK_ACQUIRED);
+  Py_END_ALLOW_THREADS;
+
+  if (!lock_ok) {
+    if (b) { Py_DECREF(b_new); }
+    if (c) { Py_DECREF(c_new); }
+    return none_with_error("Failed to acquire instance lock");
+  }
+
+  /* Check workspace under lock to avoid TOCTOU race with SCS_finish */
+  if (!self->work) {
+    PyThread_release_lock(self->lock);
+    if (b) { Py_DECREF(b_new); }
+    if (c) { Py_DECREF(c_new); }
+    return none_with_error("Workspace not initialized!");
+  }
+
+  /* Release the GIL, run the update, then release the instance lock before
+   * re-acquiring the GIL. This avoids holding the instance lock while
+   * waiting for the GIL, which could deadlock if another GIL-holding thread
+   * is waiting on this lock. SCS_solve uses a different order (release lock
+   * after Py_END_ALLOW_THREADS) because it must copy results out of sol
+   * while still holding the lock. */
   Py_BEGIN_ALLOW_THREADS;
-  PyThread_acquire_lock(self->lock, WAIT_LOCK);
   scs_update(self->work, b, c);
   PyThread_release_lock(self->lock);
   Py_END_ALLOW_THREADS;
@@ -883,7 +948,10 @@ PyObject *SCS_update(SCS *self, PyObject *args) {
 /* Deallocate SCS object */
 static scs_int SCS_finish(SCS *self) {
   if (self->work) {
-    /* Acquire lock to ensure no concurrent solve/update is in progress. */
+    /* Acquire lock to ensure no concurrent solve/update is in progress.
+     * We don't check the return value here because this is the dealloc
+     * path — the object must be cleaned up regardless of lock status,
+     * and there is no Python caller to return an error to. */
     if (self->lock) {
       PyThread_acquire_lock(self->lock, WAIT_LOCK);
     }
diff --git a/test/conftest.py b/test/conftest.py
new file mode 100644
index 00000000..91e7ebbb
--- /dev/null
+++ b/test/conftest.py
@@ -0,0 +1,30 @@
+"""Shared pytest configuration for SCS tests."""
+import sys
+
+import pytest
+
+# Detect free-threaded build (3.13t+)
+NOGIL_BUILD = hasattr(sys, "_is_gil_enabled")
+gil_enabled_at_start = True
+if NOGIL_BUILD:
+    gil_enabled_at_start = sys._is_gil_enabled()
+
+
+def pytest_terminal_summary(terminalreporter, exitstatus, config):
+    """Fail the test run if the GIL was re-enabled during tests.
+
+    Modeled after NumPy's conftest.py — if a C extension module is missing
+    the Py_MOD_GIL_NOT_USED declaration, CPython will silently re-enable
+    the GIL at import time, defeating free-threading. This hook catches that.
+    """
+    if NOGIL_BUILD and not gil_enabled_at_start and sys._is_gil_enabled():
+        tr = terminalreporter
+        tr.ensure_newline()
+        tr.section("GIL re-enabled", sep="=", red=True, bold=True)
+        tr.line("The GIL was re-enabled at runtime during the tests.")
+        tr.line("This can happen with no test failures if the RuntimeWarning")
+        tr.line("raised by Python when this happens is filtered by a test.")
+        tr.line("")
+        tr.line("Please ensure all C extension modules declare")
+        tr.line("Py_MOD_GIL_NOT_USED via PyUnstable_Module_SetGIL().")
+        pytest.exit("GIL re-enabled during tests", returncode=1)
diff --git a/test/gen_random_cone_prob.py b/test/gen_random_cone_prob.py
index b1707b23..4e9528a6 100644
--- a/test/gen_random_cone_prob.py
+++ b/test/gen_random_cone_prob.py
@@ -6,15 +6,17 @@
 #############################################
 
 
-def gen_feasible(K, n, density):
+def gen_feasible(K, n, density, rng=None):
+    if rng is None:
+        rng = np.random
     m = get_scs_cone_dims(K)
-    z = np.random.randn(m)
+    z = rng.randn(m)
     y = proj_dual_cone(z, K)  # y = s - z;
     s = y - z  # s = proj_cone(z,K)
 
-    A = sparse.rand(m, n, density, format="csc")
-    A.data = np.random.randn(A.nnz)
-    x = np.random.randn(n)
+    A = sparse.rand(m, n, density, format="csc", random_state=rng)
+    A.data = rng.randn(A.nnz)
+    x = rng.randn(n)
     c = -np.transpose(A).dot(y)
     b = A.dot(x) + s
 
@@ -22,36 +24,40 @@ def gen_feasible(K, n, density):
     return data, np.dot(c, x)
 
 
-def gen_infeasible(K, n):
+def gen_infeasible(K, n, rng=None):
+    if rng is None:
+        rng = np.random
     m = get_scs_cone_dims(K)
 
-    z = np.random.randn(m)
+    z = rng.randn(m)
     y = proj_dual_cone(z, K)  # y = s - z;
-    A = np.random.randn(m, n)
+    A = rng.randn(m, n)
     A = (
         A - np.outer(y, np.transpose(A).dot(y)) / np.linalg.norm(y) ** 2
     )  # dense...
 
-    b = np.random.randn(m)
+    b = rng.randn(m)
     b = -b / np.dot(b, y)
 
-    data = {"A": sparse.csc_matrix(A), "b": b, "c": np.random.randn(n)}
+    data = {"A": sparse.csc_matrix(A), "b": b, "c": rng.randn(n)}
     return data
 
 
-def gen_unbounded(K, n):
+def gen_unbounded(K, n, rng=None):
+    if rng is None:
+        rng = np.random
     m = get_scs_cone_dims(K)
 
-    z = np.random.randn(m)
+    z = rng.randn(m)
     s = proj_cone(z, K)
-    A = np.random.randn(m, n)
-    x = np.random.randn(n)
+    A = rng.randn(m, n)
+    x = rng.randn(n)
     A = A - np.outer(s + A.dot(x), x) / np.linalg.norm(x) ** 2
     # dense...
-    c = np.random.randn(n)
+    c = rng.randn(n)
     c = -c / np.dot(c, x)
 
-    data = {"A": sparse.csc_matrix(A), "b": np.random.randn(m), "c": c}
+    data = {"A": sparse.csc_matrix(A), "b": rng.randn(m), "c": c}
     return data
 
 
diff --git a/test/lsan-suppressions.txt b/test/lsan-suppressions.txt
new file mode 100644
index 00000000..c25eebed
--- /dev/null
+++ b/test/lsan-suppressions.txt
@@ -0,0 +1,35 @@
+# Leak Sanitizer suppressions for scs-python
+#
+# CPython does not free all memory at exit (interned strings, singletons,
+# pymalloc arenas, etc.). These are not real leaks.
+
+# CPython memory allocator internals
+leak:_PyMem
+leak:_PyObject
+leak:PyMem
+leak:PyObject_Malloc
+leak:PyObject_Realloc
+leak:PyObject_Calloc
+
+# CPython interpreter/module init
+leak:_Py_InitializeMain
+leak:Py_InitializeFromConfig
+leak:PyModule
+leak:PyImport
+leak:_PyImport
+leak:PyType_Ready
+leak:type_new
+leak:_PyType
+
+# CPython interned strings and caches
+leak:PyUnicode
+leak:_PyUnicode
+leak:intern_string_constants
+leak:PyCode_New
+
+# NumPy/SciPy internal allocations
+leak:numpy
+leak:scipy
+
+# Meson build probes (run during pip install)
+leak:meson
diff --git a/test/solve_random_cone_prob.py b/test/solve_random_cone_prob.py
index 60cc4c51..8de0a603 100644
--- a/test/solve_random_cone_prob.py
+++ b/test/solve_random_cone_prob.py
@@ -18,13 +18,13 @@ def main():
         pass
 
     for linear_solver in solvers:
-        np.random.seed(1)
-        solve_feasible(linear_solver)
-        solve_infeasible(linear_solver)
-        solve_unbounded(linear_solver)
+        rng = np.random.RandomState(3000)
+        solve_feasible(linear_solver, rng)
+        solve_infeasible(linear_solver, rng)
+        solve_unbounded(linear_solver, rng)
 
 
-def solve_feasible(linear_solver):
+def solve_feasible(linear_solver, rng):
     # cone:
     K = {
         "z": 10,
@@ -36,7 +36,7 @@ def solve_feasible(linear_solver):
         "p": [-0.25, 0.5, 0.75, -0.33],
     }
     m = tools.get_scs_cone_dims(K)
-    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.01)
+    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.01, rng=rng)
     params = {"normalize": True, "scale": 5}
 
     sol = scs.solve(data, K, linear_solver=linear_solver, **params)
@@ -47,7 +47,7 @@ def solve_feasible(linear_solver):
     print("dual error = ", (-np.dot(data["b"], y) - p_star) / p_star)
 
 
-def solve_infeasible(linear_solver):
+def solve_infeasible(linear_solver, rng):
     K = {
         "z": 10,
         "l": 15,
@@ -58,12 +58,12 @@ def solve_infeasible(linear_solver):
         "p": [-0.25, 0.5, 0.75, -0.33],
     }
     m = tools.get_scs_cone_dims(K)
-    data = tools.gen_infeasible(K, n=m // 3)
+    data = tools.gen_infeasible(K, n=m // 3, rng=rng)
     params = {"normalize": True, "scale": 0.5}
     sol = scs.solve(data, K, linear_solver=linear_solver, **params)
 
 
-def solve_unbounded(linear_solver):
+def solve_unbounded(linear_solver, rng):
     K = {
         "z": 10,
         "l": 15,
@@ -74,7 +74,7 @@ def solve_unbounded(linear_solver):
         "p": [-0.25, 0.5, 0.75, -0.33],
     }
     m = tools.get_scs_cone_dims(K)
-    data = tools.gen_unbounded(K, n=m // 3)
+    data = tools.gen_unbounded(K, n=m // 3, rng=rng)
     params = {"normalize": True, "scale": 0.5}
     sol = scs.solve(data, K, linear_solver=linear_solver, **params)
 
diff --git a/test/test_free_threading.py b/test/test_free_threading.py
new file mode 100644
index 00000000..95fd849c
--- /dev/null
+++ b/test/test_free_threading.py
@@ -0,0 +1,985 @@
+"""Tests for free-threading (no-GIL) support.
+
+These tests verify that scs works correctly when called concurrently from
+multiple threads. On GIL-enabled builds, the tests still run but the
+concurrency is serialized by the GIL. On free-threaded builds (3.13t+),
+these tests exercise true parallel execution.
+"""
+import sys
+import threading
+from concurrent.futures import ThreadPoolExecutor, as_completed
+
+import numpy as np
+import pytest
+import scipy.sparse as sp
+from numpy.testing import assert_almost_equal
+
+import scs
+
+
+# ---------------------------------------------------------------------------
+# Helpers
+# ---------------------------------------------------------------------------
+
+def _make_simple_lp():
+    """Simple LP: max x s.t. 0 <= x <= 1. Optimal x=1."""
+    A = sp.csc_matrix([1.0, -1.0]).T.tocsc()
+    b = np.array([1.0, 0.0])
+    c = np.array([-1.0])
+    data = {"A": A, "b": b, "c": c}
+    cone = {"l": 2}
+    return data, cone, 1.0  # expected x[0]
+
+
+def _make_socp():
+    """SOCP: max x s.t. |x| <= 1 (SOC constraint). Optimal x=0.5."""
+    A = sp.csc_matrix([1.0, -1.0]).T.tocsc()
+    b = np.array([1.0, 0.0])
+    c = np.array([-1.0])
+    data = {"A": A, "b": b, "c": c}
+    cone = {"q": [2]}
+    return data, cone, 0.5  # expected x[0]
+
+
+def _make_larger_lp(n=20, seed=42):
+    """A larger LP with n variables for more substantial computation.
+
+    max  c'x
+    s.t. Ax <= b, x >= 0
+
+    Constructed so optimal solution has x[0] > 0.
+    """
+    rng = np.random.RandomState(seed)
+    m = 3 * n
+    # Feasible LP: random A, b = A @ ones(n) + slack
+    A_dense = rng.randn(m, n)
+    x_feas = np.abs(rng.randn(n)) + 0.1
+    b_ineq = A_dense @ x_feas + np.abs(rng.randn(m)) + 0.1
+
+    # Stack: inequality constraints (Ax <= b) and nonnegativity (x >= 0)
+    A_full = sp.vstack([
+        sp.csc_matrix(A_dense),
+        sp.eye(n, format="csc") * -1.0,
+    ], format="csc")
+    b_full = np.concatenate([b_ineq, np.zeros(n)])
+    c = -np.abs(rng.randn(n))  # negative for maximization
+
+    data = {"A": A_full, "b": b_full, "c": c}
+    cone = {"l": m + n}
+    return data, cone
+
+
+def _solve_and_check(solver, expected_x0, decimal=2):
+    """Call solve() and check result. Returns the solution dict."""
+    sol = solver.solve()
+    assert sol["info"]["status_val"] in (1, 2), (
+        f"Unexpected status: {sol['info']['status']}"
+    )
+    if expected_x0 is not None:
+        assert_almost_equal(sol["x"][0], expected_x0, decimal=decimal)
+    return sol
+
+
+# ---------------------------------------------------------------------------
+# Test: concurrent solves on INDEPENDENT instances (the primary use case)
+# ---------------------------------------------------------------------------
+
+NUM_THREADS = 8
+
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestConcurrentIndependentInstances:
+    """Each thread creates and solves its own SCS instance.
+
+    This is the main use case: different optimization problems solved in
+    parallel, each with its own data and workspace.
+    """
+
+    def test_concurrent_lp_solves(self):
+        """Multiple threads each solve the same simple LP independently."""
+        data, cone, expected = _make_simple_lp()
+
+        def worker():
+            solver = scs.SCS(data, cone, verbose=False)
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return sol["x"][0]
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == NUM_THREADS
+        for r in results:
+            assert_almost_equal(r, expected, decimal=2)
+
+    def test_concurrent_socp_solves(self):
+        """Multiple threads each solve the same SOCP independently."""
+        data, cone, expected = _make_socp()
+
+        def worker():
+            solver = scs.SCS(data, cone, verbose=False)
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return sol["x"][0]
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == NUM_THREADS
+
+    def test_concurrent_mixed_problems(self):
+        """Threads solve different problem types concurrently."""
+        lp_data, lp_cone, lp_expected = _make_simple_lp()
+        socp_data, socp_cone, socp_expected = _make_socp()
+
+        def lp_worker(i):
+            solver = scs.SCS(lp_data, lp_cone, verbose=False)
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], lp_expected, decimal=2)
+            return ("lp", i, sol["x"][0])
+
+        def socp_worker(i):
+            solver = scs.SCS(socp_data, socp_cone, verbose=False)
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], socp_expected, decimal=2)
+            return ("socp", i, sol["x"][0])
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = []
+            for i in range(NUM_THREADS):
+                if i % 2 == 0:
+                    futures.append(pool.submit(lp_worker, i))
+                else:
+                    futures.append(pool.submit(socp_worker, i))
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == NUM_THREADS
+
+    def test_concurrent_direct_and_indirect(self):
+        """Threads use different solver backends (direct vs indirect)."""
+        data, cone, expected = _make_simple_lp()
+
+        def worker(linear_solver):
+            solver = scs.SCS(
+                data, cone, linear_solver=linear_solver, verbose=False
+            )
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return sol["x"][0]
+
+        backends = [scs.LinearSolver.QDLDL, scs.LinearSolver.INDIRECT]
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [
+                pool.submit(worker, backends[i % 2]) for i in range(NUM_THREADS)
+            ]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == NUM_THREADS
+
+    def test_concurrent_larger_problems(self):
+        """Threads solve larger problems that take more computation."""
+        data, cone = _make_larger_lp(n=20)
+
+        def worker(seed):
+            solver = scs.SCS(
+                data, cone, verbose=False, max_iters=5000
+            )
+            sol = solver.solve()
+            assert sol["info"]["status_val"] in (1, 2)
+            return sol["info"]["status_val"]
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker, i) for i in range(NUM_THREADS)]
+            results = [f.result(timeout=60) for f in futures]
+
+        assert len(results) == NUM_THREADS
+
+    def test_many_sequential_solves_per_thread(self):
+        """Each thread creates and solves multiple problems sequentially."""
+        data, cone, expected = _make_simple_lp()
+        solves_per_thread = 10
+
+        def worker():
+            results = []
+            for _ in range(solves_per_thread):
+                solver = scs.SCS(data, cone, verbose=False)
+                sol = solver.solve()
+                assert sol["info"]["status_val"] == 1
+                results.append(sol["x"][0])
+            return results
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            all_results = [f.result(timeout=60) for f in futures]
+
+        assert len(all_results) == NUM_THREADS
+        for thread_results in all_results:
+            assert len(thread_results) == solves_per_thread
+            for r in thread_results:
+                assert_almost_equal(r, expected, decimal=2)
+
+
+# ---------------------------------------------------------------------------
+# Test: concurrent solves on a SHARED instance
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestConcurrentSharedInstance:
+    """Multiple threads call solve() on the same SCS object.
+
+    The per-instance PyMutex (in free-threaded builds) serializes access
+    to the underlying C workspace, so results should still be correct.
+    """
+
+    def test_shared_instance_concurrent_solve(self):
+        """Multiple threads solve on the same SCS instance."""
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+
+        def worker():
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return sol["x"][0]
+
+        with ThreadPoolExecutor(max_workers=4) as pool:
+            futures = [pool.submit(worker) for _ in range(4)]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == 4
+        for r in results:
+            assert_almost_equal(r, expected, decimal=2)
+
+    def test_shared_instance_repeated_concurrent_solve(self):
+        """Repeated rounds of concurrent solves on a shared instance."""
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+
+        for round_num in range(5):
+            def worker():
+                sol = solver.solve()
+                assert sol["info"]["status_val"] == 1
+                return sol["x"][0]
+
+            with ThreadPoolExecutor(max_workers=4) as pool:
+                futures = [pool.submit(worker) for _ in range(4)]
+                results = [f.result(timeout=30) for f in futures]
+
+            for r in results:
+                assert_almost_equal(r, expected, decimal=2)
+
+
+# ---------------------------------------------------------------------------
+# Test: concurrent solve + update sequences
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestConcurrentSolveUpdate:
+    """Concurrent solve and update operations."""
+
+    def test_shared_instance_concurrent_solve_and_update(self):
+        """Threads simultaneously call solve() and update() on the same instance.
+
+        The per-instance lock serializes access, so no data race should occur.
+        We can't predict which operation runs first, but none should crash or
+        corrupt state.
+        """
+        data, cone, _ = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+        barrier = threading.Barrier(NUM_THREADS)
+        errors = []
+        lock = threading.Lock()
+
+        def solve_worker():
+            try:
+                barrier.wait(timeout=10)
+                sol = solver.solve()
+                assert sol["info"]["status_val"] in (1, 2)
+            except Exception as e:
+                with lock:
+                    errors.append(f"solve: {e}")
+
+        def update_worker():
+            try:
+                barrier.wait(timeout=10)
+                # Update with slightly perturbed b
+                solver.update(b=np.array([1.0, 0.0]))
+            except Exception as e:
+                with lock:
+                    errors.append(f"update: {e}")
+
+        threads = []
+        for i in range(NUM_THREADS):
+            if i % 2 == 0:
+                t = threading.Thread(target=solve_worker)
+            else:
+                t = threading.Thread(target=update_worker)
+            threads.append(t)
+            t.start()
+
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive(), "Thread did not complete"
+
+        assert len(errors) == 0, f"Errors in threads: {errors}"
+
+        # Solver should still be usable after the concurrent barrage
+        sol = solver.solve()
+        assert sol["info"]["status_val"] in (1, 2)
+
+    def test_concurrent_update_sequences(self):
+        """Multiple threads each do solve -> update -> solve independently."""
+        data, cone, _ = _make_simple_lp()
+
+        def worker():
+            solver = scs.SCS(data, cone, verbose=False)
+
+            # Initial solve: max x s.t. 0 <= x <= 1 => x=1
+            sol1 = solver.solve()
+            assert sol1["info"]["status_val"] == 1
+            assert_almost_equal(sol1["x"][0], 1.0, decimal=2)
+
+            # Update c to minimize: min x s.t. 0 <= x <= 1 => x=0
+            solver.update(c=np.array([1.0]))
+            sol2 = solver.solve()
+            assert sol2["info"]["status_val"] == 1
+            assert_almost_equal(sol2["x"][0], 0.0, decimal=2)
+
+            # Update b: max x s.t. -1 <= x <= 1 => x=-1 (still minimizing)
+            solver.update(b=np.array([1.0, 1.0]))
+            sol3 = solver.solve()
+            assert sol3["info"]["status_val"] == 1
+            assert_almost_equal(sol3["x"][0], -1.0, decimal=2)
+
+            return True
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=60) for f in futures]
+
+        assert all(results)
+        assert len(results) == NUM_THREADS
+
+    def test_concurrent_warm_start(self):
+        """Multiple threads each do warm-started solves independently."""
+        data, cone, expected = _make_simple_lp()
+
+        def worker():
+            solver = scs.SCS(data, cone, verbose=False)
+
+            # Cold solve
+            sol1 = solver.solve()
+            assert_almost_equal(sol1["x"][0], expected, decimal=2)
+
+            # Warm-started solve (uses previous solution)
+            sol2 = solver.solve()
+            assert_almost_equal(sol2["x"][0], expected, decimal=2)
+            # Warm start should converge faster
+            assert sol2["info"]["iter"] <= sol1["info"]["iter"]
+
+            # Explicit warm start vectors
+            sol3 = solver.solve(
+                x=np.array([0.9]),
+                y=sol2["y"],
+                s=sol2["s"],
+            )
+            assert_almost_equal(sol3["x"][0], expected, decimal=2)
+
+            return True
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=60) for f in futures]
+
+        assert all(results)
+
+
+# ---------------------------------------------------------------------------
+# Test: legacy solve() function under concurrency
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestConcurrentLegacySolve:
+    """Test the legacy scs.solve() function from multiple threads."""
+
+    def test_concurrent_legacy_solve(self):
+        """Multiple threads use the legacy solve() interface."""
+        data, cone, expected = _make_simple_lp()
+
+        def worker():
+            sol = scs.solve(data, cone, verbose=False)
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return sol["x"][0]
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert len(results) == NUM_THREADS
+
+
+# ---------------------------------------------------------------------------
+# Test: thread creation and destruction stress test
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestThreadStress:
+    """Stress tests with many short-lived threads."""
+
+    def test_rapid_thread_creation(self):
+        """Rapidly create and destroy threads, each doing a quick solve."""
+        data, cone, expected = _make_simple_lp()
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                solver = scs.SCS(data, cone, verbose=False)
+                sol = solver.solve()
+                assert_almost_equal(sol["x"][0], expected, decimal=2)
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = []
+        for _ in range(50):
+            t = threading.Thread(target=worker)
+            threads.append(t)
+            t.start()
+
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive(), "Thread did not complete"
+
+        assert len(errors) == 0, f"Errors in threads: {errors}"
+
+    def test_concurrent_construction_and_solve(self):
+        """Threads concurrently construct SCS objects and solve."""
+        data, cone, expected = _make_simple_lp()
+
+        barrier = threading.Barrier(NUM_THREADS)
+
+        def worker():
+            # All threads start constructing at roughly the same time
+            barrier.wait(timeout=10)
+            solver = scs.SCS(data, cone, verbose=False)
+            sol = solver.solve()
+            assert sol["info"]["status_val"] == 1
+            assert_almost_equal(sol["x"][0], expected, decimal=2)
+            return True
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker) for _ in range(NUM_THREADS)]
+            results = [f.result(timeout=30) for f in futures]
+
+        assert all(results)
+
+
+# ---------------------------------------------------------------------------
+# Test: result isolation (no cross-thread data corruption)
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestResultIsolation:
+    """Verify that results from concurrent solves don't get mixed up."""
+
+    def test_result_vectors_are_independent(self):
+        """Each thread's result arrays are independent memory."""
+        data, cone, expected = _make_simple_lp()
+        results_lock = threading.Lock()
+        all_results = []
+
+        def worker(thread_id):
+            solver = scs.SCS(data, cone, verbose=False)
+            sol = solver.solve()
+            # Store the actual array objects
+            with results_lock:
+                all_results.append({
+                    "id": thread_id,
+                    "x": sol["x"].copy(),
+                    "y": sol["y"].copy(),
+                    "s": sol["s"].copy(),
+                })
+            return True
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = [pool.submit(worker, i) for i in range(NUM_THREADS)]
+            [f.result(timeout=30) for f in futures]
+
+        # All results should be correct and identical in value
+        assert len(all_results) == NUM_THREADS
+        for r in all_results:
+            assert_almost_equal(r["x"][0], expected, decimal=2)
+
+    def test_different_problems_correct_results(self):
+        """Threads solving different problems get their own correct answers."""
+        # Problem 1: max x, 0 <= x <= 1 => x=1
+        data1, cone1, exp1 = _make_simple_lp()
+        # Problem 2: max x, |x| <= 1 (SOC) => x=0.5
+        data2, cone2, exp2 = _make_socp()
+
+        def worker_lp():
+            solver = scs.SCS(data1, cone1, verbose=False)
+            sol = solver.solve()
+            assert_almost_equal(sol["x"][0], exp1, decimal=2)
+            return ("lp", sol["x"][0])
+
+        def worker_socp():
+            solver = scs.SCS(data2, cone2, verbose=False)
+            sol = solver.solve()
+            assert_almost_equal(sol["x"][0], exp2, decimal=2)
+            return ("socp", sol["x"][0])
+
+        with ThreadPoolExecutor(max_workers=NUM_THREADS) as pool:
+            futures = []
+            for i in range(NUM_THREADS):
+                if i % 2 == 0:
+                    futures.append(pool.submit(worker_lp))
+                else:
+                    futures.append(pool.submit(worker_socp))
+
+            for f in as_completed(futures, timeout=30):
+                kind, val = f.result()
+                if kind == "lp":
+                    assert_almost_equal(val, exp1, decimal=2)
+                else:
+                    assert_almost_equal(val, exp2, decimal=2)
+
+
+# ---------------------------------------------------------------------------
+# Stress tests targeting specific thread-safety bugs
+# ---------------------------------------------------------------------------
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestBorrowedRefSafety:
+    """Tests exercising PyDict/PyList ref safety during concurrent init.
+
+    Before the fix, PyDict_GetItemString / PyList_GetItem returned borrowed
+    references that could be invalidated if another thread mutated the
+    container concurrently. These tests share mutable containers across
+    threads to exercise the strong-ref (PyDict_GetItemStringRef /
+    PyList_GetItemRef) code paths.
+    """
+
+    def test_shared_cone_dict_concurrent_init(self):
+        """Multiple threads construct SCS instances from the same cone dict.
+
+        The cone dict is read by the C extension during __init__. With
+        borrowed refs, a concurrent mutation to the dict could invalidate
+        a pointer mid-parse. Strong refs prevent this.
+        """
+        A = sp.csc_matrix([1.0, -1.0]).T.tocsc()
+        b = np.array([1.0, 0.0])
+        c = np.array([-1.0])
+        data = {"A": A, "b": b, "c": c}
+        cone = {"l": 2}  # shared across all threads
+
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                solver = scs.SCS(data, cone, verbose=False)
+                sol = solver.solve()
+                assert sol["info"]["status_val"] == 1
+                assert_almost_equal(sol["x"][0], 1.0, decimal=2)
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = [threading.Thread(target=worker) for _ in range(20)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+    def test_shared_cone_with_list_values_concurrent_init(self):
+        """Cone dict with list values (SOC dimensions) shared across threads.
+
+        Exercises PyList_GetItemRef in get_cone_arr_dim — the list elements
+        are read one-by-one during parsing.
+        """
+        # Feasible SOCP: minimize c'x s.t. ||x|| <= t (one SOC of size 3)
+        # A = -I (3x2 + slack), b = 0, cone q=[3]
+        # With an extra nonneg cone to bound things
+        A = sp.vstack([
+            -sp.eye(3, n=2, format="csc"),  # SOC: (s0, s1, s2) = b - Ax
+            sp.eye(2, format="csc"),         # nonneg: x >= 0
+        ], format="csc")
+        b = np.array([1.0, 0.0, 0.0, 0.0, 0.0])
+        c = np.array([-1.0, 0.0])
+        data = {"A": A, "b": b, "c": c}
+        cone = {"q": [3], "l": 2}  # list value — exercises PyList_GetItemRef
+
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                solver = scs.SCS(data, cone, verbose=False)
+                sol = solver.solve()
+                assert sol["info"]["status_val"] in (1, 2)
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = [threading.Thread(target=worker) for _ in range(20)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+    def test_shared_cone_with_float_list_values_concurrent_init(self):
+        """Cone dict with float list values (power cone exponents) shared
+        across threads.
+
+        Exercises PyList_GetItemRef in get_cone_float_arr — the float list
+        elements are read one-by-one during parsing.
+        """
+        # Power cone: (x1, x2, x3) s.t. |x3| <= x1^p * x2^(1-p), x1,x2 >= 0
+        # Each power cone triple is size 3, with exponent p
+        m = 3
+        n = 3
+        A = -sp.eye(m, n=n, format="csc")
+        b = np.array([1.0, 1.0, 0.0])
+        c = np.array([0.0, 0.0, -1.0])
+        data = {"A": A, "b": b, "c": c}
+        cone = {"p": [0.5]}  # float list — exercises PyList_GetItemRef in get_cone_float_arr
+
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                solver = scs.SCS(data, cone, verbose=False, max_iters=10000)
+                sol = solver.solve()
+                # Any valid status is fine — we're testing thread safety not convergence
+                assert sol["info"]["status_val"] in (1, 2, -1, -2, -7)
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = [threading.Thread(target=worker) for _ in range(20)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestTOCTOURaceSafety:
+    """Tests exercising the TOCTOU fix on self->work.
+
+    Before the fix, SCS_solve and SCS_update checked self->work without
+    holding the lock. A concurrent dealloc (SCS_finish) could set
+    self->work = NULL between the check and the lock acquisition.
+    """
+
+    def test_rapid_create_solve_destroy(self):
+        """Rapidly create, solve, and let SCS instances be garbage collected.
+
+        Many short-lived instances stress the init/finish paths. On
+        free-threaded builds, the dealloc can race with in-progress solves
+        if the TOCTOU check isn't under the lock.
+        """
+        import gc
+        data, cone, expected = _make_simple_lp()
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                for _ in range(50):
+                    solver = scs.SCS(data, cone, verbose=False)
+                    sol = solver.solve()
+                    assert sol["info"]["status_val"] == 1
+                    assert_almost_equal(sol["x"][0], expected, decimal=2)
+                    del solver
+                gc.collect()
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = [threading.Thread(target=worker) for _ in range(NUM_THREADS)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=60)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+    def test_solve_after_del_raises_or_succeeds(self):
+        """Verify solve doesn't crash if workspace is gone.
+
+        This is a single-threaded sanity check: calling solve on a
+        properly-constructed instance should work, and the object shouldn't
+        crash during cleanup.
+        """
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+        sol = solver.solve()
+        assert sol["info"]["status_val"] == 1
+        assert_almost_equal(sol["x"][0], expected, decimal=2)
+        # Second solve should also work (warm start)
+        sol2 = solver.solve()
+        assert sol2["info"]["status_val"] == 1
+
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestConcurrentSolveUpdateStress:
+    """Heavy stress tests for concurrent solve + update on shared instances."""
+
+    def test_solve_update_barrage_shared_instance(self):
+        """Many threads hammer solve() and update() on the same instance.
+
+        This is a more aggressive version of
+        test_shared_instance_concurrent_solve_and_update, with more threads
+        and repeated rounds.
+        """
+        data, cone, _ = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+        errors = []
+        lock = threading.Lock()
+
+        def solve_worker():
+            try:
+                for _ in range(10):
+                    sol = solver.solve()
+                    assert sol["info"]["status_val"] in (1, 2)
+            except Exception as e:
+                with lock:
+                    errors.append(f"solve: {e}")
+
+        def update_worker():
+            try:
+                for _ in range(10):
+                    solver.update(b=np.array([1.0, 0.0]))
+            except Exception as e:
+                with lock:
+                    errors.append(f"update: {e}")
+
+        threads = []
+        for i in range(16):
+            if i % 2 == 0:
+                threads.append(threading.Thread(target=solve_worker))
+            else:
+                threads.append(threading.Thread(target=update_worker))
+
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=60)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+        # Instance should still be usable
+        sol = solver.solve()
+        assert sol["info"]["status_val"] in (1, 2)
+
+    def test_concurrent_update_different_data(self):
+        """Threads update with different b/c vectors concurrently.
+
+        After the barrage, we verify the solver still produces valid results
+        (not corrupted state).
+        """
+        data, cone, _ = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+        errors = []
+        lock = threading.Lock()
+        rng = np.random.RandomState(42)
+
+        # Pre-generate update vectors (each thread gets its own)
+        b_updates = [np.array([rng.uniform(0.5, 2.0), 0.0]) for _ in range(NUM_THREADS)]
+        c_updates = [np.array([rng.uniform(-2.0, -0.5)]) for _ in range(NUM_THREADS)]
+
+        barrier = threading.Barrier(NUM_THREADS)
+
+        def worker(tid):
+            try:
+                barrier.wait(timeout=10)
+                for _ in range(5):
+                    solver.update(b=b_updates[tid], c=c_updates[tid])
+                    sol = solver.solve()
+                    # Just check it didn't crash or return garbage status
+                    assert sol["info"]["status_val"] in (1, 2, -2, -7)
+            except Exception as e:
+                with lock:
+                    errors.append(f"thread {tid}: {e}")
+
+        threads = [threading.Thread(target=worker, args=(i,)) for i in range(NUM_THREADS)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=60)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+    def test_warm_start_under_contention(self):
+        """Concurrent warm-started solves on a shared instance.
+
+        Warm starts read and write the sol struct, which is protected by the
+        per-instance lock. This test stresses that path.
+        """
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+
+        # Initial cold solve to populate warm-start data
+        sol0 = solver.solve()
+        assert sol0["info"]["status_val"] == 1
+
+        errors = []
+        lock = threading.Lock()
+
+        def worker():
+            try:
+                for _ in range(10):
+                    sol = solver.solve(
+                        x=np.array([0.9]),
+                        y=sol0["y"],
+                        s=sol0["s"],
+                    )
+                    assert sol["info"]["status_val"] == 1
+                    assert_almost_equal(sol["x"][0], expected, decimal=2)
+            except Exception as e:
+                with lock:
+                    errors.append(str(e))
+
+        threads = [threading.Thread(target=worker) for _ in range(NUM_THREADS)]
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=60)
+            assert not t.is_alive()
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
+class TestErrorPathContention:
+    """Tests that lock is properly released on error paths.
+
+    If a thread hits an error (e.g. bad warm-start dimensions) inside
+    the locked section, it must release the lock before returning. If it
+    doesn't, subsequent calls from other threads will deadlock.
+    """
+
+    def test_bad_warm_start_does_not_deadlock(self):
+        """One thread passes wrong-dimension warm-start vectors while others
+        solve normally. The error path must release the lock.
+        """
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+
+        # Do an initial solve so warm-start data is populated
+        solver.solve()
+
+        errors = []
+        lock = threading.Lock()
+        barrier = threading.Barrier(NUM_THREADS)
+
+        def good_worker():
+            try:
+                barrier.wait(timeout=10)
+                for _ in range(10):
+                    sol = solver.solve()
+                    assert sol["info"]["status_val"] == 1
+                    assert_almost_equal(sol["x"][0], expected, decimal=2)
+            except Exception as e:
+                with lock:
+                    errors.append(f"good: {e}")
+
+        def bad_worker():
+            try:
+                barrier.wait(timeout=10)
+                for _ in range(10):
+                    try:
+                        # Wrong dimension x -- should raise ValueError
+                        solver.solve(x=np.array([1.0, 2.0, 3.0]))
+                    except ValueError:
+                        pass  # expected -- lock must have been released
+            except Exception as e:
+                with lock:
+                    errors.append(f"bad: {e}")
+
+        threads = []
+        for i in range(NUM_THREADS):
+            if i == 0:
+                threads.append(threading.Thread(target=bad_worker))
+            else:
+                threads.append(threading.Thread(target=good_worker))
+
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive(), "Thread deadlocked -- lock not released on error?"
+
+        assert len(errors) == 0, f"Errors: {errors}"
+
+    def test_bad_update_does_not_deadlock(self):
+        """One thread passes wrong-dimension update vectors while others
+        solve normally. The error path in SCS_update must release the lock.
+        """
+        data, cone, expected = _make_simple_lp()
+        solver = scs.SCS(data, cone, verbose=False)
+        errors = []
+        lock = threading.Lock()
+        barrier = threading.Barrier(NUM_THREADS)
+
+        def good_worker():
+            try:
+                barrier.wait(timeout=10)
+                for _ in range(10):
+                    sol = solver.solve()
+                    assert sol["info"]["status_val"] in (1, 2)
+            except Exception as e:
+                with lock:
+                    errors.append(f"good: {e}")
+
+        def bad_worker():
+            try:
+                barrier.wait(timeout=10)
+                for _ in range(10):
+                    try:
+                        # Wrong dimension b -- should raise ValueError
+                        solver.update(b=np.array([1.0, 2.0, 3.0]))
+                    except ValueError:
+                        pass  # expected
+            except Exception as e:
+                with lock:
+                    errors.append(f"bad: {e}")
+
+        threads = []
+        for i in range(NUM_THREADS):
+            if i == 0:
+                threads.append(threading.Thread(target=bad_worker))
+            else:
+                threads.append(threading.Thread(target=good_worker))
+
+        for t in threads:
+            t.start()
+        for t in threads:
+            t.join(timeout=30)
+            assert not t.is_alive(), "Thread deadlocked -- lock not released on error?"
+
+        assert len(errors) == 0, f"Errors: {errors}"
diff --git a/test/test_mix_sd_csd_cone.py b/test/test_mix_sd_csd_cone.py
index 23556ca0..f918e791 100644
--- a/test/test_mix_sd_csd_cone.py
+++ b/test/test_mix_sd_csd_cone.py
@@ -2,12 +2,12 @@
 import scs
 import scipy
 import pytest
-    
-def gen_feasible(m, n, p_scale = 0.1):
-    P = p_scale * scipy.sparse.eye(n, format="csc")   
-    A = scipy.sparse.random(m, n, density=0.05, format="csc")
-    c = np.random.randn(n)
-    b = np.random.randn(m)
+
+def gen_feasible(m, n, rng, p_scale = 0.1):
+    P = p_scale * scipy.sparse.eye(n, format="csc")
+    A = scipy.sparse.random(m, n, density=0.05, format="csc", random_state=rng)
+    c = rng.randn(n)
+    b = rng.randn(m)
 
     return (P, A, b, c)
 
@@ -30,12 +30,11 @@ def gen_feasible(m, n, p_scale = 0.1):
 
 @pytest.mark.parametrize("solver_opts", _solver_configs)
 def test_mix_sd_csd_cones(solver_opts):
-    seed = 1234
-    np.random.seed(seed)
+    rng = np.random.RandomState(1234)
     cone = dict(z=1, l=2, s=[3, 4], cs=[5, 4])
     m = int(cone['z'] + cone['l'] + sum([j * (j+1) / 2 for j in cone['s']]) + sum([j * j for j in cone['cs']]))
     n = m
-    (P, A, b, c) = gen_feasible(m, n)
+    (P, A, b, c) = gen_feasible(m, n, rng)
     probdata = dict(P=P, A=A, b=b, c=c)
     sol = scs.solve(probdata, cone, **solver_opts)
     np.testing.assert_equal(sol['info']['status'], 'solved')
diff --git a/test/test_scs_coverage.py b/test/test_scs_coverage.py
index 61519ff3..1b724a52 100644
--- a/test/test_scs_coverage.py
+++ b/test/test_scs_coverage.py
@@ -347,12 +347,12 @@ def test_time_limit_secs_terminates_early():
     sufficiently large or ill-conditioned problem."""
     # Build a slightly larger problem (random feasible LP) where a near-zero
     # time limit causes early termination.
-    np.random.seed(42)
+    rng = np.random.RandomState(42)
     m, n = 50, 30
-    A_large = sp.random(m, n, density=0.3, format="csc")
-    A_large.data = np.random.randn(A_large.nnz)
-    b_large = np.abs(np.random.randn(m)) + 1.0
-    c_large = np.random.randn(n)
+    A_large = sp.random(m, n, density=0.3, format="csc", random_state=rng)
+    A_large.data = rng.randn(A_large.nnz)
+    b_large = np.abs(rng.randn(m)) + 1.0
+    c_large = rng.randn(n)
     data_large = {"A": A_large, "b": b_large, "c": c_large}
     cone_large = {"l": m}
 
@@ -1463,11 +1463,11 @@ def test_solution_y_s_shapes():
 def test_solution_shapes_match_problem_dimensions():
     """x, y, s shapes must match n and m from the problem data."""
     n, m = 3, 5
-    np.random.seed(99)
-    A_r = sp.random(m, n, density=0.8, format="csc")
-    A_r.data = np.random.randn(A_r.nnz) * 0.1
-    b_r = np.abs(np.random.randn(m)) + 1.0
-    c_r = np.random.randn(n)
+    rng = np.random.RandomState(99)
+    A_r = sp.random(m, n, density=0.8, format="csc", random_state=rng)
+    A_r.data = rng.randn(A_r.nnz) * 0.1
+    b_r = np.abs(rng.randn(m)) + 1.0
+    c_r = rng.randn(n)
     sol = scs.SCS(
         {"A": A_r, "b": b_r, "c": c_r}, {"l": m}, verbose=False
     ).solve()
@@ -1592,10 +1592,10 @@ def test_dual_exp_cone_random_feasible():
     sys.path.insert(0, "test")
     import gen_random_cone_prob as tools
 
-    np.random.seed(7)
+    rng = np.random.RandomState(7)
     K_ed = {"z": 3, "l": 5, "q": [4], "s": [], "ep": 0, "ed": 1, "p": []}
     m_ed = tools.get_scs_cone_dims(K_ed)
-    data, pstar = tools.gen_feasible(K_ed, n=m_ed // 3, density=0.2)
+    data, pstar = tools.gen_feasible(K_ed, n=m_ed // 3, density=0.2, rng=rng)
 
     solver = scs.SCS(
         data, K_ed,
@@ -1627,12 +1627,12 @@ def test_mixed_soc_sdp_random_feasible():
     sys.path.insert(0, "test")
     import gen_random_cone_prob as tools
 
-    np.random.seed(42)
+    rng = np.random.RandomState(42)
     K_mix = {
         "z": 5, "l": 5, "q": [4], "s": [3], "ep": 0, "ed": 0, "p": []
     }
     m_mix = tools.get_scs_cone_dims(K_mix)
-    data, pstar = tools.gen_feasible(K_mix, n=m_mix // 3, density=0.2)
+    data, pstar = tools.gen_feasible(K_mix, n=m_mix // 3, density=0.2, rng=rng)
 
     solver = scs.SCS(
         data, K_mix,
@@ -1767,10 +1767,10 @@ def test_large_random_lp():
     sys.path.insert(0, "test")
     import gen_random_cone_prob as tools
 
-    np.random.seed(3)
+    rng = np.random.RandomState(3)
     K_lp = {"z": 0, "l": 40, "q": [], "s": [], "ep": 0, "ed": 0, "p": []}
     m_lp = tools.get_scs_cone_dims(K_lp)
-    data, pstar = tools.gen_feasible(K_lp, n=20, density=0.3)
+    data, pstar = tools.gen_feasible(K_lp, n=20, density=0.3, rng=rng)
 
     solver = scs.SCS(
         data, K_lp,
@@ -1850,7 +1850,7 @@ def test_cs_cone_mixed_with_real_cones():
     Uses a randomly generated feasible problem (P=ε·I ensures bounded objective).
     Checks the solver returns a valid status.
     """
-    np.random.seed(1234)
+    rng = np.random.RandomState(1234)
     cone_cs = {"z": 1, "l": 2, "s": [3, 4], "cs": [5, 4]}
     # dimension: z=1, l=2, s=[3,4]->(3+4)→6+10=16, cs=[5,4]->25+16=41 => total=60
     m_cs = int(
@@ -1861,10 +1861,10 @@ def test_cs_cone_mixed_with_real_cones():
     )
     n_cs = m_cs
     P_cs = 0.1 * sp.eye(n_cs, format="csc")
-    A_cs = sp.random(m_cs, n_cs, density=0.05, format="csc")
-    A_cs.data = np.random.randn(A_cs.nnz)
-    b_cs = np.random.randn(m_cs)
-    c_cs = np.random.randn(n_cs)
+    A_cs = sp.random(m_cs, n_cs, density=0.05, format="csc", random_state=rng)
+    A_cs.data = rng.randn(A_cs.nnz)
+    b_cs = rng.randn(m_cs)
+    c_cs = rng.randn(n_cs)
 
     solver = scs.SCS(
         {"P": P_cs, "A": A_cs, "b": b_cs, "c": c_cs},
@@ -2078,14 +2078,14 @@ def test_mixed_ep_and_power_cone():
     sys.path.insert(0, "test")
     import gen_random_cone_prob as tools
 
-    np.random.seed(77)
+    rng = np.random.RandomState(77)
     K_mix = {
         "z": 3, "l": 5, "q": [], "s": [],
         "ep": 1, "ed": 0,
         "p": [0.4, -0.6],
     }
     m_mix = tools.get_scs_cone_dims(K_mix)
-    data, pstar = tools.gen_feasible(K_mix, n=m_mix // 3, density=0.2)
+    data, pstar = tools.gen_feasible(K_mix, n=m_mix // 3, density=0.2, rng=rng)
 
     solver = scs.SCS(
         data, K_mix,
@@ -2115,10 +2115,10 @@ def test_two_instances_same_problem_identical_result():
     sys.path.insert(0, "test")
     import gen_random_cone_prob as tools
 
-    np.random.seed(11)
+    rng = np.random.RandomState(11)
     K = {"z": 3, "l": 5, "q": [4], "s": [], "ep": 0, "ed": 0, "p": []}
     m = tools.get_scs_cone_dims(K)
-    data, _ = tools.gen_feasible(K, n=m // 2, density=0.2)
+    data, _ = tools.gen_feasible(K, n=m // 2, density=0.2, rng=rng)
 
     sol1 = scs.SCS(data, K, verbose=False).solve()
     sol2 = scs.SCS(data, K, verbose=False).solve()
@@ -2772,17 +2772,17 @@ def test_P_only_lower_triangular():
 
 def test_all_cone_types_simultaneously():
     """Problem using z, l, q, s, ep, p cone types all at once."""
-    np.random.seed(42)
+    rng = np.random.RandomState(42)
     # Cone dims: z=1, l=2, q=[3], s=[2] (vec dim=3), ep=1 (dim 3), p=[0.5] (dim 3)
     # Total rows: 1 + 2 + 3 + 3 + 3 + 3 = 15
     cone = {"z": 1, "l": 2, "q": [3], "s": [2], "ep": 1, "p": [0.5]}
     m = 15
     n = m
     P = 0.1 * sp.eye(n, format="csc")
-    A = sp.random(m, n, density=0.1, format="csc")
-    A.data = np.random.randn(A.nnz)
-    b = np.random.randn(m)
-    c = np.random.randn(n)
+    A = sp.random(m, n, density=0.1, format="csc", random_state=rng)
+    A.data = rng.randn(A.nnz)
+    b = rng.randn(m)
+    c = rng.randn(n)
     data = {"P": P, "A": A, "b": b, "c": c}
     sol = scs.solve(data, cone, verbose=False, max_iters=50000)
     assert sol["info"]["status"] in ("solved", "solved_inaccurate")
@@ -2897,6 +2897,7 @@ def test_linear_solver_auto_string():
 # ===========================================================================
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_falls_back_to_qdldl():
     """AUTO should fall back to QDLDL when platform-preferred module is missing."""
     from unittest.mock import patch
@@ -2910,6 +2911,7 @@ def fail_import(name):
     assert module is _scs_direct
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_darwin_tries_accelerate():
     """On macOS, AUTO should try _scs_accelerate first."""
     from unittest.mock import patch, MagicMock
@@ -2929,6 +2931,7 @@ def mock_load(name):
     assert module is fake_accel
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_linux_tries_mkl():
     """On Linux, AUTO should try _scs_mkl first."""
     from unittest.mock import patch, MagicMock
@@ -2948,6 +2951,7 @@ def mock_load(name):
     assert module is fake_mkl
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_windows_tries_mkl():
     """On Windows, AUTO should try _scs_mkl first."""
     from unittest.mock import patch, MagicMock
@@ -2967,6 +2971,7 @@ def mock_load(name):
     assert module is fake_mkl
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_darwin_fallback_when_no_accelerate():
     """On macOS, AUTO should fall back to QDLDL when Accelerate is missing."""
     from unittest.mock import patch
@@ -2982,6 +2987,7 @@ def fail_import(name):
     assert module is _scs_direct
 
 
+@pytest.mark.thread_unsafe(reason="patches module-level scs._load_module / scs.sys")
 def test_resolve_auto_linux_fallback_when_no_mkl():
     """On Linux, AUTO should fall back to QDLDL when MKL is missing."""
     from unittest.mock import patch
diff --git a/test/test_scs_rand.py b/test/test_scs_rand.py
index e5433af1..30909437 100644
--- a/test/test_scs_rand.py
+++ b/test/test_scs_rand.py
@@ -52,7 +52,6 @@ def check_unbounded(sol):
     assert_(sol["info"]["status"], "unbounded")
 
 
-np.random.seed(0)
 num_feas = 50
 num_unb = 10
 num_infeas = 10
@@ -92,8 +91,9 @@ def check_unbounded(sol):
 
 @pytest.mark.parametrize("solver_opts", _solver_configs)
 def test_feasible(solver_opts):
+    rng = np.random.RandomState(1000)
     for i in range(num_feas):
-        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         assert_almost_equal(np.dot(data["c"], sol["x"]), p_star, decimal=2)
@@ -102,8 +102,9 @@ def test_feasible(solver_opts):
 
 @pytest.mark.parametrize("solver_opts", _solver_configs)
 def test_infeasible(solver_opts):
+    rng = np.random.RandomState(1001)
     for i in range(num_infeas):
-        data = tools.gen_infeasible(K, n=m // 2)
+        data = tools.gen_infeasible(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         check_infeasible(sol)
@@ -117,8 +118,9 @@ def test_infeasible(solver_opts):
 
 @pytest.mark.parametrize("solver_opts", _unbounded_configs)
 def test_unbounded(solver_opts):
+    rng = np.random.RandomState(1002)
     for i in range(num_unb):
-        data = tools.gen_unbounded(K, n=m // 2)
+        data = tools.gen_unbounded(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         check_unbounded(sol)
diff --git a/test/test_scs_sdp.py b/test/test_scs_sdp.py
index 7c27301d..b0967b46 100644
--- a/test/test_scs_sdp.py
+++ b/test/test_scs_sdp.py
@@ -52,7 +52,6 @@ def check_unbounded(sol):
     assert_(sol["info"]["status"], "unbounded")
 
 
-np.random.seed(0)
 num_feas = 50
 num_unb = 10
 num_infeas = 10
@@ -92,8 +91,9 @@ def check_unbounded(sol):
 
 @pytest.mark.parametrize("solver_opts", _solver_configs)
 def test_feasible(solver_opts):
+    rng = np.random.RandomState(2000)
     for i in range(num_feas):
-        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         assert_almost_equal(np.dot(data["c"], sol["x"]), p_star, decimal=2)
@@ -102,8 +102,9 @@ def test_feasible(solver_opts):
 
 @pytest.mark.parametrize("solver_opts", _solver_configs)
 def test_infeasible(solver_opts):
+    rng = np.random.RandomState(2001)
     for i in range(num_infeas):
-        data = tools.gen_infeasible(K, n=m // 2)
+        data = tools.gen_infeasible(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         check_infeasible(sol)
@@ -117,8 +118,9 @@ def test_infeasible(solver_opts):
 
 @pytest.mark.parametrize("solver_opts", _unbounded_configs)
 def test_unbounded(solver_opts):
+    rng = np.random.RandomState(2002)
     for i in range(num_unb):
-        data = tools.gen_unbounded(K, n=m // 2)
+        data = tools.gen_unbounded(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, **solver_opts, **opts)
         sol = solver.solve()
         check_unbounded(sol)
diff --git a/test/test_solve_random_cone_prob.py b/test/test_solve_random_cone_prob.py
index 2a27b37e..a3882e12 100644
--- a/test/test_solve_random_cone_prob.py
+++ b/test/test_solve_random_cone_prob.py
@@ -29,8 +29,6 @@ def import_error(msg):
 except ImportError:
     pass
 
-np.random.seed(1)
-
 # cone:
 K = {
     "z": 10,
@@ -47,7 +45,8 @@ def import_error(msg):
 
 @pytest.mark.parametrize("linear_solver", solvers)
 def test_solve_feasible(linear_solver):
-    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+    rng = np.random.RandomState(3000)
+    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
     solver = scs.SCS(data, K, linear_solver=linear_solver, **params)
     sol = solver.solve()
     x = sol["x"]
@@ -68,7 +67,8 @@ def test_solve_feasible(linear_solver):
 
 @pytest.mark.parametrize("linear_solver", solvers)
 def test_solve_infeasible(linear_solver):
-    data = tools.gen_infeasible(K, n=m // 2)
+    rng = np.random.RandomState(3001)
+    data = tools.gen_infeasible(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=linear_solver, **params)
     sol = solver.solve()
     y = sol["y"]
@@ -80,7 +80,8 @@ def test_solve_infeasible(linear_solver):
 # TODO: indirect solver has trouble in this test, so disable for now
 @pytest.mark.parametrize("linear_solver", [scs.LinearSolver.QDLDL])
 def test_solve_unbounded(linear_solver):
-    data = tools.gen_unbounded(K, n=m // 2)
+    rng = np.random.RandomState(3002)
+    data = tools.gen_unbounded(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=linear_solver, **params)
     sol = solver.solve()
     x = sol["x"]
diff --git a/test/test_solve_random_cone_prob_accelerate.py b/test/test_solve_random_cone_prob_accelerate.py
index f6faf800..1c27a58b 100644
--- a/test/test_solve_random_cone_prob_accelerate.py
+++ b/test/test_solve_random_cone_prob_accelerate.py
@@ -17,8 +17,6 @@
 
 from scs import _scs_accelerate  # noqa: E402
 
-np.random.seed(1)
-
 # cone:
 K = {
     "z": 10,
@@ -34,7 +32,8 @@
 
 
 def test_solve_feasible():
-    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+    rng = np.random.RandomState(3000)
+    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.ACCELERATE, **params)
     sol = solver.solve()
     x = sol["x"]
@@ -53,7 +52,8 @@ def test_solve_feasible():
 
 
 def test_solve_infeasible():
-    data = tools.gen_infeasible(K, n=m // 2)
+    rng = np.random.RandomState(3001)
+    data = tools.gen_infeasible(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.ACCELERATE, **params)
     sol = solver.solve()
     y = sol["y"]
@@ -63,7 +63,8 @@ def test_solve_infeasible():
 
 
 def test_solve_unbounded():
-    data = tools.gen_unbounded(K, n=m // 2)
+    rng = np.random.RandomState(3002)
+    data = tools.gen_unbounded(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.ACCELERATE, **params)
     sol = solver.solve()
     x = sol["x"]
diff --git a/test/test_solve_random_cone_prob_cudss.py b/test/test_solve_random_cone_prob_cudss.py
index fcf280eb..664be696 100644
--- a/test/test_solve_random_cone_prob_cudss.py
+++ b/test/test_solve_random_cone_prob_cudss.py
@@ -21,8 +21,6 @@ def import_error(msg):
     import_error("Please install pytest to run tests.")
     raise
 
-np.random.seed(1)
-
 # cone:
 K = {
     "z": 10,
@@ -41,7 +39,8 @@ def import_error(msg):
     from scs import _scs_cudss
 
     def test_solve_feasible():
-        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+        rng = np.random.RandomState(3000)
+        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.CUDSS, **params)
         sol = solver.solve()
         x = sol["x"]
@@ -60,7 +59,8 @@ def test_solve_feasible():
         np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
     def test_solve_infeasible():
-        data = tools.gen_infeasible(K, n=m // 2)
+        rng = np.random.RandomState(3001)
+        data = tools.gen_infeasible(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.CUDSS, **params)
         sol = solver.solve()
         y = sol["y"]
@@ -69,7 +69,8 @@ def test_solve_infeasible():
         np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
     def test_solve_unbounded():
-        data = tools.gen_unbounded(K, n=m // 2)
+        rng = np.random.RandomState(3002)
+        data = tools.gen_unbounded(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.CUDSS, **params)
         sol = solver.solve()
         x = sol["x"]
diff --git a/test/test_solve_random_cone_prob_dense.py b/test/test_solve_random_cone_prob_dense.py
index 327d1a09..c95d1ef3 100644
--- a/test/test_solve_random_cone_prob_dense.py
+++ b/test/test_solve_random_cone_prob_dense.py
@@ -21,8 +21,6 @@ def import_error(msg):
     import_error("Please install pytest to run tests.")
     raise
 
-np.random.seed(1)
-
 # cone:
 K = {
     "z": 10,
@@ -41,7 +39,8 @@ def import_error(msg):
     from scs import _scs_dense
 
     def test_solve_feasible():
-        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+        rng = np.random.RandomState(3000)
+        data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.DENSE, **params)
         sol = solver.solve()
         x = sol["x"]
@@ -60,7 +59,8 @@ def test_solve_feasible():
         np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
     def test_solve_infeasible():
-        data = tools.gen_infeasible(K, n=m // 2)
+        rng = np.random.RandomState(3001)
+        data = tools.gen_infeasible(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.DENSE, **params)
         sol = solver.solve()
         y = sol["y"]
@@ -69,7 +69,8 @@ def test_solve_infeasible():
         np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
     def test_solve_unbounded():
-        data = tools.gen_unbounded(K, n=m // 2)
+        rng = np.random.RandomState(3002)
+        data = tools.gen_unbounded(K, n=m // 2, rng=rng)
         solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.DENSE, **params)
         sol = solver.solve()
         x = sol["x"]
diff --git a/test/test_solve_random_cone_prob_mkl.py b/test/test_solve_random_cone_prob_mkl.py
index 6aef9033..efb343ed 100644
--- a/test/test_solve_random_cone_prob_mkl.py
+++ b/test/test_solve_random_cone_prob_mkl.py
@@ -27,8 +27,6 @@
     # musllinux x86_64 ships openblas, not MKL
     pytest.skip("MKL module not installed", allow_module_level=True)
 
-np.random.seed(1)
-
 # cone:
 K = {
     "z": 10,
@@ -53,7 +51,8 @@ def test_mkl_module_links_intel_openmp():
 
 
 def test_solve_feasible():
-    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
+    rng = np.random.RandomState(3000)
+    data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.MKL, **params)
     sol = solver.solve()
     x = sol["x"]
@@ -72,7 +71,8 @@ def test_solve_feasible():
 
 
 def test_solve_infeasible():
-    data = tools.gen_infeasible(K, n=m // 2)
+    rng = np.random.RandomState(3001)
+    data = tools.gen_infeasible(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.MKL, **params)
     sol = solver.solve()
     y = sol["y"]
@@ -82,7 +82,8 @@ def test_solve_infeasible():
 
 
 def test_solve_unbounded():
-    data = tools.gen_unbounded(K, n=m // 2)
+    rng = np.random.RandomState(3002)
+    data = tools.gen_unbounded(K, n=m // 2, rng=rng)
     solver = scs.SCS(data, K, linear_solver=scs.LinearSolver.MKL, **params)
     sol = solver.solve()
     x = sol["x"]
diff --git a/test/test_spectral_and_complex_cones.py b/test/test_spectral_and_complex_cones.py
index 92a13fd7..0b73725c 100644
--- a/test/test_spectral_and_complex_cones.py
+++ b/test/test_spectral_and_complex_cones.py
@@ -51,20 +51,22 @@ def _cone_dims(cone):
     return m
 
 
-def _gen_feasible_qp(cone, n=None, p_scale=1.0):
+def _gen_feasible_qp(cone, n=None, p_scale=1.0, rng=None):
     """Generate a random feasible QP over the given cone.
 
     Uses P = p_scale * I (strong regularisation) and dense-ish A to
     guarantee boundedness and feasibility.
     """
+    if rng is None:
+        rng = np.random
     m = _cone_dims(cone)
     if n is None:
         n = m
     P = p_scale * sp.eye(n, format="csc")
-    A = sp.random(m, n, density=0.5, format="csc")
-    A.data = np.random.randn(A.nnz)
-    c = np.random.randn(n)
-    b = A @ np.random.randn(n) + np.abs(np.random.randn(m))
+    A = sp.random(m, n, density=0.5, format="csc", random_state=rng)
+    A.data = rng.randn(A.nnz)
+    c = rng.randn(n)
+    b = A @ rng.randn(n) + np.abs(rng.randn(m))
     return dict(P=P, A=A, b=b, c=c)
 
 
@@ -120,25 +122,25 @@ class TestComplexPSDCone:
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_complex_psd_feasibility(self, solver_opts):
-        np.random.seed(42)
+        rng = np.random.RandomState(42)
         cone = {"cs": [3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_complex_psd_multiple(self, solver_opts):
-        np.random.seed(123)
+        rng = np.random.RandomState(123)
         cone = {"cs": [2, 3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_mixed_real_complex_psd(self, solver_opts):
-        np.random.seed(456)
+        rng = np.random.RandomState(456)
         cone = dict(z=1, l=2, s=[3], cs=[3])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
@@ -153,33 +155,33 @@ class TestEll1Cone:
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_ell1_simple(self, solver_opts):
-        np.random.seed(10)
+        rng = np.random.RandomState(10)
         cone = {"ell1": [4]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_ell1_multiple(self, solver_opts):
-        np.random.seed(20)
+        rng = np.random.RandomState(20)
         cone = {"ell1": [3, 5]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_ell1_mixed_with_standard_cones(self, solver_opts):
-        np.random.seed(30)
+        rng = np.random.RandomState(30)
         cone = dict(z=1, l=2, q=[3], ell1=[4])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     def test_ell1_norm_bound(self):
         """Verify s in ell1 cone satisfies s[0] >= ||s[1:]||_1."""
-        np.random.seed(35)
+        rng = np.random.RandomState(35)
         cone = {"ell1": [5]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False, eps_abs=1e-9, eps_rel=1e-9)
         if sol["info"]["status_val"] == 1:
             s = sol["s"]
@@ -196,25 +198,25 @@ class TestNuclearNormCone:
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_nuclear_simple(self, solver_opts):
-        np.random.seed(40)
+        rng = np.random.RandomState(40)
         cone = {"nuc_m": [3], "nuc_n": [2]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_nuclear_multiple(self, solver_opts):
-        np.random.seed(50)
+        rng = np.random.RandomState(50)
         cone = {"nuc_m": [3, 4], "nuc_n": [2, 3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_nuclear_mixed(self, solver_opts):
-        np.random.seed(60)
+        rng = np.random.RandomState(60)
         cone = dict(l=3, nuc_m=[3], nuc_n=[2])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False)
         _assert_solved(sol)
 
@@ -227,10 +229,10 @@ def test_nuclear_mismatched_lengths_raises(self):
 
     def test_nuclear_norm_bound(self):
         """Verify s in nuclear norm cone satisfies ||X||_* <= t."""
-        np.random.seed(65)
+        rng = np.random.RandomState(65)
         rows, cols = 4, 3
         cone = {"nuc_m": [rows], "nuc_n": [cols]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False, eps_abs=1e-9, eps_rel=1e-9)
         if sol["info"]["status_val"] == 1:
             s = sol["s"]
@@ -251,25 +253,25 @@ class TestLogdetCone:
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_logdet_simple(self, solver_opts):
-        np.random.seed(70)
+        rng = np.random.RandomState(70)
         cone = {"d": [3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_logdet_multiple(self, solver_opts):
-        np.random.seed(80)
+        rng = np.random.RandomState(80)
         cone = {"d": [2, 3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_logdet_mixed(self, solver_opts):
-        np.random.seed(90)
+        rng = np.random.RandomState(90)
         cone = dict(l=2, d=[2])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
@@ -284,25 +286,25 @@ class TestSumLargestCone:
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_sum_largest_simple(self, solver_opts):
-        np.random.seed(100)
+        rng = np.random.RandomState(100)
         cone = {"sl_n": [4], "sl_k": [2]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_sum_largest_multiple(self, solver_opts):
-        np.random.seed(110)
+        rng = np.random.RandomState(110)
         cone = {"sl_n": [3, 4], "sl_k": [1, 2]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_sum_largest_mixed(self, solver_opts):
-        np.random.seed(120)
+        rng = np.random.RandomState(120)
         cone = dict(z=1, l=2, sl_n=[3], sl_k=[1])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
@@ -325,7 +327,7 @@ class TestAllConesCombined:
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_kitchen_sink(self, solver_opts):
         """Problem with every cone type present."""
-        np.random.seed(200)
+        rng = np.random.RandomState(200)
         cone = dict(
             z=1,
             l=2,
@@ -342,15 +344,15 @@ def test_kitchen_sink(self, solver_opts):
             sl_n=[3],
             sl_k=[1],
         )
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     @pytest.mark.parametrize("solver_opts", _solver_configs)
     def test_spectral_and_complex_psd(self, solver_opts):
-        np.random.seed(210)
+        rng = np.random.RandomState(210)
         cone = dict(cs=[3], ell1=[5], nuc_m=[4], nuc_n=[3])
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, **solver_opts, verbose=False, max_iters=5000)
         _assert_solved(sol)
 
@@ -364,37 +366,37 @@ def test_spectral_and_complex_psd(self, solver_opts):
 class TestEdgeCases:
 
     def test_ell1_dim_one(self):
-        np.random.seed(300)
+        rng = np.random.RandomState(300)
         cone = {"ell1": [1]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False)
         _assert_solved(sol)
 
     def test_nuclear_square(self):
-        np.random.seed(310)
+        rng = np.random.RandomState(310)
         cone = {"nuc_m": [3], "nuc_n": [3]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False)
         _assert_solved(sol)
 
     def test_logdet_dim_one(self):
-        np.random.seed(320)
+        rng = np.random.RandomState(320)
         cone = {"d": [1]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     def test_sum_largest_k_one(self):
-        np.random.seed(330)
+        rng = np.random.RandomState(330)
         cone = {"sl_n": [3], "sl_k": [1]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         sol = scs.solve(data, cone, verbose=False, max_iters=10000)
         _assert_solved(sol)
 
     def test_warm_start_with_spectral(self):
-        np.random.seed(340)
+        rng = np.random.RandomState(340)
         cone = {"ell1": [4]}
-        data = _gen_feasible_qp(cone)
+        data = _gen_feasible_qp(cone, rng=rng)
         solver = scs.SCS(data, cone, verbose=False, max_iters=500)
         sol1 = solver.solve()
         sol2 = solver.solve(warm_start=True)
diff --git a/test/test_thread_safety.py b/test/test_thread_safety.py
index 3f5c3f12..dfbc3966 100644
--- a/test/test_thread_safety.py
+++ b/test/test_thread_safety.py
@@ -8,6 +8,7 @@
 from concurrent.futures import ThreadPoolExecutor
 
 import numpy as np
+import pytest
 import scipy.sparse as sp
 from numpy.testing import assert_almost_equal
 
@@ -27,6 +28,7 @@ def _make_simple_lp():
     return data, cone, 1.0  # expected x[0]
 
 
+@pytest.mark.thread_unsafe(reason="creates its own threads internally")
 class TestThreadSafety:
     """Multiple threads call solve() on the same SCS object.
 
diff --git a/test/tsan-suppressions.txt b/test/tsan-suppressions.txt
new file mode 100644
index 00000000..9c6ad019
--- /dev/null
+++ b/test/tsan-suppressions.txt
@@ -0,0 +1,25 @@
+# Thread Sanitizer suppressions for scs-python
+#
+# These suppress known CPython-internal races that are not caused by scs.
+# Based on: https://github.com/python/cpython/blob/3.14/Tools/tsan/suppressions_free_threading.txt
+
+# CPython internal races (free-threaded build)
+race_top:assign_version_tag
+race_top:_Py_slot_tp_getattr_hook
+race_top:dump_traceback
+race_top:fatal_error
+race_top:_PyFrame_GetCode
+race_top:_PyFrame_Initialize
+race_top:_PyObject_TryGetInstanceAttribute
+race_top:PyUnstable_InterpreterFrame_GetLine
+race_top:type_modified_unlocked
+race_top:write_thread_id
+race_top:rangeiter_next
+race_top:mi_block_set_nextx
+race_top:update_one_slot
+race_top:set_tp_bases
+race_top:type_set_bases_unlocked
+thread:pthread_create
+race:list_ass_slice_lock_held
+race:list_inplace_repeat_lock_held
+race:PyObject_Realloc