From 596fc3f3be4f1caaf9b5584a7c6e9b0a7410dc14 Mon Sep 17 00:00:00 2001
From: Andrew Reisner <areisner@lanl.gov>
Date: Tue, 11 Aug 2026 12:04:09 -0600
Subject: [PATCH 1/4] Begin adding support for futures as scalars in vector
 operations.

This begins adding preliminary support for scalars in vector
operations.  The basic machinery is added along with a proof of
concept with axpy.
---
 flecsolve/util/CMakeLists.txt                |   1 +
 flecsolve/util/future.hh                     |  27 ++-
 flecsolve/util/scalar_ops.hh                 |  58 +++++
 flecsolve/vectors/core.hh                    |  23 +-
 flecsolve/vectors/operations/multi.hh        |  16 ++
 flecsolve/vectors/operations/topo_tasks.hh   | 219 +++++++++++++------
 flecsolve/vectors/operations/topo_view.hh    | 105 ++++++---
 flecsolve/vectors/test/flecsi_multivector.cc |   3 +
 flecsolve/vectors/test/flecsi_vector.cc      | 132 +++++++----
 9 files changed, 437 insertions(+), 147 deletions(-)
 create mode 100644 flecsolve/util/scalar_ops.hh

diff --git a/flecsolve/util/CMakeLists.txt b/flecsolve/util/CMakeLists.txt
index beffa42..2b3641c 100644
--- a/flecsolve/util/CMakeLists.txt
+++ b/flecsolve/util/CMakeLists.txt
@@ -2,4 +2,5 @@ flecsolve_headers(
   config.hh
   future.hh
   traits.hh
+  scalar_ops.hh
   test/mesh.hh)
diff --git a/flecsolve/util/future.hh b/flecsolve/util/future.hh
index 7620466..53560c6 100644
--- a/flecsolve/util/future.hh
+++ b/flecsolve/util/future.hh
@@ -15,11 +15,16 @@ to do so.
 */
 #pragma once
 
+#include <flecsi/execution.hh>
+
 #include <tuple>
 #include <functional>
-
 namespace flecsolve {
 
+// for now, only consider scalar futures
+template<class T>
+using future = flecsi::future<T, flecsi::exec::launch_type_t::single>;
+
 template<class... Futures>
 struct future_vector {
 	decltype(auto) get() {
@@ -57,4 +62,24 @@ struct future_transform {
 template<class Future, class F>
 future_transform(Future, F)
 	-> future_transform<Future, F>; // automatic in C++20
+
+template<class T, class Scalar>
+struct is_scalar_or_future_impl : std::false_type {
+};
+
+template<class T, class Scalar>
+struct is_scalar_or_future_impl<future<T>, Scalar>
+	: std::bool_constant<std::is_convertible_v<T, Scalar>> {
+};
+
+template<class T, class F, class Scalar>
+struct is_scalar_or_future_impl<future_transform<future<T>, F>, Scalar>
+	: std::bool_constant<std::is_convertible_v<T, Scalar>> {
+};
+
+template<class T, class Scalar>
+inline constexpr bool is_scalar_or_future_v =
+	std::is_convertible_v<T, Scalar> ||
+	is_scalar_or_future_impl<T, Scalar>::value;
+
 }
diff --git a/flecsolve/util/scalar_ops.hh b/flecsolve/util/scalar_ops.hh
new file mode 100644
index 0000000..55229d1
--- /dev/null
+++ b/flecsolve/util/scalar_ops.hh
@@ -0,0 +1,58 @@
+/*
+© 2025. Triad National Security, LLC. All rights reserved.
+
+This program was produced under U.S. Government contract
+89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+operated by Triad National Security, LLC for the U.S. Department of
+Energy/National Nuclear Security Administration. All rights in the
+program are reserved by Triad National Security, LLC, and the U.S.
+Department of Energy/National Nuclear Security Administration. The
+Government is granted for itself and others acting on its behalf a
+nonexclusive, paid-up, irrevocable worldwide license in this material
+to reproduce, prepare. derivative works, distribute copies to the
+public, perform publicly and display publicly, and to permit others
+to do so.
+*/
+#ifndef FLECSI_LINALG_UTIL_SCALAR_OPS_HH
+#define FLECSI_LINALG_UTIL_SCALAR_OPS_HH
+
+#include <flecsi/execution.hh>
+
+namespace flecsolve::scalar_ops {
+
+template<class T>
+struct axpy {
+	constexpr T operator()(T x) const noexcept { return a * x + y; }
+	T a{0}, y{0};
+};
+
+template<class T>
+struct identity {
+	constexpr T operator()(T v) const noexcept { return v; }
+};
+
+template<class F, class G>
+struct compose {
+	template<class T>
+	constexpr auto operator()(T v) const noexcept {
+		return g(f(v));
+	}
+	F f;
+	G g;
+};
+}
+
+namespace flecsolve {
+template<class T>
+auto operator*(T v, flecsi::future<T> f) {
+	return future_transform{std::move(f), scalar_ops::axpy<T>{v, 0}};
+}
+
+template<class T>
+auto operator*(flecsi::future<T> f, T v) {
+	return future_transform{std::move(f), scalar_ops::axpy<T>{v, 0}};
+}
+
+}
+
+#endif
diff --git a/flecsolve/vectors/core.hh b/flecsolve/vectors/core.hh
index f1bc561..6fb3448 100644
--- a/flecsolve/vectors/core.hh
+++ b/flecsolve/vectors/core.hh
@@ -19,6 +19,9 @@ to do so.
 #include <random>
 
 #include "flecsolve/util/traits.hh"
+#include "flecsolve/util/future.hh"
+#include "flecsolve/util/scalar_ops.hh"
+
 #include "variable.hh"
 
 namespace flecsolve::vec {
@@ -182,12 +185,15 @@ struct core : Config {
 	 *
 	 * \f$\mathit{this}_i = alpha x_i + y_i\f$
 	 */
-	template<class V1,
+	template<class Alpha,
+	         class V1,
 	         class V2,
 	         std::enable_if_t<is_vector_v<V1>, bool> = true,
 	         std::enable_if_t<is_vector_v<V2>, bool> = true>
-	void axpy(scalar alpha, const V1 & x, const V2 & y) {
-		ops::axpy(alpha, x.data, y.data, data);
+	void axpy(Alpha alpha, const V1 & x, const V2 & y) {
+		static_assert(is_scalar_or_future_v<Alpha, scalar>,
+		              "axpy: alpha must be scalar or future");
+		ops::axpy(defer(alpha), x.data, y.data, data);
 	}
 
 	/**
@@ -350,6 +356,17 @@ struct core : Config {
 	}
 
 	data_t data;
+
+private:
+	scalar defer(scalar v) { return v; }
+	template<class T>
+	auto defer(future<T> f) {
+		return future_transform{std::move(f), scalar_ops::identity<T>{}};
+	}
+	template<class T, class F>
+	auto defer(future_transform<future<T>, F> f) {
+		return f;
+	}
 };
 
 template<template<class> class Data, template<class> class Ops, class Config>
diff --git a/flecsolve/vectors/operations/multi.hh b/flecsolve/vectors/operations/multi.hh
index 5c32468..409dc9a 100644
--- a/flecsolve/vectors/operations/multi.hh
+++ b/flecsolve/vectors/operations/multi.hh
@@ -135,6 +135,22 @@ struct multi {
 		      y.components);
 	}
 
+	template<class F, class T, class T0, class T1>
+	static void axpy(future_transform<future<T>, F> alpha,
+	                 const T0 & x,
+	                 const T1 & y,
+	                 vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "axpy: future type must be convertible to scalar");
+		apply([&alpha](auto & z,
+		               const auto & x,
+		               const auto & y) { z.axpy(alpha, x, y); },
+		      make_is(),
+		      z.components,
+		      x.components,
+		      y.components);
+	}
+
 	template<class T>
 	static void axpby(scalar alpha, scalar beta, const T & x, vec_data & z) {
 		apply([alpha, beta](auto & z,
diff --git a/flecsolve/vectors/operations/topo_tasks.hh b/flecsolve/vectors/operations/topo_tasks.hh
index d9faf99..17085a4 100644
--- a/flecsolve/vectors/operations/topo_tasks.hh
+++ b/flecsolve/vectors/operations/topo_tasks.hh
@@ -49,13 +49,15 @@ struct topo_tasks {
 
 	using util = typename VecData::util;
 
-	static scalar scalar_prod(flecsi::exec::accelerator s, topo_acc m, acc<ro> x, acc<ro> y) noexcept {
-		scalar res = s.executor().named("scalar_prod")
-			.reduceall(dof,
-			           up,
-			           util::dofs(m),
-			           flecsi::exec::fold::sum,
-			           scalar) {
+	static scalar scalar_prod(flecsi::exec::accelerator s,
+	                          topo_acc m,
+	                          acc<ro> x,
+	                          acc<ro> y) noexcept {
+		scalar res =
+			s.executor()
+				.named("scalar_prod")
+				.reduceall(
+					dof, up, util::dofs(m), flecsi::exec::fold::sum, scalar) {
 			if constexpr (is_complex)
 				up(std::conj(y[dof]) * x[dof]);
 			else
@@ -65,38 +67,65 @@ struct topo_tasks {
 		return res;
 	}
 
-	static void set_to_scalar(flecsi::exec::accelerator s, topo_acc m, acc<wo> x, scalar val) noexcept {
+	static void set_to_scalar(flecsi::exec::accelerator s,
+	                          topo_acc m,
+	                          acc<wo> x,
+	                          scalar val) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = val; };
 	}
 
-	static void scale_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x, scalar val) noexcept {
+	static void scale_self(flecsi::exec::accelerator s,
+	                       topo_acc m,
+	                       acc<rw> x,
+	                       scalar val) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] *= val; };
 	}
 
-	static void scale(flecsi::exec::accelerator s, topo_acc m, acc<ro> x, acc<wo> y, scalar val) noexcept {
+	static void scale(flecsi::exec::accelerator s,
+	                  topo_acc m,
+	                  acc<ro> x,
+	                  acc<wo> y,
+	                  scalar val) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { y[dof] = x[dof] * val; };
 	}
 
 	template<class OtherAcc>
-	static void copy(flecsi::exec::accelerator s, topo_acc, acc_all<wo> xa, OtherAcc ya) noexcept {
+	static void copy(flecsi::exec::accelerator s,
+	                 topo_acc,
+	                 acc_all<wo> xa,
+	                 OtherAcc ya) noexcept {
 		flecsi::util::iota_view<std::size_t> v(0, xa.span().size());
 		s.executor().forall(i, v) { xa[i] = ya[i]; };
 	}
 
-	static void add_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> z, acc<ro> x) noexcept {
+	static void add_self(flecsi::exec::accelerator s,
+	                     topo_acc m,
+	                     acc<rw> z,
+	                     acc<ro> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { z[dof] = z[dof] + x[dof]; };
 	}
 
-	static void add(flecsi::exec::accelerator s, topo_acc m, acc<wo> z, acc<ro> x, acc<ro> y) noexcept {
+	static void add(flecsi::exec::accelerator s,
+	                topo_acc m,
+	                acc<wo> z,
+	                acc<ro> x,
+	                acc<ro> y) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { z[dof] = x[dof] + y[dof]; };
 	}
 
-	static void subtract(flecsi::exec::accelerator s, topo_acc m, acc<wo> x, acc<ro> a, acc<ro> b) noexcept {
+	static void subtract(flecsi::exec::accelerator s,
+	                     topo_acc m,
+	                     acc<wo> x,
+	                     acc<ro> a,
+	                     acc<ro> b) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = a[dof] - b[dof]; };
 	}
 
 	template<bool inv>
-	static void subtract_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x, acc<ro> b) noexcept {
+	static void subtract_self(flecsi::exec::accelerator s,
+	                          topo_acc m,
+	                          acc<rw> x,
+	                          acc<ro> b) noexcept {
 		s.executor().forall(dof, util::dofs(m)) {
 			if constexpr (inv) {
 				x[dof] = x[dof] - b[dof];
@@ -107,18 +136,26 @@ struct topo_tasks {
 		};
 	}
 
-	static void multiply(flecsi::exec::accelerator s, topo_acc m, acc<wo> z, acc<ro> x, acc<ro> y) noexcept {
+	static void multiply(flecsi::exec::accelerator s,
+	                     topo_acc m,
+	                     acc<wo> z,
+	                     acc<ro> x,
+	                     acc<ro> y) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { z[dof] = x[dof] * y[dof]; };
 	}
 
-	static void multiply_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x, acc<ro> y) noexcept {
-		s.executor().forall(dof, util::dofs(m)) {
-			x[dof] = x[dof] * y[dof];
-		};
+	static void multiply_self(flecsi::exec::accelerator s,
+	                          topo_acc m,
+	                          acc<rw> x,
+	                          acc<ro> y) noexcept {
+		s.executor().forall(dof, util::dofs(m)) { x[dof] = x[dof] * y[dof]; };
 	}
 
 	template<bool inv>
-	static void divide_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> z, acc<ro> x) noexcept {
+	static void divide_self(flecsi::exec::accelerator s,
+	                        topo_acc m,
+	                        acc<rw> z,
+	                        acc<ro> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) {
 			if constexpr (inv) {
 				z[dof] = z[dof] / x[dof];
@@ -129,15 +166,24 @@ struct topo_tasks {
 		};
 	}
 
-	static void divide(flecsi::exec::accelerator s, topo_acc m, acc<wo> z, acc<ro> x, acc<ro> y) noexcept {
+	static void divide(flecsi::exec::accelerator s,
+	                   topo_acc m,
+	                   acc<wo> z,
+	                   acc<ro> x,
+	                   acc<ro> y) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { z[dof] = x[dof] / y[dof]; };
 	}
 
-	static void reciprocal_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x) noexcept {
+	static void reciprocal_self(flecsi::exec::accelerator s,
+	                            topo_acc m,
+	                            acc<rw> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = 1.0 / x[dof]; };
 	}
 
-	static void reciprocal(flecsi::exec::accelerator s, topo_acc m, acc<wo> x, acc<ro> y) noexcept {
+	static void reciprocal(flecsi::exec::accelerator s,
+	                       topo_acc m,
+	                       acc<wo> x,
+	                       acc<ro> y) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = 1.0 / y[dof]; };
 	}
 
@@ -171,15 +217,35 @@ struct topo_tasks {
 		};
 	}
 
-	static void
-	axpy(flecsi::exec::accelerator s, topo_acc m, acc<wo> z, scalar alpha, acc<ro> x, acc<ro> y) noexcept {
+	static void axpy(flecsi::exec::accelerator s,
+	                 topo_acc m,
+	                 acc<wo> z,
+	                 scalar alpha,
+	                 acc<ro> x,
+	                 acc<ro> y) noexcept {
 		s.executor().forall(dof, util::dofs(m)) {
 			z[dof] = alpha * x[dof] + y[dof];
 		};
 	}
 
+	template<class F, class T>
+	static void axpy_future(flecsi::exec::accelerator s,
+	                        topo_acc m,
+	                        acc<wo> z,
+	                        future<T> alpha_fut,
+	                        F f,
+	                        acc<ro> x,
+	                        acc<ro> y) noexcept {
+		scalar alpha = f(alpha_fut.get());
+		axpy(s, m, z, alpha, x, y);
+	}
+
 	template<bool inv>
-	static void axpy_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> z, acc<ro> x, scalar alpha) noexcept {
+	static void axpy_self(flecsi::exec::accelerator s,
+	                      topo_acc m,
+	                      acc<rw> z,
+	                      acc<ro> x,
+	                      scalar alpha) noexcept {
 		s.executor().forall(dof, util::dofs(m)) {
 			if constexpr (inv) {
 				z[dof] = alpha * z[dof] + x[dof];
@@ -190,42 +256,69 @@ struct topo_tasks {
 		};
 	}
 
-	static void
-	axpby(flecsi::exec::accelerator s, topo_acc m, acc<rw> y, acc<ro> x, scalar alpha, scalar beta) noexcept {
+	template<bool inv, class F, class T>
+	static void axpy_self_future(flecsi::exec::accelerator s,
+	                             topo_acc m,
+	                             acc<rw> z,
+	                             acc<ro> x,
+	                             future<T> alpha_fut,
+	                             F f) noexcept {
+		scalar alpha = f(alpha_fut.get());
+		axpy_self<inv>(s, m, z, x, alpha);
+	}
+
+	static void axpby(flecsi::exec::accelerator s,
+	                  topo_acc m,
+	                  acc<rw> y,
+	                  acc<ro> x,
+	                  scalar alpha,
+	                  scalar beta) noexcept {
 		s.executor().forall(dof, util::dofs(m)) {
 			y[dof] = alpha * x[dof] + beta * y[dof];
 		};
 	}
 
-	static void abs_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x) noexcept {
+	static void
+	abs_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = std::abs(x[dof]); };
 	}
 
-	static void abs(flecsi::exec::accelerator s, topo_acc m, acc<wo> y, acc<ro> x) noexcept {
+	static void abs(flecsi::exec::accelerator s,
+	                topo_acc m,
+	                acc<wo> y,
+	                acc<ro> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { y[dof] = std::abs(x[dof]); };
 	}
 
-	static void add_scalar_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x, scalar alpha) noexcept {
+	static void add_scalar_self(flecsi::exec::accelerator s,
+	                            topo_acc m,
+	                            acc<rw> x,
+	                            scalar alpha) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] += alpha; };
 	}
 
-	static void add_scalar(flecsi::exec::accelerator s, topo_acc m, acc<wo> y, acc<ro> x, scalar alpha) noexcept {
+	static void add_scalar(flecsi::exec::accelerator s,
+	                       topo_acc m,
+	                       acc<wo> y,
+	                       acc<ro> x,
+	                       scalar alpha) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { y[dof] = x[dof] + alpha; };
 	}
 
-	static real lp_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> u, int p) noexcept {
-		real ret = s.executor().reduceall(dof,
-		                                  up,
-		                                  util::dofs(m),
-		                                  flecsi::exec::fold::sum,
-		                                  real) {
+	static real lp_norm_local(flecsi::exec::accelerator s,
+	                          topo_acc m,
+	                          acc<ro> u,
+	                          int p) noexcept {
+		real ret = s.executor().reduceall(
+			dof, up, util::dofs(m), flecsi::exec::fold::sum, real) {
 			up(std::pow(u[dof], p));
 		};
 
 		return ret;
 	}
 
-	static real l1_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
+	static real
+	l1_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
 		real ret = s.executor().reduceall(
 			dof, up, util::dofs(m), flecsi::exec::fold::sum, real) {
 			up(std::abs(u[dof]));
@@ -234,7 +327,8 @@ struct topo_tasks {
 		return ret;
 	}
 
-	static real l2_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept  {
+	static real
+	l2_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
 		auto ret = scalar_prod(s, m, u, u);
 		if constexpr (is_complex)
 			return ret.real();
@@ -242,53 +336,45 @@ struct topo_tasks {
 			return ret;
 	}
 
-	static real local_max(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
+	static real
+	local_max(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
 		if constexpr (is_complex) {
-			auto ret = s.executor().reduceall(dof,
-			                                  up,
-			                                  util::dofs(m),
-			                                  flecsi::exec::fold::max,
-			                                  real) {
+			auto ret = s.executor().reduceall(
+				dof, up, util::dofs(m), flecsi::exec::fold::max, real) {
 				up(u[dof].real());
 			};
 			return ret;
 		}
 		else {
-			auto ret = s.executor().reduceall(dof,
-			                                  up,
-			                                  util::dofs(m),
-			                                  flecsi::exec::fold::max,
-			                                  real) {
+			auto ret = s.executor().reduceall(
+				dof, up, util::dofs(m), flecsi::exec::fold::max, real) {
 				up(u[dof]);
 			};
 			return ret;
 		}
 	}
 
-	static real local_min(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
+	static real
+	local_min(flecsi::exec::accelerator s, topo_acc m, acc<ro> u) noexcept {
 		if constexpr (is_complex) {
-			auto ret = s.executor().reduceall(dof,
-			                                  up,
-			                                  util::dofs(m),
-			                                  flecsi::exec::fold::min,
-			                                  real) {
+			auto ret = s.executor().reduceall(
+				dof, up, util::dofs(m), flecsi::exec::fold::min, real) {
 				up(u[dof].real());
 			};
 			return ret;
 		}
 		else {
-			auto ret = s.executor().reduceall(dof,
-			                                  up,
-			                                  util::dofs(m),
-			                                  flecsi::exec::fold::min,
-			                                  real) {
+			auto ret = s.executor().reduceall(
+				dof, up, util::dofs(m), flecsi::exec::fold::min, real) {
 				up(u[dof]);
 			};
 			return ret;
 		}
 	}
 
-	static real inf_norm_local(flecsi::exec::accelerator s, topo_acc m, acc<ro> x) noexcept {
+	static real inf_norm_local(flecsi::exec::accelerator s,
+	                           topo_acc m,
+	                           acc<ro> x) noexcept {
 		auto ret = s.executor().reduceall(
 			dof, up, util::dofs(m), flecsi::exec::fold::max, real) {
 			up(std::abs(x[dof]));
@@ -313,7 +399,10 @@ struct topo_tasks {
 		};
 	}
 
-	static void dump(flecsi::exec::accelerator s, std::string_view pre, topo_acc m, acc<ro> x) {
+	static void dump(flecsi::exec::accelerator s,
+	                 std::string_view pre,
+	                 topo_acc m,
+	                 acc<ro> x) {
 		std::string fname{pre};
 		fname += "-" + std::to_string(s.launch().index);
 		std::ofstream ofile(fname);
diff --git a/flecsolve/vectors/operations/topo_view.hh b/flecsolve/vectors/operations/topo_view.hh
index b3946a5..ab5472d 100644
--- a/flecsolve/vectors/operations/topo_view.hh
+++ b/flecsolve/vectors/operations/topo_view.hh
@@ -44,8 +44,9 @@ struct topo_view {
 		static_assert(
 			std::is_same_v<typename Other::topo_t, typename vec_data::topo_t>);
 		static_assert(Other::space == vec_data::space);
-		scheduler().execute<
-			tasks::template copy<typename Other::template acc_all<flecsi::ro>>>(
+		scheduler()
+			.execute<tasks::template copy<
+				typename Other::template acc_all<flecsi::ro>>>(
 				flecsi::exec::on, x.topo(), z.ref(), x.ref());
 	}
 
@@ -160,8 +161,13 @@ struct topo_view {
 				flecsi::exec::on, z.topo(), x.ref(), z.ref(), alpha, beta);
 		}
 		else {
-			scheduler().execute<tasks::linear_sum>(
-				flecsi::exec::on, z.topo(), z.ref(), alpha, x.ref(), beta, y.ref());
+			scheduler().execute<tasks::linear_sum>(flecsi::exec::on,
+			                                       z.topo(),
+			                                       z.ref(),
+			                                       alpha,
+			                                       x.ref(),
+			                                       beta,
+			                                       y.ref());
 		}
 	}
 
@@ -169,11 +175,11 @@ struct topo_view {
 	axpy(scalar alpha, const vec_data & x, const vec_data & y, vec_data & z) {
 		if (z.fid() == x.fid()) {
 			scheduler().execute<tasks::template axpy_self<true>>(
-				                flecsi::exec::on, z.topo(), z.ref(), y.ref(), alpha);
+				flecsi::exec::on, z.topo(), z.ref(), y.ref(), alpha);
 		}
 		else if (z.fid() == y.fid()) {
 			scheduler().execute<tasks::template axpy_self<false>>(
-				                flecsi::exec::on, z.topo(), z.ref(), x.ref(), alpha);
+				flecsi::exec::on, z.topo(), z.ref(), x.ref(), alpha);
 		}
 		else {
 			scheduler().execute<tasks::axpy>(
@@ -181,6 +187,43 @@ struct topo_view {
 		}
 	}
 
+	template<class F, class T>
+	static void axpy(future_transform<future<T>, F> alpha,
+	                 const vec_data & x,
+	                 const vec_data & y,
+	                 vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "axpy: future type must be convertible to scalar");
+		if (z.fid() == x.fid()) {
+			scheduler().execute<tasks::template axpy_self_future<true, F, T>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				y.ref(),
+				alpha.fut,
+				alpha.f);
+		}
+		else if (z.fid() == y.fid()) {
+			scheduler().execute<tasks::template axpy_self_future<false, F, T>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				x.ref(),
+				alpha.fut,
+				alpha.f);
+		}
+		else {
+			scheduler().execute<tasks::template axpy_future<F, T>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				alpha.fut,
+				alpha.f,
+				x.ref(),
+				y.ref());
+		}
+	}
+
 	static void
 	axpby(scalar alpha, scalar beta, const vec_data & x, vec_data & z) {
 		scheduler().execute<tasks::axpby>(
@@ -210,33 +253,32 @@ struct topo_view {
 	}
 
 	static auto min(const vec_data & x) {
-		return scheduler().reduce<tasks::local_min,
-		                      flecsi::exec::fold::min>(
-			                      flecsi::exec::on, x.topo(), x.ref());
+		return scheduler().reduce<tasks::local_min, flecsi::exec::fold::min>(
+			flecsi::exec::on, x.topo(), x.ref());
 	}
 
 	static auto max(const vec_data & y) {
-		return scheduler().reduce<tasks::local_max,
-		                      flecsi::exec::fold::max>(
-			                      flecsi::exec::on, y.topo(), y.ref());
+		return scheduler().reduce<tasks::local_max, flecsi::exec::fold::max>(
+			flecsi::exec::on, y.topo(), y.ref());
 	}
 
 	template<unsigned short p>
 	static auto lp_norm_local(const vec_data & x) {
 		if constexpr (p == 1) {
-			return scheduler().reduce<tasks::l1_norm_local,
-			                      flecsi::exec::fold::sum>(
-				                      flecsi::exec::on, x.topo(), x.ref());
+			return scheduler()
+			    .reduce<tasks::l1_norm_local, flecsi::exec::fold::sum>(
+					flecsi::exec::on, x.topo(), x.ref());
 		}
 		else if constexpr (p == 2) {
-			return scheduler().reduce<tasks::l2_norm_local,
-			                      flecsi::exec::fold::sum>(flecsi::exec::on,
-			                                               x.topo(), x.ref());
+			return scheduler()
+			    .reduce<tasks::l2_norm_local, flecsi::exec::fold::sum>(
+					flecsi::exec::on, x.topo(), x.ref());
 		}
 		else {
-			return scheduler().reduce<tasks::lp_norm_local,
-			                      flecsi::exec::fold::sum>(
-				                      flecsi::exec::on, x.topo(), x.ref(), p).get();
+			return scheduler()
+			    .reduce<tasks::lp_norm_local, flecsi::exec::fold::sum>(
+					flecsi::exec::on, x.topo(), x.ref(), p)
+			    .get();
 		}
 	}
 
@@ -257,31 +299,32 @@ struct topo_view {
 	}
 
 	static auto inf_norm(const vec_data & x) {
-		return scheduler().reduce<tasks::inf_norm_local,
-		                      flecsi::exec::fold::max>(
-			                      flecsi::exec::on, x.topo(), x.ref());
+		return scheduler()
+		    .reduce<tasks::inf_norm_local, flecsi::exec::fold::max>(
+				flecsi::exec::on, x.topo(), x.ref());
 	}
 
 	static auto dot(const vec_data & x, const vec_data & y) {
-		return scheduler().reduce<tasks::scalar_prod,
-		                      flecsi::exec::fold::sum>(
-			                      flecsi::exec::on, x.topo(), x.ref(), y.ref());
+		return scheduler().reduce<tasks::scalar_prod, flecsi::exec::fold::sum>(
+			flecsi::exec::on, x.topo(), x.ref(), y.ref());
 	}
 
 	static auto global_size(const vec_data & x) {
-		return scheduler().reduce<tasks::local_size,
-		                      flecsi::exec::fold::sum>(x.topo());
+		return scheduler().reduce<tasks::local_size, flecsi::exec::fold::sum>(
+			x.topo());
 	}
 
 	static len_t local_size(const vec_data & x) {
 		len_t length;
-		scheduler().execute<tasks::get_local_size, flecsi::mpi>(x.topo(), &length);
+		scheduler().execute<tasks::get_local_size, flecsi::mpi>(x.topo(),
+		                                                        &length);
 		return length;
 	}
 
 	static void dump(std::string_view pre, const vec_data & x) {
 		// TODO: update for multiaccessor
-		scheduler().execute<tasks::dump>(flecsi::exec::on, pre, x.topo(), x.ref());
+		scheduler().execute<tasks::dump>(
+			flecsi::exec::on, pre, x.topo(), x.ref());
 	}
 
 	template<class F, class... Vecs>
diff --git a/flecsolve/vectors/test/flecsi_multivector.cc b/flecsolve/vectors/test/flecsi_multivector.cc
index d5a0a10..db19ad9 100644
--- a/flecsolve/vectors/test/flecsi_multivector.cc
+++ b/flecsolve/vectors/test/flecsi_multivector.cc
@@ -173,6 +173,9 @@ int vectest(flecsi::scheduler & s) {
 
 		tmp.axpy(7, x, y);
 		EXPECT_TRUE(run(tmp, msh, axpy));
+		tmp.zero();
+		tmp.axpy(flecsi::make_future(7), x, y);
+		EXPECT_TRUE(run(tmp, msh, axpy));
 
 		tmp.copy(y);
 		tmp.axpby(4, 11, z);
diff --git a/flecsolve/vectors/test/flecsi_vector.cc b/flecsolve/vectors/test/flecsi_vector.cc
index 9598bff..b80724d 100644
--- a/flecsolve/vectors/test/flecsi_vector.cc
+++ b/flecsolve/vectors/test/flecsi_vector.cc
@@ -60,17 +60,19 @@ using namespace std::complex_literals;
 template<class P>
 struct expected {
 	template<bool is_complex>
-	using scalar_type = std::conditional_t<is_complex, std::complex<double>, double>;
+	using scalar_type =
+		std::conditional_t<is_complex, std::complex<double>, double>;
 	static constexpr double ftol = 1e-8;
 
 	static int real(flecsi::exec::accelerator s,
-					testmesh::accessor<ro, ro> m,
-					realf::accessor<ro, na> x) noexcept {
-		UNIT(P::name) {
-			auto res = s.executor().reduceall(dof, up,
-											  m.dofs<testmesh::cells>(),
-											  flecsi::exec::fold::sum,
-											  double) {
+	                testmesh::accessor<ro, ro> m,
+	                realf::accessor<ro, na> x) noexcept {
+		UNIT (P::name) {
+			auto res = s.executor().reduceall(dof,
+			                                  up,
+			                                  m.dofs<testmesh::cells>(),
+			                                  flecsi::exec::fold::sum,
+			                                  double) {
 				auto gid = m.global_id(dof);
 				up(std::abs(P::real_answer(gid) - x[dof]));
 			};
@@ -79,9 +81,9 @@ struct expected {
 	}
 
 	static int complex(flecsi::exec::cpu s,
-					   testmesh::accessor<ro, ro> m,
-					   compf::accessor<ro, na> x) noexcept {
-		UNIT(P::name) {
+	                   testmesh::accessor<ro, ro> m,
+	                   compf::accessor<ro, na> x) noexcept {
+		UNIT (P::name) {
 			double diff = 0;
 			for (auto dof : m.dofs<testmesh::cells>()) {
 				auto gid = m.global_id(dof);
@@ -130,8 +132,10 @@ struct scalar_add_check : expected<scalar_add_check> {
 	static constexpr const char * name = "add scalar";
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto scalar_value() {
-		if constexpr (is_complex) return 1. + 1i;
-		else return 1;
+		if constexpr (is_complex)
+			return 1. + 1i;
+		else
+			return 1;
 	}
 
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
@@ -145,25 +149,30 @@ struct scalar_add_check : expected<scalar_add_check> {
 struct div_check : expected<div_check> {
 	static constexpr const char * name = "divide";
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
-		return rconv<1>(gid) / (rconv<0>(gid) + scalar_add_check::scalar_value<false>());
+		return rconv<1>(gid) /
+		       (rconv<0>(gid) + scalar_add_check::scalar_value<false>());
 	}
 	static auto complex_answer(double gid) {
 		return cconv<1>(gid) /
-			(cconv<0>(gid) + conv(scalar_add_check::scalar_value<true>()));
+		       (cconv<0>(gid) + conv(scalar_add_check::scalar_value<true>()));
 	}
 };
 
 struct scale_check : expected<scale_check> {
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto shift_value() {
-		if constexpr (!is_complex) return -1;
-		else return -1. - 1i;
+		if constexpr (!is_complex)
+			return -1;
+		else
+			return -1. - 1i;
 	}
 
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static double scale_value() {
-		if constexpr (!is_complex) return 2;
-		else return 2.4;
+		if constexpr (!is_complex)
+			return 2;
+		else
+			return 2.4;
 	}
 
 	static constexpr const char * name = "scale";
@@ -180,8 +189,10 @@ struct recip_check : expected<recip_check> {
 
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto shift_value() {
-		if constexpr (!is_complex) return 1;
-		else return 1. + 1i;
+		if constexpr (!is_complex)
+			return 1;
+		else
+			return 1. + 1i;
 	}
 
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
@@ -196,8 +207,10 @@ struct linsum_check : expected<linsum_check> {
 	static constexpr const char * name = "linear sum";
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto shift_value() {
-		if constexpr (!is_complex) return -1;
-		else return -1. - 1i;
+		if constexpr (!is_complex)
+			return -1;
+		else
+			return -1. - 1i;
 	}
 	inline static constexpr double alpha = 8;
 	inline static constexpr double beta = 9;
@@ -213,8 +226,10 @@ struct axpy_check : expected<axpy_check> {
 	static constexpr const char * name = "axpy";
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto alpha() {
-		if constexpr (!is_complex) return 7;
-		else return 7. + 3i;
+		if constexpr (!is_complex)
+			return 7;
+		else
+			return 7. + 3i;
 	}
 
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
@@ -229,14 +244,18 @@ struct axpby_check : expected<axpby_check> {
 	static constexpr const char * name = "axpby";
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto alpha() {
-		if constexpr (!is_complex) return 4;
-		else return 4.3 + 7i;
+		if constexpr (!is_complex)
+			return 4;
+		else
+			return 4.3 + 7i;
 	}
 
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto beta() {
-		if constexpr (!is_complex) return 11;
-		else return 11.8 + 3i;
+		if constexpr (!is_complex)
+			return 11;
+		else
+			return 11.8 + 3i;
 	}
 
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
@@ -245,7 +264,7 @@ struct axpby_check : expected<axpby_check> {
 
 	static auto complex_answer(double gid) {
 		return cconv<2>(gid) * conv(alpha<true>()) +
-			cconv<1>(gid) * conv(beta<true>());
+		       cconv<1>(gid) * conv(beta<true>());
 	}
 };
 
@@ -253,8 +272,10 @@ struct abs_check : expected<abs_check> {
 	static constexpr const char * name = "abs";
 	template<bool is_complex>
 	FLECSI_INLINE_TARGET static auto shift() {
-		if constexpr (!is_complex) return -4;
-		else return -4. -4i;
+		if constexpr (!is_complex)
+			return -4;
+		else
+			return -4. - 4i;
 	}
 
 	FLECSI_INLINE_TARGET static double real_answer(double gid) {
@@ -270,21 +291,26 @@ struct dot_check : expected<dot_check> {
 	template<bool is_complex>
 	static auto answer() {
 		auto [alpha, beta] = scalars<is_complex>();
-		return scalar_type<is_complex>(32) * alpha * (is_complex ? std::conj(beta) : beta);
+		return scalar_type<is_complex>(32) * alpha *
+		       (is_complex ? std::conj(beta) : beta);
 	}
 
 	template<bool is_complex>
 	static auto scalars() {
-		if constexpr (is_complex) return std::pair{1.2 + 3i, 7. + 2.3i};
-		else return std::pair{1.5, 3.8};
+		if constexpr (is_complex)
+			return std::pair{1.2 + 3i, 7. + 2.3i};
+		else
+			return std::pair{1.5, 3.8};
 	}
 };
 
 struct l1norm_check : expected<l1norm_check> {
 	template<bool is_complex>
 	static auto scalar() {
-		if constexpr (is_complex) return 1.3 + 8.7i;
-		else return -3.141719;
+		if constexpr (is_complex)
+			return 1.3 + 8.7i;
+		else
+			return -3.141719;
 	}
 	template<bool is_complex>
 	static double answer() {
@@ -300,9 +326,9 @@ struct l2norm_check : expected<l2norm_check> {
 
 	template<bool is_complex>
 	static auto answer() {
-		return std::sqrt(scalar_type<is_complex>(32) *
-						 scalar<is_complex>() *
-						 (is_complex ? std::conj(scalar<is_complex>()) : scalar<is_complex>()));
+		return std::sqrt(scalar_type<is_complex>(32) * scalar<is_complex>() *
+		                 (is_complex ? std::conj(scalar<is_complex>())
+		                             : scalar<is_complex>()));
 	}
 };
 
@@ -321,14 +347,17 @@ int vectest(flecsi::scheduler & s) {
 		auto [x, y, z, tmp] = create(xd, yd, zd, tmpd);
 		auto [x_c, y_c, z_c, tmp_c] = create(xd_c, yd_c, zd_c, tmpd_c);
 
-		auto test_vecops = [&](auto & x, auto & y, auto & z, auto & tmp) -> int {
-			constexpr bool is_complex = num_traits<typename std::decay_t<decltype(x)>::scalar>::is_complex;
+		auto test_vecops =
+			[&](auto & x, auto & y, auto & z, auto & tmp) -> int {
+			constexpr bool is_complex = num_traits<
+				typename std::decay_t<decltype(x)>::scalar>::is_complex;
 			auto check_answer = [&](auto expect, auto & vec) {
 				if constexpr (is_complex) {
 					EXPECT_EQ(s.test<decltype(expect)::complex>(
 								  flecsi::exec::on, msh, vec.data.ref()),
 					          0);
-				} else {
+				}
+				else {
 					EXPECT_EQ(s.test<decltype(expect)::real>(
 								  flecsi::exec::on, msh, vec.data.ref()),
 					          0);
@@ -363,9 +392,15 @@ int vectest(flecsi::scheduler & s) {
 
 			tmp.axpy(axpy_check::alpha<is_complex>(), x, y);
 			check_answer(axpy_check{}, tmp);
+			tmp.zero();
+			tmp.axpy(
+				flecsi::make_future(axpy_check::alpha<is_complex>()), x, y);
+			check_answer(axpy_check{}, tmp);
 
 			tmp.copy(y);
-			tmp.axpby(axpby_check::alpha<is_complex>(), axpby_check::beta<is_complex>(), z);
+			tmp.axpby(axpby_check::alpha<is_complex>(),
+			          axpby_check::beta<is_complex>(),
+			          z);
 			check_answer(axpby_check{}, tmp);
 
 			tmp.add_scalar(y, abs_check::shift<is_complex>());
@@ -383,11 +418,14 @@ int vectest(flecsi::scheduler & s) {
 				auto [alpha, beta] = dot_check::scalars<is_complex>();
 				a.set_scalar(alpha);
 				b.set_scalar(beta);
-				EXPECT_LT(dot_check::abs_error<is_complex>(a.dot(b).get()), dot_check::ftol);
+				EXPECT_LT(dot_check::abs_error<is_complex>(a.dot(b).get()),
+				          dot_check::ftol);
 
 				a.set_scalar(l1norm_check::scalar<is_complex>());
-				EXPECT_LT(l1norm_check::abs_error<is_complex>(a.l1norm().get()), l1norm_check::ftol);
-				EXPECT_LT(l2norm_check::abs_error<is_complex>(a.l2norm().get()), l2norm_check::ftol);
+				EXPECT_LT(l1norm_check::abs_error<is_complex>(a.l1norm().get()),
+				          l1norm_check::ftol);
+				EXPECT_LT(l2norm_check::abs_error<is_complex>(a.l2norm().get()),
+				          l2norm_check::ftol);
 			};
 			test_reductions(x, y);
 

From 2bfb5477c5642f8ebee99d8eb5074999f3aad343 Mon Sep 17 00:00:00 2001
From: Andrew Reisner <areisner@lanl.gov>
Date: Tue, 11 Aug 2026 12:51:10 -0600
Subject: [PATCH 2/4] Add future support for scalar parameters in all vec ops.

---
 flecsolve/vectors/core.hh                    |  85 ++++--
 flecsolve/vectors/operations/multi.hh        |  92 ++++---
 flecsolve/vectors/operations/topo_tasks.hh   | 156 +++++++++++
 flecsolve/vectors/operations/topo_view.hh    | 257 +++++++++++++++++++
 flecsolve/vectors/seq.hh                     |  69 +++--
 flecsolve/vectors/test/flecsi_multivector.cc |  11 +-
 flecsolve/vectors/test/flecsi_vector.cc      |  24 +-
 7 files changed, 602 insertions(+), 92 deletions(-)

diff --git a/flecsolve/vectors/core.hh b/flecsolve/vectors/core.hh
index 6fb3448..f9b928c 100644
--- a/flecsolve/vectors/core.hh
+++ b/flecsolve/vectors/core.hh
@@ -75,7 +75,12 @@ struct core : Config {
 	 * \f$\mathit{this}_i = val\f$
 	 * \param[in] val scalar
 	 */
-	void set_scalar(scalar val) { ops::set_to_scalar(val, data); }
+	template<class T,
+	         std::enable_if_t<is_scalar_or_future_v<std::decay_t<T>, scalar>,
+	                          bool> = true>
+	void set_scalar(T val) {
+		ops::set_to_scalar(defer(val), data);
+	}
 
 	/**
 	 * Scale vector components.
@@ -84,9 +89,14 @@ struct core : Config {
 	 * \param[in] alpha scalar
 	 * \param[in] x vector
 	 */
-	template<class V, std::enable_if_t<is_vector_v<V>, bool> = true>
-	void scale(scalar alpha, const V & x) {
-		ops::scale(alpha, x.data, data);
+	template<
+		class Alpha,
+		class V,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_vector_v<V>, bool> = true>
+	void scale(Alpha alpha, const V & x) {
+		ops::scale(defer(alpha), x.data, data);
 	}
 
 	/**
@@ -95,7 +105,13 @@ struct core : Config {
 	 * \f$\mathit{this}_i = alpha * \mathit{this}_i\f$
 	 * \param[in] alpha scalar
 	 */
-	void scale(scalar alpha) { ops::scale(alpha, data); }
+	template<
+		class Alpha,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	void scale(Alpha alpha) {
+		ops::scale(defer(alpha), data);
+	}
 
 	/**
 	 * Component-wise addition of two vectors.
@@ -172,12 +188,19 @@ struct core : Config {
 	 *
 	 * \f$\mathit{this}_i = alpha * x_i + beta * y_i\f$
 	 */
-	template<class V1,
-	         class V2,
-	         std::enable_if_t<is_vector_v<V1>, bool> = true,
-	         std::enable_if_t<is_vector_v<V2>, bool> = true>
-	void linear_sum(scalar alpha, const V1 & x, scalar beta, const V2 & y) {
-		ops::linear_sum(alpha, x.data, beta, y.data, data);
+	template<
+		class Alpha,
+		class Beta,
+		class V1,
+		class V2,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Beta>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_vector_v<V1>, bool> = true,
+		std::enable_if_t<is_vector_v<V2>, bool> = true>
+	void linear_sum(Alpha alpha, const V1 & x, Beta beta, const V2 & y) {
+		ops::linear_sum(defer(alpha), x.data, defer(beta), y.data, data);
 	}
 
 	/**
@@ -185,14 +208,15 @@ struct core : Config {
 	 *
 	 * \f$\mathit{this}_i = alpha x_i + y_i\f$
 	 */
-	template<class Alpha,
-	         class V1,
-	         class V2,
-	         std::enable_if_t<is_vector_v<V1>, bool> = true,
-	         std::enable_if_t<is_vector_v<V2>, bool> = true>
+	template<
+		class Alpha,
+		class V1,
+		class V2,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_vector_v<V1>, bool> = true,
+		std::enable_if_t<is_vector_v<V2>, bool> = true>
 	void axpy(Alpha alpha, const V1 & x, const V2 & y) {
-		static_assert(is_scalar_or_future_v<Alpha, scalar>,
-		              "axpy: alpha must be scalar or future");
 		ops::axpy(defer(alpha), x.data, y.data, data);
 	}
 
@@ -201,9 +225,17 @@ struct core : Config {
 	 *
 	 * \f$\mathit{this}_i = alpha * x_i + beta * \mathit{this}_i\f$
 	 */
-	template<class V, std::enable_if_t<is_vector_v<V>, bool> = true>
-	void axpby(scalar alpha, scalar beta, const V & x) {
-		ops::axpby(alpha, beta, x.data, data);
+	template<
+		class Alpha,
+		class Beta,
+		class V,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Beta>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_vector_v<V>, bool> = true>
+	void axpby(Alpha alpha, Beta beta, const V & x) {
+		ops::axpby(defer(alpha), defer(beta), x.data, data);
 	}
 
 	/**
@@ -221,9 +253,14 @@ struct core : Config {
 	 *
 	 * \f$\mathit{this}_i = alpha x_i\f$
 	 */
-	template<class V, std::enable_if_t<is_vector_v<V>, bool> = true>
-	void add_scalar(const V & x, scalar alpha) {
-		ops::add_scalar(x.data, alpha, data);
+	template<
+		class V,
+		class Alpha,
+		std::enable_if_t<is_vector_v<V>, bool> = true,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	void add_scalar(const V & x, Alpha alpha) {
+		ops::add_scalar(x.data, defer(alpha), data);
 	}
 
 	/**
diff --git a/flecsolve/vectors/operations/multi.hh b/flecsolve/vectors/operations/multi.hh
index 409dc9a..60e54c6 100644
--- a/flecsolve/vectors/operations/multi.hh
+++ b/flecsolve/vectors/operations/multi.hh
@@ -44,19 +44,31 @@ struct multi {
 		apply([](auto & x) { x.zero(); }, make_is(), x.components);
 	}
 
-	static void set_to_scalar(scalar alpha, vec_data & x) {
-		apply([alpha](auto & x) { x.set_scalar(alpha); },
+	template<
+		class Alpha,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	static void set_to_scalar(Alpha alpha, vec_data & x) {
+		apply([&alpha](auto & x) { x.set_scalar(alpha); },
 		      make_is(),
 		      x.components);
 	}
 
-	static void scale(scalar alpha, vec_data & x) {
-		apply([alpha](auto & x) { x.scale(alpha); }, make_is(), x.components);
+	template<
+		class Alpha,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	static void scale(Alpha alpha, vec_data & x) {
+		apply([&alpha](auto & x) { x.scale(alpha); }, make_is(), x.components);
 	}
 
-	template<class T>
-	static void scale(scalar alpha, const T & x, vec_data & y) {
-		apply([alpha](auto & v, const auto & y) { v.scale(alpha, y); },
+	template<
+		class Alpha,
+		class T,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	static void scale(Alpha alpha, const T & x, vec_data & y) {
+		apply([&alpha](auto & v, const auto & y) { v.scale(alpha, y); },
 		      make_is(),
 		      y.components,
 		      x.components);
@@ -108,14 +120,22 @@ struct multi {
 		      x.components);
 	}
 
-	template<class T0, class T1>
-	static void linear_sum(scalar alpha,
+	template<
+		class Alpha,
+		class Beta,
+		class T0,
+		class T1,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Beta>, scalar>,
+	                     bool> = true>
+	static void linear_sum(Alpha alpha,
 	                       const T0 & x,
-	                       scalar beta,
+	                       Beta beta,
 	                       const T1 & y,
 	                       vec_data & z) {
 		apply(
-			[alpha, beta](auto & z, const auto & x, const auto & y) {
+			[&alpha, &beta](auto & z, const auto & x, const auto & y) {
 				z.linear_sum(alpha, x, beta, y);
 			},
 			make_is(),
@@ -124,24 +144,13 @@ struct multi {
 			y.components);
 	}
 
-	template<class T0, class T1>
-	static void axpy(scalar alpha, const T0 & x, const T1 & y, vec_data & z) {
-		apply([alpha](auto & z,
-		              const auto & x,
-		              const auto & y) { z.axpy(alpha, x, y); },
-		      make_is(),
-		      z.components,
-		      x.components,
-		      y.components);
-	}
-
-	template<class F, class T, class T0, class T1>
-	static void axpy(future_transform<future<T>, F> alpha,
-	                 const T0 & x,
-	                 const T1 & y,
-	                 vec_data & z) {
-		static_assert(std::is_convertible_v<T, scalar>,
-		              "axpy: future type must be convertible to scalar");
+	template<
+		class Alpha,
+		class T0,
+		class T1,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	static void axpy(Alpha alpha, const T0 & x, const T1 & y, vec_data & z) {
 		apply([&alpha](auto & z,
 		               const auto & x,
 		               const auto & y) { z.axpy(alpha, x, y); },
@@ -151,10 +160,17 @@ struct multi {
 		      y.components);
 	}
 
-	template<class T>
-	static void axpby(scalar alpha, scalar beta, const T & x, vec_data & z) {
-		apply([alpha, beta](auto & z,
-		                    const auto & x) { z.axpby(alpha, beta, x); },
+	template<
+		class Alpha,
+		class Beta,
+		class T,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Beta>, scalar>,
+	                     bool> = true>
+	static void axpby(Alpha alpha, Beta beta, const T & x, vec_data & z) {
+		apply([&alpha, &beta](auto & z,
+		                      const auto & x) { z.axpby(alpha, beta, x); },
 		      make_is(),
 		      z.components,
 		      x.components);
@@ -168,9 +184,13 @@ struct multi {
 		      x.components);
 	}
 
-	template<class T>
-	static void add_scalar(const T & x, scalar alpha, vec_data & y) {
-		apply([alpha](auto & z, const auto & x) { z.add_scalar(x, alpha); },
+	template<
+		class T,
+		class Alpha,
+		std::enable_if_t<is_scalar_or_future_v<std::decay_t<Alpha>, scalar>,
+	                     bool> = true>
+	static void add_scalar(const T & x, Alpha alpha, vec_data & y) {
+		apply([&alpha](auto & z, const auto & x) { z.add_scalar(x, alpha); },
 		      make_is(),
 		      y.components,
 		      x.components);
diff --git a/flecsolve/vectors/operations/topo_tasks.hh b/flecsolve/vectors/operations/topo_tasks.hh
index 17085a4..a01f68a 100644
--- a/flecsolve/vectors/operations/topo_tasks.hh
+++ b/flecsolve/vectors/operations/topo_tasks.hh
@@ -23,6 +23,7 @@ to do so.
 #include <flecsi/execution.hh>
 #include <flecsi/data.hh>
 
+#include "flecsolve/util/future.hh"
 #include "flecsolve/util/traits.hh"
 #include "flecsolve/vectors/data/topo_view.hh"
 
@@ -74,6 +75,16 @@ struct topo_tasks {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = val; };
 	}
 
+	template<class F, class T>
+	static void set_to_scalar_future(flecsi::exec::accelerator s,
+	                                 topo_acc m,
+	                                 acc<wo> x,
+	                                 future<T> fut,
+	                                 F f) noexcept {
+		scalar val = f(fut.get());
+		set_to_scalar(s, m, x, val);
+	}
+
 	static void scale_self(flecsi::exec::accelerator s,
 	                       topo_acc m,
 	                       acc<rw> x,
@@ -81,6 +92,15 @@ struct topo_tasks {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] *= val; };
 	}
 
+	template<class F, class T>
+	static void scale_self_future(flecsi::exec::accelerator s,
+	                              topo_acc m,
+	                              acc<rw> x,
+	                              future<T> fut,
+	                              F f) noexcept {
+		scale_self(s, m, x, f(fut.get()));
+	}
+
 	static void scale(flecsi::exec::accelerator s,
 	                  topo_acc m,
 	                  acc<ro> x,
@@ -89,6 +109,16 @@ struct topo_tasks {
 		s.executor().forall(dof, util::dofs(m)) { y[dof] = x[dof] * val; };
 	}
 
+	template<class F, class T>
+	static void scale_future(flecsi::exec::accelerator s,
+	                         topo_acc m,
+	                         acc<ro> x,
+	                         acc<wo> y,
+	                         future<T> fut,
+	                         F f) noexcept {
+		scale(s, m, x, y, f(fut.get()));
+	}
+
 	template<class OtherAcc>
 	static void copy(flecsi::exec::accelerator s,
 	                 topo_acc,
@@ -200,6 +230,44 @@ struct topo_tasks {
 		};
 	}
 
+	template<class F, class T>
+	static void linear_sum_alpha_future(flecsi::exec::accelerator s,
+	                                    topo_acc m,
+	                                    acc<wo> z,
+	                                    future<T> alpha_fut,
+	                                    F alpha_f,
+	                                    acc<ro> x,
+	                                    scalar beta,
+	                                    acc<ro> y) noexcept {
+		linear_sum(s, m, z, alpha_f(alpha_fut.get()), x, beta, y);
+	}
+
+	template<class F, class T>
+	static void linear_sum_beta_future(flecsi::exec::accelerator s,
+	                                   topo_acc m,
+	                                   acc<wo> z,
+	                                   scalar alpha,
+	                                   acc<ro> x,
+	                                   future<T> beta_fut,
+	                                   F beta_f,
+	                                   acc<ro> y) noexcept {
+		linear_sum(s, m, z, alpha, x, beta_f(beta_fut.get()), y);
+	}
+
+	template<class AlphaF, class AlphaT, class BetaF, class BetaT>
+	static void linear_sum_future(flecsi::exec::accelerator s,
+	                              topo_acc m,
+	                              acc<wo> z,
+	                              future<AlphaT> alpha_fut,
+	                              AlphaF alpha_f,
+	                              acc<ro> x,
+	                              future<BetaT> beta_fut,
+	                              BetaF beta_f,
+	                              acc<ro> y) noexcept {
+		linear_sum(
+			s, m, z, alpha_f(alpha_fut.get()), x, beta_f(beta_fut.get()), y);
+	}
+
 	template<bool inv>
 	static void linear_sum_self(flecsi::exec::accelerator s,
 	                            topo_acc m,
@@ -217,6 +285,41 @@ struct topo_tasks {
 		};
 	}
 
+	template<bool inv, class F, class T>
+	static void linear_sum_self_alpha_future(flecsi::exec::accelerator s,
+	                                         topo_acc m,
+	                                         acc<rw> z,
+	                                         acc<ro> x,
+	                                         future<T> alpha_fut,
+	                                         F alpha_f,
+	                                         scalar beta) noexcept {
+		linear_sum_self<inv>(s, m, z, x, alpha_f(alpha_fut.get()), beta);
+	}
+
+	template<bool inv, class F, class T>
+	static void linear_sum_self_beta_future(flecsi::exec::accelerator s,
+	                                        topo_acc m,
+	                                        acc<rw> z,
+	                                        acc<ro> x,
+	                                        scalar alpha,
+	                                        future<T> beta_fut,
+	                                        F beta_f) noexcept {
+		linear_sum_self<inv>(s, m, z, x, alpha, beta_f(beta_fut.get()));
+	}
+
+	template<bool inv, class AlphaF, class AlphaT, class BetaF, class BetaT>
+	static void linear_sum_self_future(flecsi::exec::accelerator s,
+	                                   topo_acc m,
+	                                   acc<rw> z,
+	                                   acc<ro> x,
+	                                   future<AlphaT> alpha_fut,
+	                                   AlphaF alpha_f,
+	                                   future<BetaT> beta_fut,
+	                                   BetaF beta_f) noexcept {
+		linear_sum_self<inv>(
+			s, m, z, x, alpha_f(alpha_fut.get()), beta_f(beta_fut.get()));
+	}
+
 	static void axpy(flecsi::exec::accelerator s,
 	                 topo_acc m,
 	                 acc<wo> z,
@@ -278,6 +381,40 @@ struct topo_tasks {
 		};
 	}
 
+	template<class F, class T>
+	static void axpby_alpha_future(flecsi::exec::accelerator s,
+	                               topo_acc m,
+	                               acc<rw> y,
+	                               acc<ro> x,
+	                               future<T> alpha_fut,
+	                               F alpha_f,
+	                               scalar beta) noexcept {
+		axpby(s, m, y, x, alpha_f(alpha_fut.get()), beta);
+	}
+
+	template<class F, class T>
+	static void axpby_beta_future(flecsi::exec::accelerator s,
+	                              topo_acc m,
+	                              acc<rw> y,
+	                              acc<ro> x,
+	                              scalar alpha,
+	                              future<T> beta_fut,
+	                              F beta_f) noexcept {
+		axpby(s, m, y, x, alpha, beta_f(beta_fut.get()));
+	}
+
+	template<class AlphaF, class AlphaT, class BetaF, class BetaT>
+	static void axpby_future(flecsi::exec::accelerator s,
+	                         topo_acc m,
+	                         acc<rw> y,
+	                         acc<ro> x,
+	                         future<AlphaT> alpha_fut,
+	                         AlphaF alpha_f,
+	                         future<BetaT> beta_fut,
+	                         BetaF beta_f) noexcept {
+		axpby(s, m, y, x, alpha_f(alpha_fut.get()), beta_f(beta_fut.get()));
+	}
+
 	static void
 	abs_self(flecsi::exec::accelerator s, topo_acc m, acc<rw> x) noexcept {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] = std::abs(x[dof]); };
@@ -297,6 +434,15 @@ struct topo_tasks {
 		s.executor().forall(dof, util::dofs(m)) { x[dof] += alpha; };
 	}
 
+	template<class F, class T>
+	static void add_scalar_self_future(flecsi::exec::accelerator s,
+	                                   topo_acc m,
+	                                   acc<rw> x,
+	                                   future<T> alpha_fut,
+	                                   F alpha_f) noexcept {
+		add_scalar_self(s, m, x, alpha_f(alpha_fut.get()));
+	}
+
 	static void add_scalar(flecsi::exec::accelerator s,
 	                       topo_acc m,
 	                       acc<wo> y,
@@ -305,6 +451,16 @@ struct topo_tasks {
 		s.executor().forall(dof, util::dofs(m)) { y[dof] = x[dof] + alpha; };
 	}
 
+	template<class F, class T>
+	static void add_scalar_future(flecsi::exec::accelerator s,
+	                              topo_acc m,
+	                              acc<wo> y,
+	                              acc<ro> x,
+	                              future<T> alpha_fut,
+	                              F alpha_f) noexcept {
+		add_scalar(s, m, y, x, alpha_f(alpha_fut.get()));
+	}
+
 	static real lp_norm_local(flecsi::exec::accelerator s,
 	                          topo_acc m,
 	                          acc<ro> u,
diff --git a/flecsolve/vectors/operations/topo_view.hh b/flecsolve/vectors/operations/topo_view.hh
index ab5472d..498b166 100644
--- a/flecsolve/vectors/operations/topo_view.hh
+++ b/flecsolve/vectors/operations/topo_view.hh
@@ -64,11 +64,29 @@ struct topo_view {
 			flecsi::exec::on, x.topo(), x.ref(), alpha);
 	}
 
+	template<class F, class T>
+	static void set_to_scalar(future_transform<future<T>, F> alpha,
+	                          vec_data & x) {
+		static_assert(
+			std::is_convertible_v<T, scalar>,
+			"set_to_scalar: future type must be convertible to scalar");
+		scheduler().execute<tasks::template set_to_scalar_future<F, T>>(
+			flecsi::exec::on, x.topo(), x.ref(), alpha.fut, alpha.f);
+	}
+
 	static void scale(scalar alpha, vec_data & x) {
 		scheduler().execute<tasks::scale_self>(
 			flecsi::exec::on, x.topo(), x.ref(), alpha);
 	}
 
+	template<class F, class T>
+	static void scale(future_transform<future<T>, F> alpha, vec_data & x) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "scale: future type must be convertible to scalar");
+		scheduler().execute<tasks::template scale_self_future<F, T>>(
+			flecsi::exec::on, x.topo(), x.ref(), alpha.fut, alpha.f);
+	}
+
 	static void scale(scalar alpha, const vec_data & x, vec_data & y) {
 		flog_assert(x.fid() != y.fid(),
 		            "scale operation: vector data cannot be the same");
@@ -76,6 +94,18 @@ struct topo_view {
 			flecsi::exec::on, x.topo(), x.ref(), y.ref(), alpha);
 	}
 
+	template<class F, class T>
+	static void scale(future_transform<future<T>, F> alpha,
+	                  const vec_data & x,
+	                  vec_data & y) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "scale: future type must be convertible to scalar");
+		flog_assert(x.fid() != y.fid(),
+		            "scale operation: vector data cannot be the same");
+		scheduler().execute<tasks::template scale_future<F, T>>(
+			flecsi::exec::on, x.topo(), x.ref(), y.ref(), alpha.fut, alpha.f);
+	}
+
 	static void add(const vec_data & x, const vec_data & y, vec_data & z) {
 		if (x.fid() == z.fid()) {
 			scheduler().execute<tasks::add_self>(
@@ -171,6 +201,156 @@ struct topo_view {
 		}
 	}
 
+	template<class F, class T>
+	static void linear_sum(future_transform<future<T>, F> alpha,
+	                       const vec_data & x,
+	                       scalar beta,
+	                       const vec_data & y,
+	                       vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "linear_sum: future type must be convertible to scalar");
+		if (z.fid() == x.fid()) {
+			scheduler()
+				.execute<
+					tasks::template linear_sum_self_alpha_future<true, F, T>>(
+					flecsi::exec::on,
+					z.topo(),
+					z.ref(),
+					y.ref(),
+					alpha.fut,
+					alpha.f,
+					beta);
+		}
+		else if (z.fid() == y.fid()) {
+			scheduler()
+				.execute<
+					tasks::template linear_sum_self_alpha_future<false, F, T>>(
+					flecsi::exec::on,
+					z.topo(),
+					x.ref(),
+					z.ref(),
+					alpha.fut,
+					alpha.f,
+					beta);
+		}
+		else {
+			scheduler().execute<tasks::template linear_sum_alpha_future<F, T>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				alpha.fut,
+				alpha.f,
+				x.ref(),
+				beta,
+				y.ref());
+		}
+	}
+
+	template<class F, class T>
+	static void linear_sum(scalar alpha,
+	                       const vec_data & x,
+	                       future_transform<future<T>, F> beta,
+	                       const vec_data & y,
+	                       vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "linear_sum: future type must be convertible to scalar");
+		if (z.fid() == x.fid()) {
+			scheduler()
+				.execute<
+					tasks::template linear_sum_self_beta_future<true, F, T>>(
+					flecsi::exec::on,
+					z.topo(),
+					z.ref(),
+					y.ref(),
+					alpha,
+					beta.fut,
+					beta.f);
+		}
+		else if (z.fid() == y.fid()) {
+			scheduler()
+				.execute<
+					tasks::template linear_sum_self_beta_future<false, F, T>>(
+					flecsi::exec::on,
+					z.topo(),
+					x.ref(),
+					z.ref(),
+					alpha,
+					beta.fut,
+					beta.f);
+		}
+		else {
+			scheduler().execute<tasks::template linear_sum_beta_future<F, T>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				alpha,
+				x.ref(),
+				beta.fut,
+				beta.f,
+				y.ref());
+		}
+	}
+
+	template<class AlphaF, class AlphaT, class BetaF, class BetaT>
+	static void linear_sum(future_transform<future<AlphaT>, AlphaF> alpha,
+	                       const vec_data & x,
+	                       future_transform<future<BetaT>, BetaF> beta,
+	                       const vec_data & y,
+	                       vec_data & z) {
+		static_assert(std::is_convertible_v<AlphaT, scalar>,
+		              "linear_sum: future type must be convertible to scalar");
+		static_assert(std::is_convertible_v<BetaT, scalar>,
+		              "linear_sum: future type must be convertible to scalar");
+		if (z.fid() == x.fid()) {
+			scheduler()
+				.execute<tasks::template linear_sum_self_future<true,
+			                                                    AlphaF,
+			                                                    AlphaT,
+			                                                    BetaF,
+			                                                    BetaT>>(
+					flecsi::exec::on,
+					z.topo(),
+					z.ref(),
+					y.ref(),
+					alpha.fut,
+					alpha.f,
+					beta.fut,
+					beta.f);
+		}
+		else if (z.fid() == y.fid()) {
+			scheduler()
+				.execute<tasks::template linear_sum_self_future<false,
+			                                                    AlphaF,
+			                                                    AlphaT,
+			                                                    BetaF,
+			                                                    BetaT>>(
+					flecsi::exec::on,
+					z.topo(),
+					x.ref(),
+					z.ref(),
+					alpha.fut,
+					alpha.f,
+					beta.fut,
+					beta.f);
+		}
+		else {
+			scheduler()
+				.execute<tasks::template linear_sum_future<AlphaF,
+			                                               AlphaT,
+			                                               BetaF,
+			                                               BetaT>>(
+					flecsi::exec::on,
+					z.topo(),
+					z.ref(),
+					alpha.fut,
+					alpha.f,
+					x.ref(),
+					beta.fut,
+					beta.f,
+					y.ref());
+		}
+	}
+
 	static void
 	axpy(scalar alpha, const vec_data & x, const vec_data & y, vec_data & z) {
 		if (z.fid() == x.fid()) {
@@ -230,6 +410,62 @@ struct topo_view {
 			flecsi::exec::on, z.topo(), z.ref(), x.ref(), alpha, beta);
 	}
 
+	template<class F, class T>
+	static void axpby(future_transform<future<T>, F> alpha,
+	                  scalar beta,
+	                  const vec_data & x,
+	                  vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "axpby: future type must be convertible to scalar");
+		scheduler().execute<tasks::template axpby_alpha_future<F, T>>(
+			flecsi::exec::on,
+			z.topo(),
+			z.ref(),
+			x.ref(),
+			alpha.fut,
+			alpha.f,
+			beta);
+	}
+
+	template<class F, class T>
+	static void axpby(scalar alpha,
+	                  future_transform<future<T>, F> beta,
+	                  const vec_data & x,
+	                  vec_data & z) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "axpby: future type must be convertible to scalar");
+		scheduler().execute<tasks::template axpby_beta_future<F, T>>(
+			flecsi::exec::on,
+			z.topo(),
+			z.ref(),
+			x.ref(),
+			alpha,
+			beta.fut,
+			beta.f);
+	}
+
+	template<class AlphaF, class AlphaT, class BetaF, class BetaT>
+	static void axpby(future_transform<future<AlphaT>, AlphaF> alpha,
+	                  future_transform<future<BetaT>, BetaF> beta,
+	                  const vec_data & x,
+	                  vec_data & z) {
+		static_assert(std::is_convertible_v<AlphaT, scalar>,
+		              "axpby: future type must be convertible to scalar");
+		static_assert(std::is_convertible_v<BetaT, scalar>,
+		              "axpby: future type must be convertible to scalar");
+		scheduler()
+			.execute<
+				tasks::template axpby_future<AlphaF, AlphaT, BetaF, BetaT>>(
+				flecsi::exec::on,
+				z.topo(),
+				z.ref(),
+				x.ref(),
+				alpha.fut,
+				alpha.f,
+				beta.fut,
+				beta.f);
+	}
+
 	static void abs(const vec_data & x, vec_data & y) {
 		if (y.fid() == x.fid()) {
 			scheduler().execute<tasks::abs_self>(
@@ -252,6 +488,27 @@ struct topo_view {
 		}
 	}
 
+	template<class F, class T>
+	static void add_scalar(const vec_data & x,
+	                       future_transform<future<T>, F> alpha,
+	                       vec_data & y) {
+		static_assert(std::is_convertible_v<T, scalar>,
+		              "add_scalar: future type must be convertible to scalar");
+		if (x.fid() == y.fid()) {
+			scheduler().execute<tasks::template add_scalar_self_future<F, T>>(
+				flecsi::exec::on, y.topo(), y.ref(), alpha.fut, alpha.f);
+		}
+		else {
+			scheduler().execute<tasks::template add_scalar_future<F, T>>(
+				flecsi::exec::on,
+				y.topo(),
+				y.ref(),
+				x.ref(),
+				alpha.fut,
+				alpha.f);
+		}
+	}
+
 	static auto min(const vec_data & x) {
 		return scheduler().reduce<tasks::local_min, flecsi::exec::fold::min>(
 			flecsi::exec::on, x.topo(), x.ref());
diff --git a/flecsolve/vectors/seq.hh b/flecsolve/vectors/seq.hh
index 898faee..8dd3e2e 100644
--- a/flecsolve/vectors/seq.hh
+++ b/flecsolve/vectors/seq.hh
@@ -114,19 +114,22 @@ struct seq_ops {
 			x);
 	}
 
-	template<class D>
-	static void set_to_scalar(scalar alpha, D & x) {
-		fordofs([&](len_t i) { x[i] = alpha; }, x);
+	template<class Alpha, class D>
+	static void set_to_scalar(Alpha alpha, D & x) {
+		const scalar value = scalar_value(std::move(alpha));
+		fordofs([&](len_t i) { x[i] = value; }, x);
 	}
 
-	template<class D>
-	static void scale(scalar alpha, D & x) {
-		fordofs([&](len_t i) { x[i] *= alpha; }, x);
+	template<class Alpha, class D>
+	static void scale(Alpha alpha, D & x) {
+		const scalar value = scalar_value(std::move(alpha));
+		fordofs([&](len_t i) { x[i] *= value; }, x);
 	}
 
-	template<class D0, class D1>
-	static void scale(scalar alpha, const D0 & x, D1 & y) {
-		fordofs([&](len_t i) { y[i] = alpha * x[i]; }, x, y);
+	template<class Alpha, class D0, class D1>
+	static void scale(Alpha alpha, const D0 & x, D1 & y) {
+		const scalar value = scalar_value(std::move(alpha));
+		fordofs([&](len_t i) { y[i] = value * x[i]; }, x, y);
 	}
 
 	template<class X, class Y, class Z>
@@ -157,25 +160,36 @@ struct seq_ops {
 	template<class X>
 	static void dump(std::string_view pre, const X & x) {
 		std::string fname{pre};
-		fname += "-" + std::to_string(flecsi::run::context::instance().process());
+		fname +=
+			"-" + std::to_string(flecsi::run::context::instance().process());
 		std::ofstream ofile{fname};
 		fordofs([&](len_t i) { ofile << x[i] << '\n'; }, x);
 	}
 
-	template<class X, class Y, class Z>
+	template<class Alpha, class Beta, class X, class Y, class Z>
 	static void
-	linear_sum(scalar alpha, const X & x, scalar beta, const Y & y, Z & z) {
-		fordofs([&](len_t i) { z[i] = alpha * x[i] + beta * y[i]; }, x, y, z);
+	linear_sum(Alpha alpha, const X & x, Beta beta, const Y & y, Z & z) {
+		const scalar alpha_value = scalar_value(std::move(alpha));
+		const scalar beta_value = scalar_value(std::move(beta));
+		fordofs([&](len_t i) { z[i] = alpha_value * x[i] + beta_value * y[i]; },
+		        x,
+		        y,
+		        z);
 	}
 
-	template<class X, class Y, class Z>
-	static void axpy(scalar alpha, const X & x, const Y & y, Z & z) {
-		fordofs([&](len_t i) { z[i] = alpha * x[i] + y[i]; }, x, y, z);
+	template<class Alpha, class X, class Y, class Z>
+	static void axpy(Alpha alpha, const X & x, const Y & y, Z & z) {
+		const scalar value = scalar_value(std::move(alpha));
+		fordofs([&](len_t i) { z[i] = value * x[i] + y[i]; }, x, y, z);
 	}
 
-	template<class X, class Z>
-	static void axpby(scalar alpha, scalar beta, const X & x, Z & z) {
-		fordofs([&](len_t i) { z[i] = alpha * x[i] + beta * z[i]; }, x, z);
+	template<class Alpha, class Beta, class X, class Z>
+	static void axpby(Alpha alpha, Beta beta, const X & x, Z & z) {
+		const scalar alpha_value = scalar_value(std::move(alpha));
+		const scalar beta_value = scalar_value(std::move(beta));
+		fordofs([&](len_t i) { z[i] = alpha_value * x[i] + beta_value * z[i]; },
+		        x,
+		        z);
 	}
 
 	template<class X, class Y>
@@ -183,9 +197,10 @@ struct seq_ops {
 		fordofs([&](len_t i) { y[i] = std::abs(x[i]); }, x, y);
 	}
 
-	template<class X, class Y>
-	static void add_scalar(const X & x, scalar alpha, Y & y) {
-		fordofs([&](len_t i) { y[i] = x[i] + alpha; }, x, y);
+	template<class X, class Alpha, class Y>
+	static void add_scalar(const X & x, Alpha alpha, Y & y) {
+		const scalar value = scalar_value(std::move(alpha));
+		fordofs([&](len_t i) { y[i] = x[i] + value; }, x, y);
 	}
 
 	template<class X>
@@ -279,6 +294,16 @@ struct seq_ops {
 	}
 
 protected:
+	template<class T>
+	static scalar scalar_value(T value) {
+		return value;
+	}
+
+	template<class F, class T>
+	static scalar scalar_value(future_transform<future<T>, F> value) {
+		return value.get();
+	}
+
 	template<class X, class Y>
 	static scalar scalar_prod(const X & x, const Y & y) {
 		scalar res = 0.0;
diff --git a/flecsolve/vectors/test/flecsi_multivector.cc b/flecsolve/vectors/test/flecsi_multivector.cc
index db19ad9..12f4a74 100644
--- a/flecsolve/vectors/test/flecsi_multivector.cc
+++ b/flecsolve/vectors/test/flecsi_multivector.cc
@@ -152,7 +152,7 @@ int vectest(flecsi::scheduler & s) {
 		tmp.multiply(x, z);
 		EXPECT_TRUE(run(tmp, msh, mult));
 
-		x.add_scalar(x, 1);
+		x.add_scalar(x, flecsi::make_future(1));
 		EXPECT_TRUE(run(x, msh, scalar_add));
 
 		tmp.divide(y, x);
@@ -160,7 +160,10 @@ int vectest(flecsi::scheduler & s) {
 
 		x.add_scalar(x, -1);
 
-		tmp.scale(3.5, x);
+		tmp.scale(flecsi::make_future(3.5), x);
+		EXPECT_TRUE(run(tmp, msh, scale));
+		tmp.copy(x);
+		tmp.scale(flecsi::make_future(3.5));
 		EXPECT_TRUE(run(tmp, msh, scale));
 
 		y.add_scalar(y, 2);
@@ -168,7 +171,7 @@ int vectest(flecsi::scheduler & s) {
 		EXPECT_TRUE(run(tmp, msh, recip));
 		y.add_scalar(y, -2);
 
-		tmp.linear_sum(8, y, 9, z);
+		tmp.linear_sum(flecsi::make_future(8), y, flecsi::make_future(9), z);
 		EXPECT_TRUE(run(tmp, msh, linsum));
 
 		tmp.axpy(7, x, y);
@@ -178,7 +181,7 @@ int vectest(flecsi::scheduler & s) {
 		EXPECT_TRUE(run(tmp, msh, axpy));
 
 		tmp.copy(y);
-		tmp.axpby(4, 11, z);
+		tmp.axpby(flecsi::make_future(4), flecsi::make_future(11), z);
 		EXPECT_TRUE(run(tmp, msh, axpby));
 
 		tmp.add_scalar(y, -4.3);
diff --git a/flecsolve/vectors/test/flecsi_vector.cc b/flecsolve/vectors/test/flecsi_vector.cc
index b80724d..5b4e835 100644
--- a/flecsolve/vectors/test/flecsi_vector.cc
+++ b/flecsolve/vectors/test/flecsi_vector.cc
@@ -372,6 +372,10 @@ int vectest(flecsi::scheduler & s) {
 			tmp.multiply(x, z);
 			check_answer(mult_check{}, tmp);
 
+			tmp.add_scalar(x,
+			               flecsi::make_future(
+							   scalar_add_check::scalar_value<is_complex>()));
+			check_answer(scalar_add_check{}, tmp);
 			x.add_scalar(x, scalar_add_check::scalar_value<is_complex>());
 			check_answer(scalar_add_check{}, x);
 
@@ -379,7 +383,12 @@ int vectest(flecsi::scheduler & s) {
 			check_answer(div_check{}, tmp);
 
 			x.add_scalar(x, scale_check::shift_value<is_complex>());
-			tmp.scale(scale_check::scale_value<is_complex>(), x);
+			tmp.scale(
+				flecsi::make_future(scale_check::scale_value<is_complex>()), x);
+			check_answer(scale_check{}, tmp);
+			tmp.copy(x);
+			tmp.scale(
+				flecsi::make_future(scale_check::scale_value<is_complex>()));
 			check_answer(scale_check{}, tmp);
 
 			y.add_scalar(y, recip_check::shift_value<is_complex>());
@@ -387,7 +396,10 @@ int vectest(flecsi::scheduler & s) {
 			check_answer(recip_check{}, tmp);
 
 			y.add_scalar(y, linsum_check::shift_value<is_complex>());
-			tmp.linear_sum(linsum_check::alpha, y, linsum_check::beta, z);
+			tmp.linear_sum(flecsi::make_future(linsum_check::alpha),
+			               y,
+			               flecsi::make_future(linsum_check::beta),
+			               z);
 			check_answer(linsum_check{}, tmp);
 
 			tmp.axpy(axpy_check::alpha<is_complex>(), x, y);
@@ -398,8 +410,8 @@ int vectest(flecsi::scheduler & s) {
 			check_answer(axpy_check{}, tmp);
 
 			tmp.copy(y);
-			tmp.axpby(axpby_check::alpha<is_complex>(),
-			          axpby_check::beta<is_complex>(),
+			tmp.axpby(flecsi::make_future(axpby_check::alpha<is_complex>()),
+			          flecsi::make_future(axpby_check::beta<is_complex>()),
 			          z);
 			check_answer(axpby_check{}, tmp);
 
@@ -416,8 +428,8 @@ int vectest(flecsi::scheduler & s) {
 
 			auto test_reductions = [&](auto & a, auto & b) {
 				auto [alpha, beta] = dot_check::scalars<is_complex>();
-				a.set_scalar(alpha);
-				b.set_scalar(beta);
+				a.set_scalar(flecsi::make_future(alpha));
+				b.set_scalar(flecsi::make_future(beta));
 				EXPECT_LT(dot_check::abs_error<is_complex>(a.dot(b).get()),
 				          dot_check::ftol);
 

From 8c4c228710a5c77364bbfd304ac6de636eddd302 Mon Sep 17 00:00:00 2001
From: Andrew Reisner <areisner@lanl.gov>
Date: Thu, 27 Aug 2026 09:47:13 -0600
Subject: [PATCH 3/4] First pass at supporting futures in explicit time
 integration.

This separates the time step function from the explicit time
integrators allowing them to be used as a lightweight interface to the
time integration schemes.  Support for futures is also added to this
interface.

Assisted-by: gpt-5.4
---
 flecsolve/time-integrators/rk23.hh          | 108 +++++++++++----
 flecsolve/time-integrators/rk45.hh          | 146 ++++++++++++++------
 flecsolve/time-integrators/test/explicit.cc |  42 +++++-
 3 files changed, 221 insertions(+), 75 deletions(-)

diff --git a/flecsolve/time-integrators/rk23.hh b/flecsolve/time-integrators/rk23.hh
index 07c11cc..957981b 100644
--- a/flecsolve/time-integrators/rk23.hh
+++ b/flecsolve/time-integrators/rk23.hh
@@ -16,9 +16,14 @@ to do so.
 #ifndef FLECSI_LINALG_TIME_INTEGRATOR_RK23_H
 #define FLECSI_LINALG_TIME_INTEGRATOR_RK23_H
 
+#include <functional>
+#include <type_traits>
+#include <utility>
+
 #include "flecsolve/solvers/solver_settings.hh"
 #include "flecsolve/time-integrators/parameters.hh"
 #include "flecsolve/time-integrators/base.hh"
+#include "flecsolve/util/future.hh"
 
 namespace flecsolve::time_integrator::rk23 {
 
@@ -28,6 +33,9 @@ struct settings : base_settings {
 	bool use_fixed_dt;
 };
 
+struct stepper_settings {
+};
+
 struct options : base_options {
 	using settings_type = settings;
 	explicit options(const char * pre) : base_options(pre) {}
@@ -45,16 +53,40 @@ struct options : base_options {
 		return desc;
 	}
 };
-template<class Op, class Work>
-struct parameters : time_integrator::parameters<settings, Op, Work> {
-	using base = time_integrator::parameters<settings, Op, Work>;
-
-	template<class W>
-	parameters(const settings & s, op::handle<Op> op, W && work)
+template<class Op, class Work, class Settings = settings>
+struct parameters : time_integrator::parameters<Settings, Op, Work> {
+	using base = time_integrator::parameters<Settings, Op, Work>;
+
+	template<
+		class W,
+		class S = Settings,
+		std::enable_if_t<!std::is_same_v<S, stepper_settings>, bool> = true>
+	parameters(const S & s, op::handle<Op> op, W && work)
 		: base(s, op, std::forward<W>(work)) {}
+
+	template<class W,
+	         class S = Settings,
+	         std::enable_if_t<std::is_same_v<S, stepper_settings>, bool> = true>
+	parameters(op::handle<Op> op, W && work)
+		: base(stepper_settings{}, op, std::forward<W>(work)) {}
+
+	auto & get_work() { return unwrap_work(this->work); }
+
+private:
+	template<class T>
+	static T & unwrap_work(T & w) {
+		return w;
+	}
+
+	template<class T>
+	static T & unwrap_work(std::reference_wrapper<T> w) {
+		return w.get();
+	}
 };
 template<class O, class W>
 parameters(const settings &, op::handle<O>, W &&) -> parameters<O, W>;
+template<class O, class W>
+parameters(op::handle<O>, W &&) -> parameters<O, W, stepper_settings>;
 
 enum workvecs : std::size_t { k1, k2, k3, k4, z, next, nvecs };
 
@@ -63,46 +95,74 @@ template<std::size_t Version = 0>
 using topo_work = topo_work_base<workvecs::nvecs, Version>;
 
 template<class O, class W>
-struct integrator : base<parameters<O, W>> {
-	using P = parameters<O, W>;
-	using base<P>::params;
-	using base<P>::current_dt;
-	using base<P>::current_time;
-	using base<P>::assert_can_advance;
-
-	integrator(P p) : base<P>(std::move(p)), total_step_rejects(0) {}
-
-	template<class Curr, class Out>
-	void advance(double dt, Curr & curr, Out & out) {
-		assert_can_advance();
-		current_dt = dt;
+struct stepper {
+	using P = parameters<O, W, stepper_settings>;
+
+	stepper(P p) : params(std::move(p)) {}
+
+	template<class DeltaT,
+	         class Curr,
+	         class Out,
+	         std::enable_if_t<is_scalar_or_future_v<std::decay_t<DeltaT>,
+	                                                typename Curr::scalar>,
+	                          bool> = true>
+	void advance(DeltaT dt, Curr & curr, Out & out) {
 		auto & F = params.get_operator();
-		auto & [k1, k2, k3, k4, z, next] = params.work;
+		auto & [k1, k2, k3, k4, z, next] = params.get_work();
+
+		auto dt_1_2 = 0.5 * dt;
+		auto dt_3_4 = 0.75 * dt;
+		auto dt_1_9 = (1. / 9.) * dt;
+		auto dt_1_72 = (1. / 72.) * dt;
 
 		// k1 = f(tn, un)
 		F.apply(curr, k1);
 		// u* = un + k1 * dt/2
-		next.axpy(0.5 * dt, k1, curr);
+		next.axpy(std::move(dt_1_2), k1, curr);
 		// k2 = f(t+dt/2, u*)
 		F.apply(next, k2);
 		// u* = un + 0.75 *k2 * dt
-		next.axpy(0.75 * dt, k2, curr);
+		next.axpy(std::move(dt_3_4), k2, curr);
 		// k3 = f(t + 0.75dt, u*)
 		F.apply(next, k3);
 
 		next.linear_sum(2.0, k1, 3.0, k2);
 		next.axpy(4.0, k3, next);
-		next.axpy(dt / 9.0, next, curr);
+		next.axpy(std::move(dt_1_9), next, curr);
 
 		F.apply(next, k4);
 
 		z.linear_sum(-5., k1, 6., k2);
 		z.axpy(8., k3, z);
 		z.axpy(-9., k4, z);
-		z.scale(dt / 72.);
+		z.scale(std::move(dt_1_72));
 		out.copy(next);
 	}
 
+protected:
+	P params;
+};
+template<class O, class W>
+stepper(parameters<O, W, stepper_settings>) -> stepper<O, W>;
+
+template<class O, class W>
+struct integrator : base<parameters<O, W>> {
+	using P = parameters<O, W>;
+	using base<P>::params;
+	using base<P>::current_dt;
+	using base<P>::current_time;
+	using base<P>::assert_can_advance;
+
+	integrator(P p) : base<P>(std::move(p)), total_step_rejects(0) {}
+
+	template<class Curr, class Out>
+	void advance(double dt, Curr & curr, Out & out) {
+		assert_can_advance();
+		current_dt = dt;
+		stepper(parameters{params.op, std::ref(params.work)})
+			.advance(dt, curr, out);
+	}
+
 	bool check_solution() {
 		auto & z = std::get<workvecs::z>(params.work);
 
diff --git a/flecsolve/time-integrators/rk45.hh b/flecsolve/time-integrators/rk45.hh
index 8e7c341..405b575 100644
--- a/flecsolve/time-integrators/rk45.hh
+++ b/flecsolve/time-integrators/rk45.hh
@@ -16,27 +16,35 @@ to do so.
 #ifndef FLECSI_LINALG_TIME_INTEGRATOR_RK45_H
 #define FLECSI_LINALG_TIME_INTEGRATOR_RK45_H
 
+#include <functional>
+#include <type_traits>
+#include <utility>
+
 #include "flecsolve/solvers/solver_settings.hh"
 #include "flecsolve/time-integrators/parameters.hh"
 #include "flecsolve/time-integrators/rk23.hh"
 #include "flecsolve/time-integrators/base.hh"
+#include "flecsolve/util/future.hh"
 
 namespace flecsolve::time_integrator::rk45 {
 
-struct settings : rk23::settings {};
+struct settings : rk23::settings {
+};
 struct options : rk23::options {
 	using settings_type = settings;
 	explicit options(const char * pre) : rk23::options(pre) {}
 };
-template<class Op, class Work>
-struct parameters : rk23::parameters<Op, Work> {
-	template<class W>
-	parameters(const settings & s, op::handle<Op> op, W && work)
-		: rk23::parameters<Op, Work>(s, op,
-	                                 std::forward<W>(work)) {}
+using stepper_settings = rk23::stepper_settings;
+
+template<class Op, class Work, class Settings = settings>
+struct parameters : rk23::parameters<Op, Work, Settings> {
+	using base = rk23::parameters<Op, Work, Settings>;
+	using base::base;
 };
 template<class O, class W>
 parameters(const settings &, op::handle<O>, W &&) -> parameters<O, W>;
+template<class O, class W>
+parameters(op::handle<O>, W &&) -> parameters<O, W, stepper_settings>;
 
 enum workvecs : std::size_t { k1, k2, k3, k4, k5, k6, z, next, nvecs };
 
@@ -45,67 +53,115 @@ template<std::size_t Version = 0>
 using topo_work = topo_work_base<workvecs::nvecs, Version>;
 
 template<class O, class W>
-struct integrator : base<parameters<O, W>> {
-	using P = parameters<O, W>;
-	using base<P>::params;
-	using base<P>::current_dt;
-	using base<P>::current_time;
-	using base<P>::assert_can_advance;
-
-	integrator(P p) : base<P>(std::move(p)), total_step_rejects(0) {}
-
-	template<class Curr, class Out>
-	void advance(double dt, Curr & curr, Out & out) {
-		assert_can_advance();
-		current_dt = dt;
+struct stepper {
+	using P = parameters<O, W, stepper_settings>;
+
+	stepper(P p) : params(std::move(p)) {}
+
+	template<class DeltaT,
+	         class Curr,
+	         class Out,
+	         std::enable_if_t<
+				 is_scalar_or_future_v<std::decay_t<DeltaT>, typename Curr::scalar>,
+				 bool> = true>
+	void advance(DeltaT dt, Curr & curr, Out & out) {
 		auto & F = params.get_operator();
-		auto & [k1, k2, k3, k4, k5, k6, z, next] = params.work;
+		auto & [k1, k2, k3, k4, k5, k6, z, next] = params.get_work();
+
+		auto dt_1_4 = 0.25 * dt;
+		auto dt_3_32 = (3. / 32.) * dt;
+		auto dt_9_32 = (9. / 32.) * dt;
+		auto dt_1932_2197 = (1932. / 2197.) * dt;
+		auto dt_neg_7200_2197 = (-7200. / 2197.) * dt;
+		auto dt_7296_2197 = (7296. / 2197.) * dt;
+		auto dt_439_216 = (439. / 216.) * dt;
+		auto dt_neg_8 = -8. * dt;
+		auto dt_3680_513 = (3680. / 513.) * dt;
+		auto dt_neg_845_4104 = (-845. / 4104.) * dt;
+		auto dt_neg_8_27 = (-8. / 27.) * dt;
+		auto dt_2 = 2. * dt;
+		auto dt_neg_3544_2565 = (-3544. / 2565.) * dt;
+		auto dt_1859_4104 = (1859. / 4104.) * dt;
+		auto dt_neg_11_40 = (-11. / 40.) * dt;
+		auto dt_25_216 = (25. / 216.) * dt;
+		auto dt_1408_2565 = (1408. / 2565.) * dt;
+		auto dt_2197_4104 = (2197. / 4104.) * dt;
+		auto dt_neg_1_5 = -0.2 * dt;
+		auto dt_16_135 = (16. / 135.) * dt;
+		auto dt_6656_12825 = (6656. / 12825.) * dt;
+		auto dt_28561_56430 = (28561. / 56430.) * dt;
+		auto dt_neg_9_50 = (-9. / 50.) * dt;
+		auto dt_2_55 = (2. / 55.) * dt;
 
 		F.apply(curr, k1);
-		next.axpy(0.25 * dt, k1, curr);
+		next.axpy(std::move(dt_1_4), k1, curr);
 		F.apply(next, k2);
 
-		next.axpy(3.0 * dt / 32., k1, curr);
-		next.axpy(9.0 * dt / 32., k2, next);
+		next.axpy(std::move(dt_3_32), k1, curr);
+		next.axpy(std::move(dt_9_32), k2, next);
 
 		F.apply(next, k3);
 
-		next.axpy(1932. * dt / 2197., k1, curr);
-		next.axpy(-7200. * dt / 2197., k2, next);
-		next.axpy(7296. * dt / 2197., k3, next);
+		next.axpy(std::move(dt_1932_2197), k1, curr);
+		next.axpy(std::move(dt_neg_7200_2197), k2, next);
+		next.axpy(std::move(dt_7296_2197), k3, next);
 
 		F.apply(next, k4);
 
-		next.axpy(439. * dt / 216., k1, curr);
-		next.axpy(-8. * dt, k2, next);
-		next.axpy(3680. * dt / 513., k3, next);
-		next.axpy(-845. * dt / 4104., k4, next);
+		next.axpy(std::move(dt_439_216), k1, curr);
+		next.axpy(std::move(dt_neg_8), k2, next);
+		next.axpy(std::move(dt_3680_513), k3, next);
+		next.axpy(std::move(dt_neg_845_4104), k4, next);
 
 		F.apply(next, k5);
 
-		next.axpy(-8. * dt / 27., k1, curr);
-		next.axpy(2. * dt, k2, next);
-		next.axpy(-3544. * dt / 2565., k3, next);
-		next.axpy(1859. * dt / 4104., k4, next);
-		next.axpy(-11. * dt / 40., k5, next);
+		next.axpy(std::move(dt_neg_8_27), k1, curr);
+		next.axpy(std::move(dt_2), k2, next);
+		next.axpy(std::move(dt_neg_3544_2565), k3, next);
+		next.axpy(std::move(dt_1859_4104), k4, next);
+		next.axpy(std::move(dt_neg_11_40), k5, next);
 
 		F.apply(next, k6);
 
-		z.axpy(25. * dt / 216., k1, curr);
-		z.axpy(1408. * dt / 2565., k3, z);
-		z.axpy(2197. * dt / 4104., k4, z);
-		z.axpy(-0.2 * dt, k5, z);
+		z.axpy(std::move(dt_25_216), k1, curr);
+		z.axpy(std::move(dt_1408_2565), k3, z);
+		z.axpy(std::move(dt_2197_4104), k4, z);
+		z.axpy(std::move(dt_neg_1_5), k5, z);
 
-		next.axpy(16. * dt / 135., k1, curr);
-		next.axpy(6656. * dt / 12825., k3, next);
-		next.axpy(28561. * dt / 56430., k4, next);
-		next.axpy(-9. * dt / 50., k5, next);
-		next.axpy(2. * dt / 55., k6, next);
+		next.axpy(std::move(dt_16_135), k1, curr);
+		next.axpy(std::move(dt_6656_12825), k3, next);
+		next.axpy(std::move(dt_28561_56430), k4, next);
+		next.axpy(std::move(dt_neg_9_50), k5, next);
+		next.axpy(std::move(dt_2_55), k6, next);
 
 		z.subtract(next, z);
 		out.copy(next);
 	}
 
+protected:
+	P params;
+};
+template<class O, class W>
+stepper(parameters<O, W, stepper_settings>) -> stepper<O, W>;
+
+template<class O, class W>
+struct integrator : base<parameters<O, W>> {
+	using P = parameters<O, W>;
+	using base<P>::params;
+	using base<P>::current_dt;
+	using base<P>::current_time;
+	using base<P>::assert_can_advance;
+
+	integrator(P p) : base<P>(std::move(p)), total_step_rejects(0) {}
+
+	template<class Curr, class Out>
+	void advance(double dt, Curr & curr, Out & out) {
+		assert_can_advance();
+		current_dt = dt;
+		stepper(parameters{params.op, std::ref(params.work)})
+			.advance(dt, curr, out);
+	}
+
 	bool check_solution() {
 		auto & z = std::get<workvecs::z>(params.work);
 
diff --git a/flecsolve/time-integrators/test/explicit.cc b/flecsolve/time-integrators/test/explicit.cc
index ba79855..74b5750 100644
--- a/flecsolve/time-integrators/test/explicit.cc
+++ b/flecsolve/time-integrators/test/explicit.cc
@@ -13,7 +13,9 @@ namespace flecsolve {
 
 const flecsi::field<double>::definition<testmesh, testmesh::cells> xd, xnewd;
 
-struct parameters { double lambda; };
+struct parameters {
+	double lambda;
+};
 struct rate : op::base<parameters> {
 	using base = op::base<parameters>;
 	explicit rate(double l) : base{l} {}
@@ -40,8 +42,8 @@ int extest(flecsi::scheduler & s) {
 
 		auto [ti23_var, ti23_fixed] = std::apply(
 			[&](auto &&... s) {
-				return std::make_tuple(rk23::integrator(rk23::parameters(
-					                                        s, op::ref(F), rk23::make_work(x)))...);
+				return std::make_tuple(rk23::integrator(
+					rk23::parameters(s, op::ref(F), rk23::make_work(x)))...);
 			},
 			read_config("explicit.cfg",
 		                rk23::options("variable"),
@@ -49,8 +51,8 @@ int extest(flecsi::scheduler & s) {
 
 		auto [ti45_var, ti45_fixed] = std::apply(
 			[&](auto &&... s) {
-				return std::make_tuple(rk45::integrator(rk45::parameters(
-					                                        s, op::ref(F), rk45::make_work(x)))...);
+				return std::make_tuple(rk45::integrator(
+					rk45::parameters(s, op::ref(F), rk45::make_work(x)))...);
 			},
 			read_config("explicit.cfg",
 		                rk45::options("variable"),
@@ -69,7 +71,8 @@ int extest(flecsi::scheduler & s) {
 				dt = ti.get_next_dt(good_solution);
 			}
 
-			auto sol = ic * std::exp(F.get_params().lambda * ti.get_final_time());
+			auto sol =
+				ic * std::exp(F.get_params().lambda * ti.get_final_time());
 			auto approx = x.max().get();
 			return std::tuple(ti.get_final_time(),
 			                  std::abs(sol - approx),
@@ -97,6 +100,33 @@ int extest(flecsi::scheduler & s) {
 			EXPECT_LT(err, 1e-6);
 			EXPECT_EQ(step, 6);
 		}
+		{
+			auto dt = 0.05;
+
+			rk23::stepper rk23_double(
+				rk23::parameters(op::ref(F), rk23::make_work(x)));
+			x.set_scalar(ic);
+			rk23_double.advance(dt, x, xnew);
+			auto rk23_double_ans = xnew.max().get();
+
+			rk23::stepper rk23_future(
+				rk23::parameters(op::ref(F), rk23::make_work(x)));
+			x.set_scalar(ic);
+			rk23_future.advance(flecsi::make_future(dt), x, xnew);
+			EXPECT_EQ(xnew.max().get(), rk23_double_ans);
+
+			rk45::stepper rk45_double(
+				rk45::parameters(op::ref(F), rk45::make_work(x)));
+			x.set_scalar(ic);
+			rk45_double.advance(dt, x, xnew);
+			auto rk45_double_ans = xnew.max().get();
+
+			rk45::stepper rk45_future(
+				rk45::parameters(op::ref(F), rk45::make_work(x)));
+			x.set_scalar(ic);
+			rk45_future.advance(flecsi::make_future(dt), x, xnew);
+			EXPECT_EQ(xnew.max().get(), rk45_double_ans);
+		}
 	};
 }
 

From 2cd6c1e9de68ac541058c4af03f9b84ef3cdab75 Mon Sep 17 00:00:00 2001
From: Andrew Reisner <areisner@lanl.gov>
Date: Fri, 11 Sep 2026 10:44:53 -0600
Subject: [PATCH 4/4] Separate time step function in implicit time integrators.

Assisted-by: gpt-5.4
---
 flecsolve/time-integrators/bdf.hh            | 226 ++++++++++++-------
 flecsolve/time-integrators/bdf_parameters.hh |  94 +++++---
 flecsolve/time-integrators/test/implicit.cc  |  19 +-
 flecsolve/time-integrators/test/implicit.cfg |  12 +
 4 files changed, 230 insertions(+), 121 deletions(-)

diff --git a/flecsolve/time-integrators/bdf.hh b/flecsolve/time-integrators/bdf.hh
index 4709fd7..6ab41c8 100644
--- a/flecsolve/time-integrators/bdf.hh
+++ b/flecsolve/time-integrators/bdf.hh
@@ -50,18 +50,12 @@ static inline work_factory<work_size> make_work;
 
 short order(method meth);
 
-template<class O, class W, class S>
-struct integrator : base<parameters<O, W, S>> {
-	using P = parameters<O, W, S>;
-	using base<P>::params;
-	using base<P>::current_dt;
-	using base<P>::old_dt;
-	using base<P>::current_time;
-	using base<P>::assert_can_advance;
-	using base<P>::integrator_step;
+template<class O, class W, class S, class Settings = stepper_settings>
+struct stepper {
+	using P = parameters<O, W, S, Settings>;
 
-	integrator(P p)
-		: base<P>(std::move(p)),
+	stepper(P p)
+		: params(std::move(p)), current_dt{0}, old_dt{0}, integrator_step{0},
 		  prev(memory_size(params.integrator), std::ref(params.work)),
 		  solver_success{false}, total_steprejects{0} {
 		params.validate();
@@ -74,14 +68,13 @@ struct integrator : base<parameters<O, W, S>> {
 	void advance(double dt, bool first_step, Curr & curr, Out & out) {
 		flog_assert(curr != out,
 		            "BDF integrator: curr cannot be the same as out");
-		assert_can_advance();
 
 		prev.seat(std::ref(params.work));
 
 		prev[0].solution.copy(curr);
 		current_dt = dt;
 
-		set_initial_guess(first_step, current_time, current_dt, old_dt);
+		set_initial_guess(first_step, current_dt, old_dt);
 
 		auto & rhs = getvec<workvec::rhs>();
 		auto & source = getvec<workvec::source>();
@@ -128,8 +121,6 @@ struct integrator : base<parameters<O, W, S>> {
 	}
 
 	void update() {
-		current_time += current_dt;
-		this->new_time = current_time;
 		if (params.use_predictor) {
 			if (params.predictor == predictor::ab2) {
 				auto & old_td = getvec<workvec::old_td>();
@@ -145,76 +136,10 @@ struct integrator : base<parameters<O, W, S>> {
 		++integrator_step;
 	}
 
-	double get_next_dt(bool good_solution) {
-		double tmp_dt = current_dt;
-		if (params.timestep_strategy == strategy::truncation_error) {
-			if (good_solution) {
-				current_dt = estimate_dt_with_truncation_error_estimates(
-					current_dt, good_solution);
-
-				if (params.max_dt < current_dt)
-					current_dt = params.max_dt;
-			}
-			else {
-				if (solver_success) {
-					current_dt = estimate_dt_with_truncation_error_estimates(
-						current_dt, good_solution);
-				}
-				else {
-					current_dt = params.dt_cut_lower_bound * current_dt;
-				}
-			}
-			evaluate_predictor();
-			auto & predict = getvec<workvec::predict>();
-			auto & op = params.get_operator();
-			bool valid_vector = op.is_valid(predict);
-			if (!valid_vector) {
-				int number_of_predictor_precheck_events = 10;
-				for (int i = 0; i < number_of_predictor_precheck_events; ++i) {
-					current_dt = 0.5 * current_dt;
-					evaluate_predictor();
-					valid_vector = op.is_valid(predict);
-				}
-			}
-		}
-		else {
-			if (good_solution) {
-				if (params.timestep_strategy == strategy::constant)
-					current_dt = params.initial_dt;
-				else if (params.timestep_strategy == strategy::final_constant) {
-					static int i = 1;
-					if (i < params.number_of_time_intervals) {
-						current_dt =
-							params.initial_dt +
-							((double)i) /
-								((double)params.number_of_time_intervals) *
-								(params.max_dt - params.initial_dt);
-						++i;
-					}
-					else {
-						current_dt = params.max_dt;
-					}
-				}
-				else if (params.timestep_strategy ==
-				         strategy::limit_relative_change) {
-					current_dt = estimate_dynamical_time_scale(current_dt);
-				}
-			}
-			else
-				current_dt = 0.;
-		}
-
-		// now set old dt once it has been used
-		if (good_solution)
-			old_dt = tmp_dt;
-
-		current_dt = std::min(std::min(current_dt, params.max_dt),
-		                      params.final_time - current_time);
-		return current_dt;
-	}
-
 	int num_step_rejects() const { return total_steprejects; }
 
+	int get_current_step() const { return integrator_step; }
+
 protected:
 	method get_current_integrator() const {
 		if (integrator_step == 0) {
@@ -229,10 +154,8 @@ protected:
 	}
 
 	void set_initial_guess(bool first_step,
-	                       double current_time,
 	                       double current_dt,
 	                       double old_dt) {
-		(void)current_time;
 		(void)current_dt;
 		(void)old_dt;
 
@@ -860,6 +783,11 @@ protected:
 
 	struct history;
 
+	P params;
+	double current_dt;
+	double old_dt;
+	int integrator_step;
+
 	bool prev_successive_rejects;
 	int current_steprejects;
 	double new_time;
@@ -878,10 +806,134 @@ protected:
 	int total_steprejects;
 };
 template<class O, class W, class S>
-integrator(parameters<O, W, S>) -> integrator<O, W, S>;
+stepper(parameters<O, W, S, stepper_settings>) -> stepper<O, W, S>;
+
+template<class O, class W, class S>
+struct integrator : stepper<O, W, S, settings> {
+	using P = parameters<O, W, S>;
+	using base = stepper<O, W, S, settings>;
+
+	integrator(P p)
+		: base(std::move(p)), current_time(this->params.initial_time),
+		  max_integrator_steps(this->params.max_steps) {
+		this->current_dt = this->params.initial_dt;
+		this->old_dt = this->params.initial_dt;
+	}
+
+	double get_current_time() const { return current_time; }
+
+	double get_final_time() const { return this->params.final_time; }
+
+	int get_current_step() const { return this->integrator_step; }
 
+	double get_current_dt() const { return this->current_dt; }
+
+	bool fixed_dt() const { return false; }
+
+	bool steps_remaining() const {
+		return this->integrator_step < max_integrator_steps;
+	}
+
+	template<class Curr, class Out>
+	void advance(double dt, bool first_step, Curr & curr, Out & out) {
+		assert_can_advance();
+		base::advance(dt, first_step, curr, out);
+	}
+
+	void update() {
+		current_time += this->current_dt;
+		this->new_time = current_time;
+		base::update();
+	}
+
+	double get_next_dt(bool good_solution) {
+		double tmp_dt = this->current_dt;
+		if (this->params.timestep_strategy == strategy::truncation_error) {
+			if (good_solution) {
+				this->current_dt =
+					this->estimate_dt_with_truncation_error_estimates(
+						this->current_dt, good_solution);
+
+				if (this->params.max_dt < this->current_dt)
+					this->current_dt = this->params.max_dt;
+			}
+			else {
+				if (this->solver_success) {
+					this->current_dt =
+						this->estimate_dt_with_truncation_error_estimates(
+							this->current_dt, good_solution);
+				}
+				else {
+					this->current_dt =
+						this->params.dt_cut_lower_bound * this->current_dt;
+				}
+			}
+			this->evaluate_predictor();
+			auto & predict = this->template getvec<workvec::predict>();
+			auto & op = this->params.get_operator();
+			bool valid_vector = op.is_valid(predict);
+			if (!valid_vector) {
+				int number_of_predictor_precheck_events = 10;
+				for (int i = 0; i < number_of_predictor_precheck_events; ++i) {
+					this->current_dt = 0.5 * this->current_dt;
+					this->evaluate_predictor();
+					valid_vector = op.is_valid(predict);
+				}
+			}
+		}
+		else {
+			if (good_solution) {
+				if (this->params.timestep_strategy == strategy::constant)
+					this->current_dt = this->params.initial_dt;
+				else if (this->params.timestep_strategy == strategy::final_constant) {
+					static int i = 1;
+					if (i < this->params.number_of_time_intervals) {
+						this->current_dt =
+							this->params.initial_dt +
+							((double)i) /
+								((double)this->params.number_of_time_intervals) *
+								(this->params.max_dt - this->params.initial_dt);
+						++i;
+					}
+					else {
+						this->current_dt = this->params.max_dt;
+					}
+				}
+				else if (this->params.timestep_strategy ==
+				         strategy::limit_relative_change) {
+					this->current_dt =
+						this->estimate_dynamical_time_scale(this->current_dt);
+				}
+			}
+			else
+				this->current_dt = 0.;
+		}
+
+		// now set old dt once it has been used
+		if (good_solution)
+			this->old_dt = tmp_dt;
+
+		this->current_dt =
+			std::min(std::min(this->current_dt, this->params.max_dt),
+		             this->params.final_time - current_time);
+		return this->current_dt;
+	}
+
+protected:
+	void assert_can_advance() {
+		flog_assert(steps_remaining() &&
+		                (current_time < this->params.final_time),
+		            "Time integrator: already finished integrating");
+	}
+
+	double current_time;
+	int max_integrator_steps;
+};
 template<class O, class W, class S>
-struct integrator<O, W, S>::history {
+integrator(parameters<O, W, S>) -> integrator<O, W, S>;
+
+template<class O, class W, class S, class Settings>
+struct stepper<O, W, S, Settings>::history {
 	struct vec_arr {
 		using size_type = std::size_t;
 		using vec_t = typename std::remove_reference_t<W>::value_type;
diff --git a/flecsolve/time-integrators/bdf_parameters.hh b/flecsolve/time-integrators/bdf_parameters.hh
index 790efe4..789cb06 100644
--- a/flecsolve/time-integrators/bdf_parameters.hh
+++ b/flecsolve/time-integrators/bdf_parameters.hh
@@ -46,7 +46,7 @@ std::istream & operator>>(std::istream &, method &);
 std::istream & operator>>(std::istream &, controller &);
 std::istream & operator>>(std::istream &, error_scaling &);
 
-struct settings : base_settings {
+struct stepper_settings {
 	bool use_predictor;
 	bool use_initial_predictor;
 	bool has_source_term;
@@ -112,50 +112,77 @@ struct settings : base_settings {
 	}
 };
 
+struct settings : base_settings, stepper_settings {};
+
+template<class Options>
+auto add_stepper_options(Options & opts, stepper_settings & s) {
+	po::options_description desc;
+	// clang-format off
+	desc.add_options()
+		(opts.option_label("use-predictor").c_str(), po::value<bool>(&s.use_predictor)->default_value(true), "use a predictor")
+		(opts.option_label("has-source-term").c_str(), po::value<bool>(&s.has_source_term)->default_value(false), "has source term")
+		(opts.option_label("use-initial-predictor").c_str(), po::value<bool>(&s.use_initial_predictor)->default_value(true), "use an initial predictor")
+		(opts.option_label("predictor").c_str(), po::value<bdf::predictor>(&s.predictor)->required(), "Predictor type (AB2 or Leapfrog)")
+		(opts.option_label("timestep-selection-strategy").c_str(), po::value<strategy>(&s.timestep_strategy)->required(),
+		 "Timestep selection strategy (truncation error or constant or final constant or limit relative change)")
+		// these bounds are based on the paper by Emmrich, 2008 for nonlinear evolution
+		(opts.option_label("dt-cut-lower-bound").c_str(), po::value<double>(&s.dt_cut_lower_bound)->default_value(0.58754407), "dt-cut-lower-bound")
+		(opts.option_label("dt-growth-upper-bound").c_str(), po::value<double>(&s.dt_growth_upper_bound)->default_value(1.702), "dt-growth-upper-bound")
+		(opts.option_label("number-of-time-intervals").c_str(), po::value<int>(&s.number_of_time_intervals)->default_value(100), "Number of time intervals (final constant strategy)")
+		(opts.option_label("integrator").c_str(), po::value<bdf::method>(&s.integrator)->required(), "Implicit integrator (BE or CN or BDF{2-6}")
+		(opts.option_label("starting-integrator").c_str(), po::value<bdf::method>(&s.starting_integrator)->required()->notifier([](const bdf::method & m) {
+			flog_assert(memory_size(m) == 1, "BE or CN must be used for starting integrator");}), "one step integrator (BE or CN)")
+		(opts.option_label("calculate-time-trunc-error").c_str(), po::value<bool>(&s.calculate_time_trunc_error), "Whether to calculate the time truncation error")
+		(opts.option_label("time-trunc-error-norm").c_str(), po::value<vec::norm_type>(&s.time_trunc_err_norm), "Norm time used for time truncation error")
+		(opts.option_label("target-relative-change").c_str(), po::value<double>(&s.target_relative_change), "Relative change to target in time scale")
+		(opts.option_label("use-pi-controller").c_str(), po::value<bool>(&s.use_pi_controller)->default_value(true), "Whether to use the PI controller")
+		(opts.option_label("pi-controller-type").c_str(), po::value<bdf::controller>(&s.pi_controller_type)->default_value(bdf::controller::pc4_7, "PC.4.7"), "Type of PI controller")
+		(opts.option_label("control-timestep-variation").c_str(), po::value<bool>(&s.control_timestep_variation)->default_value(false), "Control timestep variation")
+		(opts.option_label("time-error-scaling").c_str(), po::value<bdf::error_scaling>(&s.time_error_scaling)->default_value(error_scaling::fixed_scaling, "fixed-scaling"), "Time error scaling")
+		(opts.option_label("trunc-error-rtol").c_str(), po::value<double>(&s.time_rtol)->default_value(1e-9), "Relative tolerance for truncation error")
+		(opts.option_label("trunc-error-atol").c_str(), po::value<double>(&s.time_atol)->default_value(1e-15), "Absolute tolerance for truncation error")
+		(opts.option_label("combine-timestep-estimators").c_str(), po::value<bool>(&s.combine_timestep_estimators)->default_value(false), "Combine timestep estimators")
+		(opts.option_label("problem-scales").c_str(), po::value<std::vector<double>>(&s.problem_scales)->multitoken(), "Fixed scaling for problem");
+	// clang-format on
+
+	return desc;
+}
+
+struct stepper_options : with_label {
+	using settings_type = stepper_settings;
+	explicit stepper_options(const char * pre) : with_label(pre) {}
+
+	std::string option_label(const char * suffix) {
+		return label(suffix);
+	}
+
+	auto operator()(stepper_settings & s) {
+		return add_stepper_options(*this, s);
+	}
+};
+
 struct options : base_options {
 	using settings_type = settings;
 	explicit options(const char * pre) : base_options(pre) {}
 
+	std::string option_label(const char * suffix) {
+		return label(suffix);
+	}
+
 	auto operator()(settings & s) {
 		auto desc = base_options::operator()(s);
-		// clang-format off
-		desc.add_options()
-			(label("use-predictor").c_str(), po::value<bool>(&s.use_predictor)->default_value(true), "use a predictor")
-			(label("has-source-term").c_str(), po::value<bool>(&s.has_source_term)->default_value(false), "has source term")
-			(label("use-initial-predictor").c_str(), po::value<bool>(&s.use_initial_predictor)->default_value(true), "use an initial predictor")
-			(label("predictor").c_str(), po::value<bdf::predictor>(&s.predictor)->required(), "Predictor type (AB2 or Leapfrog)")
-			(label("timestep-selection-strategy").c_str(), po::value<strategy>(&s.timestep_strategy)->required(),
-			 "Timestep selection strategy (truncation error or constant or final constant or limit relative change)")
-			// these bounds are based on the paper by Emmrich, 2008 for nonlinear evolution
-			(label("dt-cut-lower-bound").c_str(), po::value<double>(&s.dt_cut_lower_bound)->default_value(0.58754407), "dt-cut-lower-bound")
-			(label("dt-growth-upper-bound").c_str(), po::value<double>(&s.dt_growth_upper_bound)->default_value(1.702), "dt-growth-upper-bound")
-			(label("number-of-time-intervals").c_str(), po::value<int>(&s.number_of_time_intervals)->default_value(100), "Number of time intervals (final constant strategy)")
-			(label("integrator").c_str(), po::value<bdf::method>(&s.integrator)->required(), "Implicit integrator (BE or CN or BDF{2-6}")
-			(label("starting-integrator").c_str(), po::value<bdf::method>(&s.starting_integrator)->required()->notifier([](const bdf::method & m) {
-				flog_assert(memory_size(m) == 1, "BE or CN must be used for starting integrator");}), "one step integrator (BE or CN)")
-			(label("calculate-time-trunc-error").c_str(), po::value<bool>(&s.calculate_time_trunc_error), "Whether to calculate the time truncation error")
-			(label("time-trunc-error-norm").c_str(), po::value<vec::norm_type>(&s.time_trunc_err_norm), "Norm time used for time truncation error")
-			(label("target-relative-change").c_str(), po::value<double>(&s.target_relative_change), "Relative change to target in time scale")
-			(label("use-pi-controller").c_str(), po::value<bool>(&s.use_pi_controller)->default_value(true), "Whether to use the PI controller")
-			(label("pi-controller-type").c_str(), po::value<bdf::controller>(&s.pi_controller_type)->default_value(bdf::controller::pc4_7, "PC.4.7"), "Type of PI controller")
-			(label("control-timestep-variation").c_str(), po::value<bool>(&s.control_timestep_variation)->default_value(false), "Control timestep variation")
-			(label("time-error-scaling").c_str(), po::value<bdf::error_scaling>(&s.time_error_scaling)->default_value(error_scaling::fixed_scaling, "fixed-scaling"), "Time error scaling")
-			(label("trunc-error-rtol").c_str(), po::value<double>(&s.time_rtol)->default_value(1e-9), "Relative tolerance for truncation error")
-			(label("trunc-error-atol").c_str(), po::value<double>(&s.time_atol)->default_value(1e-15), "Absolute tolerance for truncation error")
-			(label("combine-timestep-estimators").c_str(), po::value<bool>(&s.combine_timestep_estimators)->default_value(false), "Combine timestep estimators")
-			(label("problem-scales").c_str(), po::value<std::vector<double>>(&s.problem_scales)->multitoken(), "Fixed scaling for problem");
-		// clang-format on
+		desc.add(add_stepper_options(*this, s));
 
 		return desc;
 	}
 };
 
-template<class Op, class Work, class Solver>
-struct parameters : time_integrator::parameters<settings, Op, Work> {
-	using base = time_integrator::parameters<settings, Op, Work>;
+template<class Op, class Work, class Solver, class Settings = settings>
+struct parameters : time_integrator::parameters<Settings, Op, Work> {
+	using base = time_integrator::parameters<Settings, Op, Work>;
 
 	template<class O, class W, class S>
-	parameters(const settings & s, op::handle<O> op, W && work, op::handle<S> solver)
+	parameters(const Settings & s, op::handle<O> op, W && work, op::handle<S> solver)
 		: base(s, op, std::forward<W>(work)),
 		  solver(solver) {}
 
@@ -169,5 +196,8 @@ protected:
 
 template<class O, class W, class S>
 parameters(const settings &, op::handle<O>, W &&, op::handle<S>) -> parameters<O, W, S>;
+template<class O, class W, class S>
+parameters(const stepper_settings &, op::handle<O>, W &&, op::handle<S>)
+	-> parameters<O, W, S, stepper_settings>;
 }
 #endif
diff --git a/flecsolve/time-integrators/test/implicit.cc b/flecsolve/time-integrators/test/implicit.cc
index c2a0a64..13246f7 100644
--- a/flecsolve/time-integrators/test/implicit.cc
+++ b/flecsolve/time-integrators/test/implicit.cc
@@ -86,13 +86,17 @@ int bdftest(flecsi::scheduler & s) {
 		auto xnew = vec::make(msh)(xnewd);
 
 		auto solver = op::make_shared<rate_solver>(F);
+		auto [bdf2_settings, bdf5_settings, stepper_settings] = read_config(
+			"implicit.cfg",
+			bdf::options("bdf-2"),
+			bdf::options("bdf-5"),
+			bdf::stepper_options("bdf-stepper"));
 		auto [ti2, ti5] = std::apply(
 			[&](auto &&... s) {
 				return std::make_tuple(bdf::integrator(bdf::parameters(
 					                                       s, F, bdf::make_work(x), solver))...);
 			},
-			read_config(
-				"implicit.cfg", bdf::options("bdf-2"), bdf::options("bdf-5")));
+			std::tie(bdf2_settings, bdf5_settings));
 
 		auto run = [&](auto & ti) {
 			x.set_scalar(ic);
@@ -132,6 +136,17 @@ int bdftest(flecsi::scheduler & s) {
 			EXPECT_EQ(std::get<2>(info), 48);
 			EXPECT_EQ(std::get<3>(info), 14);
 		}
+		{
+			bdf::stepper ti(
+				bdf::parameters(stepper_settings, F, bdf::make_work(x), solver));
+			x.set_scalar(ic);
+			ti.advance(0.5, true, x, xnew);
+			ti.update();
+
+			EXPECT_EQ(ti.get_current_step(), 1);
+			EXPECT_EQ(ti.num_step_rejects(), 0);
+			EXPECT_LT(std::abs(xnew.max().get() - 2.), 1e-12);
+		}
 	};
 }
 
diff --git a/flecsolve/time-integrators/test/implicit.cfg b/flecsolve/time-integrators/test/implicit.cfg
index 904bf11..7e83ca2 100644
--- a/flecsolve/time-integrators/test/implicit.cfg
+++ b/flecsolve/time-integrators/test/implicit.cfg
@@ -35,3 +35,15 @@ problem-scales = 1.
 trunc-error-rtol = 1e-7
 trunc-error-atol = 1e-5
 time-error-scaling = fixed-resolution
+
+[bdf-stepper]
+timestep-selection-strategy = constant
+integrator = BE
+starting-integrator = BE
+calculate-time-trunc-error = false
+time-trunc-error-norm = l2
+predictor = leapfrog
+use-predictor = false
+use-initial-predictor = false
+target-relative-change = 0.
+time-error-scaling = fixed-scaling