Skip to content

Introduce additional work array to avoid data overflow #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
iskakoff merged 1 commit into from
May 20, 2026
Merged

Introduce additional work array to avoid data overflow #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion examples/symmetric_laplacian.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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) {
Expand Down
20 changes: 13 additions & 7 deletions include/arnoldi/arnoldi.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_;
Expand Down Expand Up @@ -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_);
}
}

Expand All @@ -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_);
}
}

Expand Down Expand Up @@ -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_;
Expand Down
15 changes: 8 additions & 7 deletions include/arnoldi/detail/getv0.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down
56 changes: 30 additions & 26 deletions include/arnoldi/detail/nonsym.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -756,41 +757,44 @@ 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{});
}

// ---------------------------------------------------------------------------
// naupd no-BackendRef overloads — defaults to CpuBackend for the historical
// 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).
Expand Down
Loading
Loading