tjayasinghe · tjayasinghe · Jun 30, 2026 · Jun 30, 2026 · Jun 30, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -4,6 +4,26 @@ All notable changes to cuPeriod are documented here. The format is based on
 [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to
 [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [Unreleased]
+
+### Added
+
+- **Multi-vendor GPU support via PyTorch and the Python array API.** GLS and BLS gain a
+  portable `torch` backend that runs on AMD (ROCm), Intel (XPU), and Apple (MPS) GPUs as
+  well as a real CPU path — so the accelerated code is no longer NVIDIA-only, and works
+  even with no GPU at all. Select it with `backend="torch"` (or `"torch:cpu"`,
+  `"torch:cuda"`, `"torch:mps"`, `"torch:xpu"`); `backend="auto"` now reaches a torch GPU
+  on non-NVIDIA machines after the cufinufft/cupy fast paths.
+  - GLS adds a NUFFT-free direct trig-sum path (the portable formulation; cufinufft
+    remains the NVIDIA fast path).
+  - BLS runs its vectorized box search through the array-API namespace (the cupy
+    `RawKernel` remains the NVIDIA fast path).
+  - New `device` and `precision` settings: `precision="auto"` is float64 everywhere it is
+    supported and float32 only where the device forces it (Apple MPS cannot do float64);
+    an explicit `precision="float64"` on MPS raises rather than silently downgrading.
+- `array-api-compat` is now a dependency; install the portable accelerator with the
+  `[torch]` extra (`pip install 'cuperiod[torch]'`).
+
 ## [1.0.0] — 2026-06-29
 
 First public release.

diff --git a/benchmarks/README.md b/benchmarks/README.md
@@ -18,6 +18,12 @@ The rendered results are in **[REPORT.md](REPORT.md)** with figures in `figures/
 | Period recovery | all 7 | VSX literature period | find the real period? |
 | Performance | all 7 | astropy / PyAstronomy | how much faster, and how does it scale? |
 
+The performance benchmark also times the **portable `torch` backend** (`backend="torch"`;
+the resolved device is shown in the `torch_backend` column) for the ported methods — GLS
+and BLS — alongside the CPU and CUDA paths, so the cross-vendor path (AMD/Intel/Mac/CPU) is
+tracked. A backend absent on the host (no CUDA GPU, or no torch) is recorded blank rather
+than failing the sweep.
+
 ## Data
 
 * **`dataset/light_curves.parquet`** — 72 real ASAS-SN g-band light curves

diff --git a/benchmarks/benchmark.py b/benchmarks/benchmark.py
@@ -26,6 +26,9 @@
 import cuperiod as cup  # noqa: E402
 
 from _common import RESULTS, load_dataset  # noqa: E402
+from cuperiod.core.backend import torch_available  # noqa: E402
+from cuperiod.core.errors import BackendUnavailableError  # noqa: E402
+from cuperiod.methods.base import get_method  # noqa: E402
 
 # bounded period window for the box methods (BLS/TLS) — independent of the star,
 # so the trial-period count can never blow up on a short-period target.
@@ -60,6 +63,23 @@ def best_time(fn, repeat=3):
     return min(_timed(fn) for _ in range(repeat))
 
 
+def safe_best_time(fn, repeat=3):
+    """:func:`best_time`, but a backend unavailable here records NaN instead of raising.
+
+    Lets the sweep run on machines missing a backend — no CUDA GPU (the ``gpu`` column),
+    or no torch / no torch GPU device (the ``torch`` column) — leaving that cell blank.
+    """
+    try:
+        return best_time(fn, repeat)
+    except BackendUnavailableError:
+        return float("nan")
+
+
+def supports_torch(method):
+    """Whether ``method`` has the portable torch backend and torch is importable."""
+    return torch_available() and "torch" in get_method(method).all_backends
+
+
 def _timed(fn):
     t0 = time.perf_counter(); fn(); return time.perf_counter() - t0
 
@@ -91,13 +111,23 @@ def bench_single(t, y, e):
     for m, st in FREQ_SETTINGS.items():
         be = cup.periodogram((t, y, e), m, backend="cpu", grid=grid, settings=st()).backend
         tc = best_time(lambda: cup.periodogram((t, y, e), m, backend="cpu", grid=grid, settings=st()))
-        tg = best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", grid=grid, settings=st()))
+        tg = safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", grid=grid, settings=st()))
+        has_torch = supports_torch(m)
+        tt = (safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="torch", grid=grid, settings=st()))
+              if has_torch else np.nan)
+        tbe = (cup.periodogram((t, y, e), m, backend="torch", grid=grid, settings=st()).backend
+               if has_torch and np.isfinite(tt) else "—")
         tr = best_time(reftime[m], repeat=1) if m in reftime else np.nan
         rows.append(dict(method=m, n_grid=grid.size, cpu_backend=be, cpu_s=tc, gpu_s=tg,
-                         ref_s=tr, ref=refname.get(m, "—"), gpu_speedup=tc / tg,
-                         gpu_vs_ref=(tr / tg if np.isfinite(tr) else np.nan),
+                         torch_s=tt, torch_backend=tbe,
+                         ref_s=tr, ref=refname.get(m, "—"),
+                         gpu_speedup=(tc / tg if np.isfinite(tg) else np.nan),
+                         torch_speedup=(tc / tt if np.isfinite(tt) else np.nan),
+                         gpu_vs_ref=(tr / tg if np.isfinite(tr) and np.isfinite(tg) else np.nan),
                          cpu_vs_ref=(tr / tc if np.isfinite(tr) else np.nan)))
-        print(f"  {m:12s} cpu({be})={tc:.3f}s gpu={tg:.4f}s ({tc/tg:.0f}x)", flush=True)
+        gstr = f"gpu={tg:.4f}s ({tc/tg:.0f}x)" if np.isfinite(tg) else "gpu=—"
+        tstr = f" torch({tbe})={tt:.3f}s" if np.isfinite(tt) else ""
+        print(f"  {m:12s} cpu({be})={tc:.3f}s {gstr}{tstr}", flush=True)
     # box methods on a fixed, bounded period window. cpu_s is cuPeriod's *default*
     # CPU backend: "cpu" resolves to the multicore numba box search for BLS (or
     # astropy if numba is not installed), to numpy for TLS.
@@ -107,8 +137,13 @@ def bench_single(t, y, e):
              cup.TLSSettings(min_period_days=BOX_PMIN, max_period_days=BOX_PMAX))
         be = cup.periodogram((t, y, e), m, backend="cpu", settings=S).backend
         tc = best_time(lambda: cup.periodogram((t, y, e), m, backend="cpu", settings=S), repeat=2)
-        tg = best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", settings=S), repeat=2)
-        ng = cup.periodogram((t, y, e), m, backend="gpu", settings=S).power.size
+        tg = safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", settings=S), repeat=2)
+        ng = cup.periodogram((t, y, e), m, backend="cpu", settings=S).power.size  # backend-independent
+        has_torch = supports_torch(m)
+        tt = (safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="torch", settings=S), repeat=2)
+              if has_torch else np.nan)
+        tbe = (cup.periodogram((t, y, e), m, backend="torch", settings=S).backend
+               if has_torch and np.isfinite(tt) else "—")
         # BLS reference = astropy's compiled BoxLeastSquares; also time the pure-numpy
         # GPU-parity reference port to document it is not the product path.
         ref_s = (best_time(lambda: cup.periodogram((t, y, e), "BLS", backend="astropy", settings=S), repeat=1)
@@ -117,13 +152,18 @@ def bench_single(t, y, e):
         port_s = (best_time(lambda: cup.periodogram((t, y, e), "BLS", backend="numpy", settings=S), repeat=1)
                   if m == "BLS" else np.nan)
         rows.append(dict(method=m, n_grid=int(ng), cpu_backend=be, cpu_s=tc, gpu_s=tg,
-                         ref_s=ref_s, ref=ref, gpu_speedup=tc / tg,
-                         gpu_vs_ref=(ref_s / tg if np.isfinite(ref_s) else np.nan),
+                         torch_s=tt, torch_backend=tbe,
+                         ref_s=ref_s, ref=ref,
+                         gpu_speedup=(tc / tg if np.isfinite(tg) else np.nan),
+                         torch_speedup=(tc / tt if np.isfinite(tt) else np.nan),
+                         gpu_vs_ref=(ref_s / tg if np.isfinite(ref_s) and np.isfinite(tg) else np.nan),
                          cpu_vs_ref=(ref_s / tc if np.isfinite(ref_s) else np.nan),
                          cpu_port_s=port_s))
+        gstr = f"gpu={tg:.4f}s (gpu {tc/tg:.0f}x)" if np.isfinite(tg) else "gpu=—"
+        tstr = f" torch({tbe})={tt:.3f}s" if np.isfinite(tt) else ""
         extra = (f"  [vs astropy {ref_s/tc:.0f}x faster; numpy-port {port_s:.1f}s]"
                  if m == "BLS" else "")
-        print(f"  {m:12s} cpu({be})={tc:.3f}s gpu={tg:.4f}s (gpu {tc/tg:.0f}x){extra}  [{ng:,} periods]",
+        print(f"  {m:12s} cpu({be})={tc:.3f}s {gstr}{tstr}{extra}  [{ng:,} periods]",
               flush=True)
     df = pd.DataFrame(rows)
     df.to_parquet(RESULTS / "bench_single.parquet", index=False)
@@ -138,8 +178,12 @@ def bench_scaling_npoints(t, y, e):
         for m in ["GLS", "PDM", "MHAOV"]:
             st = FREQ_SETTINGS[m]
             tc = best_time(lambda: cup.periodogram((tt, yy, ee), m, backend="cpu", grid=grid, settings=st()), repeat=1)
-            tg = best_time(lambda: cup.periodogram((tt, yy, ee), m, backend="gpu", grid=grid, settings=st()), repeat=1)
-            rows.append(dict(axis="npoints", method=m, n=n, cpu_s=tc, gpu_s=tg, speedup=tc / tg))
+            tg = safe_best_time(lambda: cup.periodogram((tt, yy, ee), m, backend="gpu", grid=grid, settings=st()), repeat=1)
+            ttor = (safe_best_time(lambda: cup.periodogram((tt, yy, ee), m, backend="torch", grid=grid, settings=st()), repeat=1)
+                    if supports_torch(m) else np.nan)
+            rows.append(dict(axis="npoints", method=m, n=n, cpu_s=tc, gpu_s=tg, torch_s=ttor,
+                             speedup=(tc / tg if np.isfinite(tg) else np.nan),
+                             torch_speedup=(tc / ttor if np.isfinite(ttor) else np.nan)))
         print(f"  N={n:>6}: done", flush=True)
     df = pd.DataFrame(rows)
     df.to_parquet(RESULTS / "bench_npoints.parquet", index=False)
@@ -153,8 +197,12 @@ def bench_scaling_grid(t, y, e):
         for m in ["GLS", "PDM", "MHAOV"]:
             st = FREQ_SETTINGS[m]
             tc = best_time(lambda: cup.periodogram((t, y, e), m, backend="cpu", grid=grid, settings=st()), repeat=1)
-            tg = best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", grid=grid, settings=st()), repeat=1)
-            rows.append(dict(axis="grid", method=m, n=n, cpu_s=tc, gpu_s=tg, speedup=tc / tg))
+            tg = safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="gpu", grid=grid, settings=st()), repeat=1)
+            ttor = (safe_best_time(lambda: cup.periodogram((t, y, e), m, backend="torch", grid=grid, settings=st()), repeat=1)
+                    if supports_torch(m) else np.nan)
+            rows.append(dict(axis="grid", method=m, n=n, cpu_s=tc, gpu_s=tg, torch_s=ttor,
+                             speedup=(tc / tg if np.isfinite(tg) else np.nan),
+                             torch_speedup=(tc / ttor if np.isfinite(ttor) else np.nan)))
         print(f"  grid={n:>7}: done", flush=True)
     df = pd.DataFrame(rows)
     df.to_parquet(RESULTS / "bench_grid.parquet", index=False)

diff --git a/benchmarks/make_report.py b/benchmarks/make_report.py
@@ -432,21 +432,27 @@ def main():
                      f"({b.cpu_s*1e3:.0f} ms vs {b.ref_s:.1f} s on this light curve), matching it "
                      f"to floating-point{par}. The GPU then adds another {b.gpu_speedup:.0f}× "
                      f"({b.ref_s/b.gpu_s:.0f}× over astropy).\n")
-        cols = ["method", "cpu_backend", "cpu_s", "gpu_s", "ref", "ref_s", "cpu_vs_ref", "gpu_speedup"]
+        cols = ["method", "cpu_backend", "cpu_s", "gpu_s", "torch_s", "torch_backend",
+                "ref", "ref_s", "cpu_vs_ref", "gpu_speedup"]
         cols = [c for c in cols if c in s.columns]
+        nan_dash = lambda fmt: (lambda v: ("—" if not np.isfinite(v) else fmt(v)))
         L.append(md_table(s, cols, {
-            "cpu_s": lambda v: f"{v:.3f}", "gpu_s": lambda v: f"{v:.4f}",
-            "gpu_speedup": lambda v: f"{v:.0f}x",
-            "cpu_vs_ref": lambda v: ("—" if not np.isfinite(v) else f"{v:.0f}x"),
-            "ref_s": lambda v: ("—" if not np.isfinite(v) else f"{v:.2f}"),
+            "cpu_s": lambda v: f"{v:.3f}",
+            "gpu_s": nan_dash(lambda v: f"{v:.4f}"),
+            "torch_s": nan_dash(lambda v: f"{v:.3f}"),
+            "gpu_speedup": nan_dash(lambda v: f"{v:.0f}x"),
+            "cpu_vs_ref": nan_dash(lambda v: f"{v:.0f}x"),
+            "ref_s": nan_dash(lambda v: f"{v:.2f}"),
             "method": ml}))
         L.append("\n*cpu_backend* = what `backend=\"cpu\"` resolves to — the fast default a user "
                  "gets: finufft (GLS), the multicore numba box search (BLS), numpy (the rest). "
                  "*ref* = the established external tool; *cpu_vs_ref* = how much faster cuPeriod's "
-                 "CPU is than that tool; *gpu_speedup* = GPU over cuPeriod's CPU. cuPeriod's CPU "
-                 "path already beats every reference tool it has (GLS, PDM, BLS) — so the GPU's "
-                 "marginal gain is small where the CPU is already fast (BLS, GLS) and large where "
-                 "it is not (PDM, MHAOV, TLS).\n")
+                 "CPU is than that tool; *gpu_speedup* = GPU over cuPeriod's CPU. *torch_s* = the "
+                 "portable PyTorch backend (device shown in *torch_backend*: cpu/cuda/mps/xpu) — "
+                 "the cross-vendor path that also runs on AMD/Intel/Mac; blank for methods not yet "
+                 "ported to it. cuPeriod's CPU path already beats every reference tool it has "
+                 "(GLS, PDM, BLS) — so the GPU's marginal gain is small where the CPU is already "
+                 "fast (BLS, GLS) and large where it is not (PDM, MHAOV, TLS).\n")
         if len(bls) and "cpu_port_s" in bls and np.isfinite(bls.cpu_port_s.iloc[0]):
             b = bls.iloc[0]
             L.append(f"\n> The pure-`numpy` BLS backend shares one array-module-generic source "

diff --git a/pyproject.toml b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "cuperiod"
-version = "1.0.0"
+version = "1.1.0.dev0"
 description = "Optimized, GPU-accelerated periodograms for astronomy"
 readme = "README.md"
 requires-python = ">=3.11"
@@ -19,6 +19,12 @@ keywords = [
     "time-series",
     "gpu",
     "cuda",
+    "rocm",
+    "amd",
+    "intel",
+    "metal",
+    "array-api",
+    "pytorch",
 ]
 classifiers = [
     "Development Status :: 5 - Production/Stable",
@@ -36,6 +42,9 @@ dependencies = [
     "scipy>=1.10",
     "astropy>=6.0",
     "finufft>=2.2",
+    # Array-API dispatch: the portable (numpy/cupy/torch) compute paths run through one
+    # standard namespace. Lightweight, pure-Python — always installed.
+    "array-api-compat>=1.9",
     "pydantic>=2.5",
     "pydantic-settings>=2.1",
     "typer>=0.12",
@@ -59,6 +68,11 @@ pandas = ["pandas>=2.0"]
 # astropy's compiled BoxLeastSquares by ~20x and is auto-selected on the CPU when
 # present (otherwise BLS falls back to astropy).
 fast = ["numba>=0.59"]
+# Portable accelerator reaching AMD (ROCm), Intel (XPU), Apple (MPS), and a fast CPU
+# path. Install the wheel matching your accelerator from https://pytorch.org/ — the
+# plain wheel is CPU-only; CUDA/ROCm/XPU builds come from PyTorch's own index, so we do
+# not pin a hardware-specific build here.
+torch = ["torch>=2.2"]
 # Documentation toolchain (Sphinx + Furo theme). pandas is included so the
 # DataFrame-ingestion examples and the autodoc of from_dataframe build cleanly.
 docs = [
@@ -77,6 +91,8 @@ dev = [
     "pandas>=2.0",
     "hypothesis>=6",
     "numba>=0.59",
+    # Exercise the portable torch-CPU path in CI (GPU devices auto-skip).
+    "torch>=2.2",
 ]
 
 [project.scripts]
@@ -112,6 +128,9 @@ src = ["src", "tests"]
 
 [tool.ruff.lint]
 select = ["E", "F", "I", "UP", "B", "SIM", "NPY"]
+# UP038 (use `X | Y` in isinstance) is deprecated upstream: the tuple form
+# `isinstance(x, (A, B))` is faster at runtime, so we keep it.
+ignore = ["UP038"]
 
 [tool.ruff.lint.per-file-ignores]
 # typer's API requires function calls (Option/Argument) in parameter defaults.
@@ -140,6 +159,10 @@ module = [
     "astropy.*",
     "scipy",
     "scipy.*",
+    "torch",
+    "torch.*",
+    "array_api_compat",
+    "array_api_compat.*",
 ]
 ignore_missing_imports = true