diff --git a/cpu b/cpu
deleted file mode 160000
index 6714939..0000000
--- a/cpu
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit 671493983e42a77c75938edd622deaa7da6ceb6e
diff --git a/cpu b/cpu
new file mode 120000
index 0000000..d9a9826
--- /dev/null
+++ b/cpu
@@ -0,0 +1 @@
+../fastfields-cpu-lib
\ No newline at end of file
diff --git a/cuda b/cuda
deleted file mode 160000
index af3b6d0..0000000
--- a/cuda
+++ /dev/null
@@ -1 +0,0 @@
-Subproject commit af3b6d0bbf09a482b7ef3f24af2ea9ec4d36c99e
diff --git a/cuda b/cuda
new file mode 120000
index 0000000..6644e18
--- /dev/null
+++ b/cuda
@@ -0,0 +1 @@
+../fastfields-cuda-lib
\ No newline at end of file
diff --git a/reg_field.cpp b/reg_field.cpp
index 31dd33c..e2e4e98 100644
--- a/reg_field.cpp
+++ b/reg_field.cpp
@@ -80,4 +80,29 @@ void field_kernel(
     throw std::invalid_argument("unsupported device");
 }
 
+void field_relax(
+          DLTensor & sol       ,
+    const DLTensor & hes       ,
+    const DLTensor & grd       ,
+    const double   * voxel_size,
+    const double   * absolute  ,
+    const double   * membrane  ,
+    const double   * bending   ,
+          int8_t     bound     ,
+          int        ndim      ,
+          int        nb_iter   ,
+          int        stream    )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(sol))
+        return FF_CUDA::field_relax(sol, hes, grd, voxel_size, absolute, membrane, bending, bound, ndim, nb_iter, stream);
+#endif
+    if (IS_CPU(sol))
+        return FF_CPU::field_relax(sol, hes, grd, voxel_size, absolute, membrane, bending, bound, ndim, nb_iter, stream);
+
+    if (IS_CUDA(sol))
+        throw std::invalid_argument("fastfields: built without CUDA support, cannot operate on CUDA tensors");
+    throw std::invalid_argument("unsupported device");
+}
+
 FF_NAMESPACE_END(FF)
diff --git a/reg_field.h b/reg_field.h
index a16023b..67f5940 100644
--- a/reg_field.h
+++ b/reg_field.h
@@ -106,6 +106,34 @@ void field_kernel(
           int        stream    = 0
 );
 
+/**
+ * @brief In-place relaxation (Gauss-Seidel) sweeps solving `(H + L) x = g`.
+ *
+ * Refines the warm-started field `sol` towards the solution of the regularised
+ * system, where `H` is the per-voxel compact-symmetric Hessian (`hes`, packed
+ * `C*(C+1)/2` last axis), `L` the field regulariser (same per-channel penalties
+ * as `field_matvec`), and `g` the gradient (`grd`, `C` last axis). Runs
+ * `nb_iter` red-black sweeps and writes the refined solution back into `sol`.
+ *
+ * @param sol        Field to refine, in/out (*batch, *spatial, C)
+ * @param hes        Compact-symmetric Hessian (*batch, *spatial, C*(C+1)/2)
+ * @param grd        Gradient (*batch, *spatial, C)
+ * @param nb_iter    Number of relaxation iterations
+ */
+void field_relax(
+          DLTensor & sol       ,
+    const DLTensor & hes       ,
+    const DLTensor & grd       ,
+    const double   * voxel_size = nullptr,
+    const double   * absolute  = nullptr,
+    const double   * membrane  = nullptr,
+    const double   * bending   = nullptr,
+          int8_t     bound     = bound_t::DCT2,
+          int        ndim      = 1,
+          int        nb_iter   = 1,
+          int        stream    = 0
+);
+
 FF_NAMESPACE_END(FF)
 
 #endif // FF_LIB_REG_FIELD