diff --git a/examples/symmetric_laplacian.cpp b/examples/symmetric_laplacian.cpp
index ecf9212..2558994 100644
--- a/examples/symmetric_laplacian.cpp
+++ b/examples/symmetric_laplacian.cpp
@@ -20,6 +20,7 @@ int main() {
     std::vector<double> resid(n, 0.0);
     std::vector<double> v(static_cast<std::size_t>(n) * static_cast<std::size_t>(ncv), 0.0);
     std::vector<double> workd(3 * n, 0.0);
+    std::vector<double> workc(ncv, 0.0);  // projection scratch (length up to ncv)
     std::vector<double> workl(lworkl, 0.0);
     std::vector<double> d(nev, 0.0);
     std::vector<double> z(static_cast<std::size_t>(n) * static_cast<std::size_t>(nev), 0.0);
@@ -36,7 +37,7 @@ int main() {
     auto op = [&](const double* x, double* y) { laplacian_matvec(n, x, y); };
 
     arnoldi::detail::saupd<double>("I", n, "SM", nev, tol, resid.data(), ncv,
-                              v.data(), n, iparam, ipntr, workd.data(),
+                              v.data(), n, iparam, ipntr, workd.data(), workc.data(),
                               workl.data(), lworkl, info, op);
 
     if (info < 0) {
diff --git a/include/arnoldi/arnoldi.hpp b/include/arnoldi/arnoldi.hpp
index 725cdbf..64184cf 100644
--- a/include/arnoldi/arnoldi.hpp
+++ b/include/arnoldi/arnoldi.hpp
@@ -106,6 +106,7 @@ namespace arnoldi {
       resid_.assign(static_cast<std::size_t>(n_), Scalar{});
       v_.assign(static_cast<std::size_t>(ldv_) * ncv_, Scalar{});
       workd_.assign(static_cast<std::size_t>(3 * n_), Scalar{});
+      workc_.assign(static_cast<std::size_t>(ncv_), Scalar{});
 
       if constexpr (K == Kind::Nonsym) {
         lworkl_ = 3 * ncv_ * ncv_ + 6 * ncv_;
@@ -253,12 +254,14 @@ namespace arnoldi {
       auto bref = backend_.ref();
       if constexpr (K == Kind::Nonsym) {
         arnoldi::detail::naupd<real_type, Backend>(bref, bmat_.c_str(), n_, which_.c_str(), nev_, tol_, resid_.data(), ncv_,
-                                                   v_.data(), ldv_, iparam_.data(), ipntr_.data(), workd_.data(), workl_.data(),
-                                                   lworkl_, info_, std::forward<Op>(op), std::forward<Bop>(bop), comm_);
+                                                   v_.data(), ldv_, iparam_.data(), ipntr_.data(), workd_.data(),
+                                                   workc_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op),
+                                                   std::forward<Bop>(bop), comm_);
       } else {  // Sym or Herm
         arnoldi::detail::saupd<Scalar, Backend>(bref, bmat_.c_str(), n_, which_.c_str(), nev_, tol_, resid_.data(), ncv_,
-                                                v_.data(), ldv_, iparam_.data(), ipntr_.data(), workd_.data(), workl_.data(),
-                                                lworkl_, info_, std::forward<Op>(op), std::forward<Bop>(bop), comm_);
+                                                v_.data(), ldv_, iparam_.data(), ipntr_.data(), workd_.data(),
+                                                workc_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op),
+                                                std::forward<Bop>(bop), comm_);
       }
     }
 
@@ -268,13 +271,13 @@ namespace arnoldi {
       if constexpr (K == Kind::Nonsym) {
         arnoldi::detail::naupd<real_type, Backend>(
             bref, bmat_.c_str(), n_, which_.c_str(), nev_, tol_, resid_.data(), ncv_, v_.data(), ldv_, iparam_.data(),
-            ipntr_.data(), workd_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op),
+            ipntr_.data(), workd_.data(), workc_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op),
             [](const real_type*, real_type*) {}, comm_);
       } else {
         arnoldi::detail::saupd<Scalar, Backend>(
             bref, bmat_.c_str(), n_, which_.c_str(), nev_, tol_, resid_.data(), ncv_, v_.data(), ldv_, iparam_.data(),
-            ipntr_.data(), workd_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op), [](const Scalar*, Scalar*) {},
-            comm_);
+            ipntr_.data(), workd_.data(), workc_.data(), workl_.data(), lworkl_, info_, std::forward<Op>(op),
+            [](const Scalar*, Scalar*) {}, comm_);
       }
     }
 
@@ -339,6 +342,9 @@ namespace arnoldi {
     typename arnoldi::detail::buffer_traits<Backend, scalar_type>::vector_type resid_;
     typename arnoldi::detail::buffer_traits<Backend, scalar_type>::vector_type v_;
     typename arnoldi::detail::buffer_traits<Backend, scalar_type>::vector_type workd_;
+    // ncv-sized Scalar scratch for the Gram-Schmidt / B-orthogonalization
+    // projection coefficients (length up to kplusp).
+    typename arnoldi::detail::buffer_traits<Backend, scalar_type>::vector_type workc_;
     // The IRA bookkeeping (Hessenberg, Ritz values/bounds, Givens scratch)
     // is small (ncv²-sized) and stays host-resident under every backend.
     std::vector<real_type>   workl_;
diff --git a/include/arnoldi/detail/getv0.hpp b/include/arnoldi/detail/getv0.hpp
index d2405a2..06760f7 100644
--- a/include/arnoldi/detail/getv0.hpp
+++ b/include/arnoldi/detail/getv0.hpp
@@ -12,7 +12,8 @@ namespace arnoldi::detail {
   // invoked directly instead of returning ido.
   template <typename Scalar, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void getv0(detail::BackendRef<Backend> bref, const char* bmat, int itry, bool initv, int n, int j, Scalar* v, int ldv,
-             Scalar* resid, detail::real_t<Scalar>& rnorm, Scalar* workd, int& ierr, OP&& op, BOP&& bop, const Comm& comm) {
+             Scalar* resid, detail::real_t<Scalar>& rnorm, Scalar* workd, Scalar* workc, int& ierr, OP&& op, BOP&& bop,
+             const Comm& comm) {
     using Real          = detail::real_t<Scalar>;
     static int iseed[4] = {1, 3, 5, 7};
     int        msglvl   = detail::debug.getv0;
@@ -50,9 +51,9 @@ namespace arnoldi::detail {
     if (j != 1) {
       for (int iter = 0;; ++iter) {
         detail::Ops<Scalar, Backend>::gemv(bref, detail::Ops<Scalar, Backend>::herm_trans(), n, j - 1, Scalar(1), v, ldv, workd, 1,
-                                           Scalar(0), &workd[n], 1);
-        detail::backend_allreduce(bref, comm, &workd[n], j - 1);
-        detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j - 1, Scalar(-1), v, ldv, &workd[n], 1, Scalar(1), resid, 1);
+                                           Scalar(0), workc, 1);
+        detail::backend_allreduce(bref, comm, workc, j - 1);
+        detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j - 1, Scalar(-1), v, ldv, workc, 1, Scalar(1), resid, 1);
 
         detail::arscnd(t2);
         rnorm = detail::bnorm<Scalar, Backend>(bref, *bmat, n, resid, workd, &workd[n], bop, comm);
@@ -85,12 +86,12 @@ namespace arnoldi::detail {
   }
 
   // getv0 without explicit BackendRef — defaults to CpuBackend, for the
-  // historical public callback API (tests, examples).
+  // public callback API (tests, examples).
   template <typename Scalar, typename OP, typename BOP, typename Comm>
   void getv0(const char* bmat, int itry, bool initv, int n, int j, Scalar* v, int ldv, Scalar* resid,
-             detail::real_t<Scalar>& rnorm, Scalar* workd, int& ierr, OP&& op, BOP&& bop, const Comm& comm) {
+             detail::real_t<Scalar>& rnorm, Scalar* workd, Scalar* workc, int& ierr, OP&& op, BOP&& bop, const Comm& comm) {
     getv0<Scalar, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, itry, initv, n, j, v, ldv, resid, rnorm,
-                                      workd, ierr, std::forward<OP>(op), std::forward<BOP>(bop), comm);
+                                      workd, workc, ierr, std::forward<OP>(op), std::forward<BOP>(bop), comm);
   }
 
 }  // namespace arnoldi::detail
diff --git a/include/arnoldi/detail/nonsym.hpp b/include/arnoldi/detail/nonsym.hpp
index 13f9015..a780391 100644
--- a/include/arnoldi/detail/nonsym.hpp
+++ b/include/arnoldi/detail/nonsym.hpp
@@ -21,7 +21,7 @@ namespace arnoldi::detail {
   // naitr — Arnoldi iteration (nonsymmetric).
   template <typename Real, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void naitr(detail::BackendRef<Backend> bref, const char* bmat, int n, int k, int np, int nb, Real* resid, Real& rnorm, Real* v,
-             int ldv, Real* h, int ldh, Real* workd, int& info, OP&& op, BOP&& bop, const Comm& comm) {
+             int ldv, Real* h, int ldh, Real* workd, Real* workc, int& info, OP&& op, BOP&& bop, const Comm& comm) {
     const int ipj = 0, irj = n, ivj = 2 * n;
 
     Real      unfl = detail::Ops<Real>::lamch("safe minimum");
@@ -56,7 +56,7 @@ namespace arnoldi::detail {
 
         int ierr = -1;
         for (int itry = 1; itry <= 3; ++itry) {
-          getv0<Real, Backend>(bref, bmat, itry, false, n, j, v, ldv, resid, rnorm, workd, ierr, op, bop, comm);
+          getv0<Real, Backend>(bref, bmat, itry, false, n, j, v, ldv, resid, rnorm, workd, workc, ierr, op, bop, comm);
           if (ierr >= 0) break;
         }
         if (ierr < 0) {
@@ -148,10 +148,10 @@ namespace arnoldi::detail {
             detail::debug.vout(j, &h[(j - 1) * ldh], "_naitr: j-th column of H");
           }
 
-          detail::Ops<Real, Backend>::gemv(bref, "T", n, j, (Real)1, v, ldv, &workd[ipj], 1, (Real)0, &workd[irj], 1);
-          detail::backend_allreduce(bref, comm, &workd[irj], j);
-          detail::Ops<Real, Backend>::gemv(bref, "N", n, j, (Real)-1, v, ldv, &workd[irj], 1, (Real)1, resid, 1);
-          detail::Ops<Real, Backend>::axpy(bref, j, (Real)1, &workd[irj], 1, &h[(j - 1) * ldh], 1);
+          detail::Ops<Real, Backend>::gemv(bref, "T", n, j, (Real)1, v, ldv, &workd[ipj], 1, (Real)0, workc, 1);
+          detail::backend_allreduce(bref, comm, workc, j);
+          detail::Ops<Real, Backend>::gemv(bref, "N", n, j, (Real)-1, v, ldv, workc, 1, (Real)1, resid, 1);
+          detail::Ops<Real, Backend>::axpy(bref, j, (Real)1, workc, 1, &h[(j - 1) * ldh], 1);
 
           detail::arscnd(t2);
           Real rnorm1;
@@ -424,8 +424,8 @@ namespace arnoldi::detail {
   template <typename Real, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void naup2(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int& nev, int& np, Real tol,
              Real* resid, int mode, int iupd, int ishift, int& mxiter, Real* v, int ldv, Real* h, int ldh, Real* ritzr,
-             Real* ritzi, Real* bounds, Real* q, int ldq, Real* workl, Real* workd, int& info, OP&& op, BOP&& bop,
-             const Comm& comm) {
+             Real* ritzi, Real* bounds, Real* q, int ldq, Real* workl, Real* workd, Real* workc, int& info, OP&& op,
+             BOP&& bop, const Comm& comm) {
     double t0, t1, t2, t3;
     detail::arscnd(t0);
     int  msglvl = detail::debug.aup2;
@@ -450,14 +450,14 @@ namespace arnoldi::detail {
       detail::stats.aup2 = t1 - t0;
     };
 
-    getv0<Real, Backend>(bref, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, workd, info, op, bop, comm);
+    getv0<Real, Backend>(bref, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, workd, workc, info, op, bop, comm);
     if (rnorm == (Real)0) {
       info = -9;
       end_naup2();
       return;
     }
 
-    naitr<Real, Backend>(bref, bmat, n, 0, nev, mode, resid, rnorm, v, ldv, h, ldh, workd, info, op, bop, comm);
+    naitr<Real, Backend>(bref, bmat, n, 0, nev, mode, resid, rnorm, v, ldv, h, ldh, workd, workc, info, op, bop, comm);
     if (info > 0) {
       np     = info;
       mxiter = iter;
@@ -479,7 +479,7 @@ namespace arnoldi::detail {
         detail::debug.ivout(1, &np, "_naup2: Extend the Arnoldi factorization by");
       }
 
-      naitr<Real, Backend>(bref, bmat, n, nev, np, mode, resid, rnorm, v, ldv, h, ldh, workd, info, op, bop, comm);
+      naitr<Real, Backend>(bref, bmat, n, nev, np, mode, resid, rnorm, v, ldv, h, ldh, workd, workc, info, op, bop, comm);
       if (info > 0) {
         np     = info;
         mxiter = iter;
@@ -639,8 +639,8 @@ namespace arnoldi::detail {
   // naupd — nonsymmetric eigensolver entry point (callback version).
   template <typename Real, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void naupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid,
-             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op,
-             BOP&& bop, const Comm& comm) {
+             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl,
+             int& info, OP&& op, BOP&& bop, const Comm& comm) {
     detail::stats.reset();
     double t0, t1;
     detail::arscnd(t0);
@@ -706,7 +706,8 @@ namespace arnoldi::detail {
     ipntr[13]  = iw;
 
     naup2<Real, Backend>(bref, bmat, n, which, nev0, np, tol, resid, mode, iupd, ishift, mxiter, v, ldv, &workl[ih], ldh,
-                         &workl[ritzr], &workl[ritzi], &workl[bounds], &workl[iq], ldq, &workl[iw], workd, info, op, bop, comm);
+                         &workl[ritzr], &workl[ritzi], &workl[bounds], &workl[iq], ldq, &workl[iw], workd, workc,
+                         info, op, bop, comm);
 
     iparam[2]  = mxiter;
     iparam[4]  = np;
@@ -756,19 +757,20 @@ namespace arnoldi::detail {
   // op + bop, no comm (defaults to SerialComm).
   template <typename Real, typename Backend = detail::CpuBackend, typename OP, typename BOP>
   void naupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid,
-             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op,
-             BOP&& bop) {
-    naupd<Real, Backend>(bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info,
+             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl,
+             int& info, OP&& op, BOP&& bop) {
+    naupd<Real, Backend>(bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workc, workl, lworkl, info,
                          std::forward<OP>(op), std::forward<BOP>(bop), SerialComm{});
   }
 
   // Standard problem (bmat='I'), no bop, no comm.
   template <typename Real, typename Backend = detail::CpuBackend, typename OP>
   void naupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid,
-             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op) {
+             int ncv, Real* v, int ldv, int* iparam, int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl,
+             int& info, OP&& op) {
     naupd<Real, Backend>(
-        bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info, std::forward<OP>(op),
-        [](const Real*, Real*) {}, SerialComm{});
+        bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workc, workl, lworkl, info,
+        std::forward<OP>(op), [](const Real*, Real*) {}, SerialComm{});
   }
 
   // ---------------------------------------------------------------------------
@@ -776,21 +778,23 @@ namespace arnoldi::detail {
   // public callback API (tests, examples).
   template <typename Real, typename OP, typename BOP, typename Comm>
   void naupd(const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid, int ncv, Real* v, int ldv, int* iparam,
-             int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op, BOP&& bop, const Comm& comm) {
+             int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl, int& info, OP&& op, BOP&& bop,
+             const Comm& comm) {
     naupd<Real, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
-                                    ipntr, workd, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop), comm);
+                                    ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop),
+                                    comm);
   }
   template <typename Real, typename OP, typename BOP>
   void naupd(const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid, int ncv, Real* v, int ldv, int* iparam,
-             int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op, BOP&& bop) {
+             int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl, int& info, OP&& op, BOP&& bop) {
     naupd<Real, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
-                                    ipntr, workd, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop));
+                                    ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop));
   }
   template <typename Real, typename OP>
   void naupd(const char* bmat, int n, const char* which, int nev, Real& tol, Real* resid, int ncv, Real* v, int ldv, int* iparam,
-             int* ipntr, Real* workd, Real* workl, int lworkl, int& info, OP&& op) {
+             int* ipntr, Real* workd, Real* workc, Real* workl, int lworkl, int& info, OP&& op) {
     naupd<Real, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam,
-                                    ipntr, workd, workl, lworkl, info, std::forward<OP>(op));
+                                    ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op));
   }
 
   // neupd — nonsymmetric eigenvector extraction (callback version).
diff --git a/include/arnoldi/detail/sym.hpp b/include/arnoldi/detail/sym.hpp
index e3ae2f3..f047694 100644
--- a/include/arnoldi/detail/sym.hpp
+++ b/include/arnoldi/detail/sym.hpp
@@ -23,8 +23,8 @@ namespace arnoldi::detail {
   // Vectors (v, resid, workd) are Scalar; tridiagonal h is detail::real_t<Scalar>.
   template <typename Scalar, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void saitr(detail::BackendRef<Backend> bref, const char* bmat, int n, int k, int np, int mode, Scalar* resid,
-             detail::real_t<Scalar>& rnorm, Scalar* v, int ldv, detail::real_t<Scalar>* h, int ldh, Scalar* workd, int& info,
-             OP&& op, BOP&& bop, const Comm& comm) {
+             detail::real_t<Scalar>& rnorm, Scalar* v, int ldv, detail::real_t<Scalar>* h, int ldh, Scalar* workd,
+             Scalar* workc, int& info, OP&& op, BOP&& bop, const Comm& comm) {
     using Real    = detail::real_t<Scalar>;
     const int ipj = 0, irj = n, ivj = 2 * n;
 
@@ -44,6 +44,7 @@ namespace arnoldi::detail {
     std::vector<Scalar> hdiag_stage(static_cast<std::size_t>(k + np));
     std::vector<Scalar> hreorth_stage(static_cast<std::size_t>(k + np) * 2);
     std::vector<int>    hreorth_cnt(static_cast<std::size_t>(k + np), 0);
+
     auto                flush_hdiag = [&](int jlast) {
       detail::Ops<Scalar, Backend>::sync(bref);
       for (int jj = k + 1; jj <= jlast; ++jj) {
@@ -77,7 +78,7 @@ namespace arnoldi::detail {
 
         int ierr = -1;
         for (int itry = 1; itry <= 3; ++itry) {
-          getv0<Scalar, Backend>(bref, bmat, itry, false, n, j, v, ldv, resid, rnorm, workd, ierr, op, bop, comm);
+          getv0<Scalar, Backend>(bref, bmat, itry, false, n, j, v, ldv, resid, rnorm, workd, workc, ierr, op, bop, comm);
           if (ierr >= 0) break;
         }
         if (ierr < 0) {
@@ -139,15 +140,15 @@ namespace arnoldi::detail {
 
       if (mode != 2) {
         detail::Ops<Scalar, Backend>::gemv(bref, detail::Ops<Scalar, Backend>::herm_trans(), n, j, Scalar(1), v, ldv, &workd[ipj],
-                                           1, Scalar(0), &workd[irj], 1);
+                                           1, Scalar(0), workc, 1);
       } else {
         detail::Ops<Scalar, Backend>::gemv(bref, detail::Ops<Scalar, Backend>::herm_trans(), n, j, Scalar(1), v, ldv, &workd[ivj],
-                                           1, Scalar(0), &workd[irj], 1);
+                                           1, Scalar(0), workc, 1);
       }
-      detail::backend_allreduce(bref, comm, &workd[irj], j);
-      detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j, Scalar(-1), v, ldv, &workd[irj], 1, Scalar(1), resid, 1);
+      detail::backend_allreduce(bref, comm, workc, j);
+      detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j, Scalar(-1), v, ldv, workc, 1, Scalar(1), resid, 1);
 
-      detail::Ops<Scalar, Backend>::read_scalar_async(bref, &hdiag_stage[j - 1], &workd[irj + j - 1]);
+      detail::Ops<Scalar, Backend>::read_scalar_async(bref, &hdiag_stage[j - 1], &workc[j - 1]);
       if (j == 1 || rstart)
         h[j - 1] = Real(0);
       else
@@ -184,15 +185,15 @@ namespace arnoldi::detail {
           }
 
           detail::Ops<Scalar, Backend>::gemv(bref, detail::Ops<Scalar, Backend>::herm_trans(), n, j, Scalar(1), v, ldv,
-                                             &workd[ipj], 1, Scalar(0), &workd[irj], 1);
-          detail::backend_allreduce(bref, comm, &workd[irj], j);
-          detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j, Scalar(-1), v, ldv, &workd[irj], 1, Scalar(1), resid, 1);
+                                             &workd[ipj], 1, Scalar(0), workc, 1);
+          detail::backend_allreduce(bref, comm, workc, j);
+          detail::Ops<Scalar, Backend>::gemv(bref, "N", n, j, Scalar(-1), v, ldv, workc, 1, Scalar(1), resid, 1);
 
           if (j == 1 || rstart) h[j - 1] = Real(0);
           {
             int& rc = hreorth_cnt[j - 1];
             detail::Ops<Scalar, Backend>::read_scalar_async(bref, &hreorth_stage[static_cast<std::size_t>(j - 1) * 2 + rc],
-                                                            &workd[irj + j - 1]);
+                                                            &workc[j - 1]);
             ++rc;
           }
 
@@ -420,8 +421,8 @@ namespace arnoldi::detail {
   void saup2(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int& nev, int& np,
              detail::real_t<Scalar> tol, Scalar* resid, int mode, int iupd, int ishift, int& mxiter, Scalar* v, int ldv,
              detail::real_t<Scalar>* h, int ldh, detail::real_t<Scalar>* ritz, detail::real_t<Scalar>* bounds,
-             detail::real_t<Scalar>* q, int ldq, detail::real_t<Scalar>* workl, Scalar* workd, int& info, OP&& op, BOP&& bop,
-             const Comm& comm) {
+             detail::real_t<Scalar>* q, int ldq, detail::real_t<Scalar>* workl, Scalar* workd, Scalar* workc, int& info,
+             OP&& op, BOP&& bop, const Comm& comm) {
     using Real = detail::real_t<Scalar>;
     double t0, t1, t2, t3;
     detail::arscnd(t0);
@@ -447,14 +448,14 @@ namespace arnoldi::detail {
       detail::stats.aup2 = t1 - t0;
     };
 
-    getv0<Scalar, Backend>(bref, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, workd, info, op, bop, comm);
+    getv0<Scalar, Backend>(bref, bmat, 1, initv, n, 1, v, ldv, resid, rnorm, workd, workc, info, op, bop, comm);
     if (rnorm == Real(0)) {
       info = -9;
       end_saup2();
       return;
     }
 
-    saitr<Scalar, Backend>(bref, bmat, n, 0, nev0, mode, resid, rnorm, v, ldv, h, ldh, workd, info, op, bop, comm);
+    saitr<Scalar, Backend>(bref, bmat, n, 0, nev0, mode, resid, rnorm, v, ldv, h, ldh, workd, workc, info, op, bop, comm);
     if (info > 0) {
       np     = info;
       mxiter = iter;
@@ -473,7 +474,7 @@ namespace arnoldi::detail {
         detail::debug.ivout(1, &np, "_saup2: Extend the Lanczos factorization by");
       }
 
-      saitr<Scalar, Backend>(bref, bmat, n, nev, np, mode, resid, rnorm, v, ldv, h, ldh, workd, info, op, bop, comm);
+      saitr<Scalar, Backend>(bref, bmat, n, nev, np, mode, resid, rnorm, v, ldv, h, ldh, workd, workc, info, op, bop, comm);
       if (info > 0) {
         np     = info;
         mxiter = iter;
@@ -632,8 +633,8 @@ namespace arnoldi::detail {
   // workd is Scalar* (3n), workl is detail::real_t<Scalar>* (ncv^2 + 8*ncv).
   template <typename Scalar, typename Backend = detail::CpuBackend, typename OP, typename BOP, typename Comm>
   void saupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol,
-             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl,
-             int lworkl, int& info, OP&& op, BOP&& bop, const Comm& comm) {
+             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc,
+             detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op, BOP&& bop, const Comm& comm) {
     using Real = detail::real_t<Scalar>;
     detail::stats.reset();
     double t0, t1;
@@ -698,7 +699,8 @@ namespace arnoldi::detail {
     ipntr[10]  = iw + 1;
 
     saup2<Scalar, Backend>(bref, bmat, n, which, nev0, np, tol, resid, mode, iupd, ishift, mxiter, v, ldv, &workl[ih], ldh,
-                           &workl[ritz], &workl[bounds], &workl[iq], ldq, &workl[iw], workd, info, op, bop, comm);
+                           &workl[ritz], &workl[bounds], &workl[iq], ldq, &workl[iw], workd, workc, info, op, bop,
+                           comm);
 
     iparam[2]  = mxiter;
     iparam[4]  = np;
@@ -733,20 +735,20 @@ namespace arnoldi::detail {
   // op + bop, no comm (defaults to SerialComm).
   template <typename Scalar, typename Backend = detail::CpuBackend, typename OP, typename BOP>
   void saupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol,
-             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl,
-             int lworkl, int& info, OP&& op, BOP&& bop) {
-    saupd<Scalar, Backend>(bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info,
+             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc,
+             detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op, BOP&& bop) {
+    saupd<Scalar, Backend>(bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workc, workl, lworkl, info,
                            std::forward<OP>(op), std::forward<BOP>(bop), SerialComm{});
   }
 
   // Standard problem (bmat='I'), no bop, no comm.
   template <typename Scalar, typename Backend = detail::CpuBackend, typename OP>
   void saupd(detail::BackendRef<Backend> bref, const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol,
-             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl,
-             int lworkl, int& info, OP&& op) {
+             Scalar* resid, int ncv, Scalar* v, int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc,
+             detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op) {
     saupd<Scalar, Backend>(
-        bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workl, lworkl, info, std::forward<OP>(op),
-        [](const Scalar*, Scalar*) {}, SerialComm{});
+        bref, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, ipntr, workd, workc, workl, lworkl, info,
+        std::forward<OP>(op), [](const Scalar*, Scalar*) {}, SerialComm{});
   }
 
   // ---------------------------------------------------------------------------
@@ -754,24 +756,26 @@ namespace arnoldi::detail {
   // public callback API (tests, examples).
   template <typename Scalar, typename OP, typename BOP, typename Comm>
   void saupd(const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol, Scalar* resid, int ncv, Scalar* v,
-             int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op,
-             BOP&& bop, const Comm& comm) {
+             int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc, detail::real_t<Scalar>* workl, int lworkl,
+             int& info, OP&& op, BOP&& bop, const Comm& comm) {
     saupd<Scalar, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv,
-                                      iparam, ipntr, workd, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop),
-                                      comm);
+                                      iparam, ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op),
+                                      std::forward<BOP>(bop), comm);
   }
   template <typename Scalar, typename OP, typename BOP>
   void saupd(const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol, Scalar* resid, int ncv, Scalar* v,
-             int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op,
-             BOP&& bop) {
+             int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc, detail::real_t<Scalar>* workl, int lworkl,
+             int& info, OP&& op, BOP&& bop) {
     saupd<Scalar, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv,
-                                      iparam, ipntr, workd, workl, lworkl, info, std::forward<OP>(op), std::forward<BOP>(bop));
+                                      iparam, ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op),
+                                      std::forward<BOP>(bop));
   }
   template <typename Scalar, typename OP>
   void saupd(const char* bmat, int n, const char* which, int nev, detail::real_t<Scalar>& tol, Scalar* resid, int ncv, Scalar* v,
-             int ldv, int* iparam, int* ipntr, Scalar* workd, detail::real_t<Scalar>* workl, int lworkl, int& info, OP&& op) {
+             int ldv, int* iparam, int* ipntr, Scalar* workd, Scalar* workc, detail::real_t<Scalar>* workl, int lworkl,
+             int& info, OP&& op) {
     saupd<Scalar, detail::CpuBackend>(detail::BackendRef<detail::CpuBackend>{}, bmat, n, which, nev, tol, resid, ncv, v, ldv,
-                                      iparam, ipntr, workd, workl, lworkl, info, std::forward<OP>(op));
+                                      iparam, ipntr, workd, workc, workl, lworkl, info, std::forward<OP>(op));
   }
 
   // seupd — symmetric/Hermitian eigenvector extraction (callback).
diff --git a/tests/test_getv0.cpp b/tests/test_getv0.cpp
index 8186475..1c3d9ff 100644
--- a/tests/test_getv0.cpp
+++ b/tests/test_getv0.cpp
@@ -25,13 +25,13 @@ static void bop_diag(const double* x, double* y) {
 
 static bool test_basic() {
     const int n = 8;
-    std::vector<double> resid(n), workd(2 * n), v(n);
+    std::vector<double> resid(n), workd(2 * n), workc(n), v(n);
     double rnorm = 0;
     int ierr = 0;
     arnoldi::SerialComm comm;
 
     arnoldi::detail::getv0<double>("I", 1, false, n, 1,
-        v.data(), n, resid.data(), rnorm, workd.data(), ierr,
+        v.data(), n, resid.data(), rnorm, workd.data(), workc.data(), ierr,
         op_identity, [](const double*, double*){}, comm);
 
     assert(ierr == 0);
@@ -43,14 +43,14 @@ static bool test_basic() {
 // itry > 1 with bmat='G' takes the else-if branch in getv0.
 static bool test_generalized_retry() {
     const int n = 8;
-    std::vector<double> resid(n), workd(2 * n), v(n);
+    std::vector<double> resid(n), workd(2 * n), workc(n), v(n);
     double rnorm = 0;
     int ierr = 0;
     arnoldi::SerialComm comm;
 
     for (int i = 0; i < n; ++i) resid[i] = 1.0;
     arnoldi::detail::getv0<double>("G", 2, true, n, 1,
-        v.data(), n, resid.data(), rnorm, workd.data(), ierr,
+        v.data(), n, resid.data(), rnorm, workd.data(), workc.data(), ierr,
         op_identity, bop_diag, comm);
 
     assert(ierr == 0);
@@ -64,7 +64,7 @@ static bool test_orthogonalization() {
     const int n = 8;
     const int j = 2;
     std::vector<double> v(n * j, 0.0);
-    std::vector<double> resid(n), workd(2 * n);
+    std::vector<double> resid(n), workd(2 * n), workc(n);
     double rnorm = 0;
     int ierr = 0;
     arnoldi::SerialComm comm;
@@ -73,7 +73,7 @@ static bool test_orthogonalization() {
     for (int i = 0; i < n; ++i) resid[i] = 1.0;
 
     arnoldi::detail::getv0<double>("I", 1, true, n, j,
-        v.data(), n, resid.data(), rnorm, workd.data(), ierr,
+        v.data(), n, resid.data(), rnorm, workd.data(), workc.data(), ierr,
         op_identity, [](const double*, double*){}, comm);
 
     assert(ierr == 0);
@@ -90,7 +90,7 @@ static bool test_orthogonalization_generalized() {
     const int n = 8;
     const int j = 2;
     std::vector<double> v(n * j, 0.0);
-    std::vector<double> resid(n), workd(2 * n);
+    std::vector<double> resid(n), workd(2 * n), workc(n);
     double rnorm = 0;
     int ierr = 0;
     arnoldi::SerialComm comm;
@@ -99,7 +99,7 @@ static bool test_orthogonalization_generalized() {
     for (int i = 0; i < n; ++i) resid[i] = 1.0;
 
     arnoldi::detail::getv0<double>("G", 1, true, n, j,
-        v.data(), n, resid.data(), rnorm, workd.data(), ierr,
+        v.data(), n, resid.data(), rnorm, workd.data(), workc.data(), ierr,
         op_identity, bop_diag, comm);
 
     assert(ierr == 0);
@@ -111,7 +111,7 @@ static bool test_orthogonalization_generalized() {
 // msglvl > 3 triggers both the rnorm and full-vector debug prints.
 static bool test_debug_output() {
     const int n = 8;
-    std::vector<double> resid(n), workd(2 * n), v(n);
+    std::vector<double> resid(n), workd(2 * n), workc(n), v(n);
     double rnorm = 0;
     int ierr = 0;
     arnoldi::SerialComm comm;
@@ -120,7 +120,7 @@ static bool test_debug_output() {
     arnoldi::detail::debug.getv0 = 4;
 
     arnoldi::detail::getv0<double>("I", 1, false, n, 1,
-        v.data(), n, resid.data(), rnorm, workd.data(), ierr,
+        v.data(), n, resid.data(), rnorm, workd.data(), workc.data(), ierr,
         op_identity, [](const double*, double*){}, comm);
 
     arnoldi::detail::debug.getv0 = saved;
diff --git a/tests/test_low_level_validation.cpp b/tests/test_low_level_validation.cpp
index c7ca4ce..aaecee8 100644
--- a/tests/test_low_level_validation.cpp
+++ b/tests/test_low_level_validation.cpp
@@ -30,7 +30,7 @@ struct SaupdHarness {
     int lworkl;
     double tol = 0.0;
     int info = 0;
-    std::vector<double> resid, v, workd, workl;
+    std::vector<double> resid, v, workd, workc, workl;
     std::array<int, 11> iparam{};
     std::array<int, 14> ipntr{};
     arnoldi::SerialComm comm;
@@ -40,6 +40,7 @@ struct SaupdHarness {
         resid.assign(n, 0.0);
         v.assign(n * ncv, 0.0);
         workd.assign(3 * n, 0.0);
+        workc.assign(ncv, 0.0);
         workl.assign(lworkl, 0.0);
         iparam[0] = 1;     // ishift
         iparam[2] = 100;   // mxiter
@@ -50,7 +51,7 @@ struct SaupdHarness {
         info = 0;
         arnoldi::detail::saupd<double>(
             "I", n, "LM", nev, tol, resid.data(), ncv, v.data(), n,
-            iparam.data(), ipntr.data(), workd.data(), workl.data(), lworkl, info,
+            iparam.data(), ipntr.data(), workd.data(), workc.data(), workl.data(), lworkl, info,
             [](const double*, double*) {}, [](const double*, double*) {}, comm);
         return info;
     }
@@ -61,7 +62,7 @@ struct NaupdHarness {
     int lworkl;
     double tol = 0.0;
     int info = 0;
-    std::vector<double> resid, v, workd, workl;
+    std::vector<double> resid, v, workd, workc, workl;
     std::array<int, 11> iparam{};
     std::array<int, 14> ipntr{};
     arnoldi::SerialComm comm;
@@ -71,6 +72,7 @@ struct NaupdHarness {
         resid.assign(n, 0.0);
         v.assign(n * ncv, 0.0);
         workd.assign(3 * n, 0.0);
+        workc.assign(ncv, 0.0);
         workl.assign(lworkl, 0.0);
         iparam[0] = 1;
         iparam[2] = 100;
@@ -81,7 +83,7 @@ struct NaupdHarness {
         info = 0;
         arnoldi::detail::naupd<double>(
             "I", n, "LM", nev, tol, resid.data(), ncv, v.data(), n,
-            iparam.data(), ipntr.data(), workd.data(), workl.data(), lworkl, info,
+            iparam.data(), ipntr.data(), workd.data(), workc.data(), workl.data(), lworkl, info,
             [](const double*, double*) {}, [](const double*, double*) {}, comm);
         return info;
     }
@@ -109,7 +111,7 @@ TEST_CASE("test_saupd_errors", "[lowlevel]") {
         SaupdHarness h; arnoldi::SerialComm c;
         arnoldi::detail::saupd<double>(
             "I", h.n, "ZZ", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("saupd which invalid -> -5", h.info == -5, h.info, -5);
     }
@@ -117,7 +119,7 @@ TEST_CASE("test_saupd_errors", "[lowlevel]") {
         SaupdHarness h; arnoldi::SerialComm c;
         arnoldi::detail::saupd<double>(
             "X", h.n, "LM", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("saupd bmat invalid -> -6", h.info == -6, h.info, -6);
     }
@@ -135,7 +137,7 @@ TEST_CASE("test_saupd_errors", "[lowlevel]") {
         h.iparam[6] = 1;
         arnoldi::detail::saupd<double>(
             "G", h.n, "LM", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("saupd mode=1 + bmat=G -> -11", h.info == -11, h.info, -11);
     }
@@ -167,7 +169,7 @@ TEST_CASE("test_naupd_errors", "[lowlevel]") {
         NaupdHarness h; arnoldi::SerialComm c;
         arnoldi::detail::naupd<double>(
             "I", h.n, "ZZ", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("naupd which invalid -> -5", h.info == -5, h.info, -5);
     }
@@ -175,7 +177,7 @@ TEST_CASE("test_naupd_errors", "[lowlevel]") {
         NaupdHarness h; arnoldi::SerialComm c;
         arnoldi::detail::naupd<double>(
             "X", h.n, "LM", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("naupd bmat invalid -> -6", h.info == -6, h.info, -6);
     }
@@ -192,7 +194,7 @@ TEST_CASE("test_naupd_errors", "[lowlevel]") {
         h.iparam[6] = 1;
         arnoldi::detail::naupd<double>(
             "G", h.n, "LM", h.nev, h.tol, h.resid.data(), h.ncv, h.v.data(), h.n,
-            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workl.data(), h.lworkl, h.info,
+            h.iparam.data(), h.ipntr.data(), h.workd.data(), h.workc.data(), h.workl.data(), h.lworkl, h.info,
             [](const double*, double*) {}, [](const double*, double*) {}, c);
         check("naupd mode=1 + bmat=G -> -11", h.info == -11, h.info, -11);
     }
@@ -391,7 +393,7 @@ TEST_CASE("test_seupd_aliased_z_v", "[lowlevel]") {
     int   lworkl = ncv * (ncv + 8);
     double tol   = 0.0;
     int   info   = 0;
-    std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workl(lworkl, 0.0);
+    std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workc(ncv, 0.0), workl(lworkl, 0.0);
     std::array<int, 11> iparam{};
     std::array<int, 14> ipntr{};
     iparam[0] = 1; iparam[2] = 200; iparam[6] = 1;
@@ -405,7 +407,7 @@ TEST_CASE("test_seupd_aliased_z_v", "[lowlevel]") {
     auto bop = [](const double*, double*) {};
 
     arnoldi::detail::saupd<double>("I", n, "LM", nev, tol, resid.data(), ncv, v.data(), n,
-                                   iparam.data(), ipntr.data(), workd.data(), workl.data(),
+                                   iparam.data(), ipntr.data(), workd.data(), workc.data(), workl.data(),
                                    lworkl, info, op, bop, comm);
     check("aliased seupd: saupd ok", info >= 0, info, 0);
     if (info < 0) return;
@@ -431,14 +433,14 @@ TEST_CASE("test_no_bop_overloads", "[lowlevel]") {
         int   lworkl = ncv * (ncv + 8);
         double tol   = 0.0;
         int   info   = 0;
-        std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workl(lworkl, 0.0);
+        std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workc(ncv, 0.0), workl(lworkl, 0.0);
         std::array<int, 11> iparam{};
         std::array<int, 14> ipntr{};
         iparam[0] = 1; iparam[2] = 200; iparam[6] = 1;
 
         arnoldi::detail::saupd<double>(
             "I", n, "LM", nev, tol, resid.data(), ncv, v.data(), n,
-            iparam.data(), ipntr.data(), workd.data(), workl.data(), lworkl, info,
+            iparam.data(), ipntr.data(), workd.data(), workc.data(), workl.data(), lworkl, info,
             [n](const double* x, double* y) {
                 y[0] = 2.0 * x[0] - x[1];
                 for (int i = 1; i < n - 1; ++i) y[i] = -x[i - 1] + 2.0 * x[i] - x[i + 1];
@@ -451,14 +453,14 @@ TEST_CASE("test_no_bop_overloads", "[lowlevel]") {
         int   lworkl = 3 * ncv * ncv + 6 * ncv;
         double tol   = 0.0;
         int   info   = 0;
-        std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workl(lworkl, 0.0);
+        std::vector<double> resid(n, 0.0), v(n * ncv, 0.0), workd(3 * n, 0.0), workc(ncv, 0.0), workl(lworkl, 0.0);
         std::array<int, 11> iparam{};
         std::array<int, 14> ipntr{};
         iparam[0] = 1; iparam[2] = 200; iparam[6] = 1;
 
         arnoldi::detail::naupd<double>(
             "I", n, "LM", nev, tol, resid.data(), ncv, v.data(), n,
-            iparam.data(), ipntr.data(), workd.data(), workl.data(), lworkl, info,
+            iparam.data(), ipntr.data(), workd.data(), workc.data(), workl.data(), lworkl, info,
             [n](const double* x, double* y) {
                 y[0] = 2.0 * x[0] - x[1];
                 for (int i = 1; i < n - 1; ++i) y[i] = -x[i - 1] + 2.0 * x[i] - x[i + 1];