[SYCL] IC/ILU NotImplemented in Ginkgo 1.10: discrepancy with Tsai 2023 paper + 34M-cell int32 overflow in ISAI

[heikogleu-dev]

Context

Pioneer CFD work: OpenFOAM 13 + OGL + Ginkgo 1.10 SYCL on
Intel Arc Pro B70 Pro (Battlemage Xe2), Ubuntu 26.04, oneAPI 2026.0.

Finding 1: SYCL IC/ILU NotImplemented

Tsai et al. 2023 (Wiley CCPE, doi:10.1002/cpe.7400) explicitly states:

"We have ported the following preconditioner functionality to the DPC++
backend: ... ParILU, ParIC, ParILUT, ParICT, ISAI"

Ginkgo 1.10 release notes state:

"New implementations of the ILU and IC factorization for
CUDA, HIP, OpenMP, and Reference backends."
→ SYCL is absent from this list.

Our measurement confirms:
preconditioner IC with executor sycl → gko::NotImplemented
preconditioner ICT with executor sycl → DEVICE_LOST

Question: Is SYCL IC/ILU in a separate development branch?
What is the roadmap for getting it into a stable release?

Finding 2: ISAI int32 overflow at 34M cells

With ISAI sparsityPower=3, scaling=-1.0 on a 34M-cell CFD mesh
(np=8, ~4.25M cells/rank):
SYCL error: Provided range and/or offset does not fit in int

This suggests the ISAI sparsity pattern generation uses int32 internally,
which overflows at this problem size. ISAI sparsityPower=1 runs but diverges.

Question: Is there a 64-bit index path for ISAI on SYCL?

Finding 3: BJ maxBlockSize > 1 OOM on SYCL

Tested BS=2,4,8,16 — all crash with gko::AllocationError.
GINKGO_JACOBI_FULL_OPTIMIZATIONS=ON does not fix this.
The workspace allocation appears to be O(N × BS²) or worse in SYCL.

Profiling Results (for KIT/Intel interest)

Direct measurement on B70 Pro, 34M cells, np=8, BJ(1):

• GPU active: 34% of wall clock
• GPU idle: 66% (MPI AllReduce wait + forced host-buffer copies)
• When active: 2528 MHz avg (90% of 2800 spec)
• Kernel launch: 5.6 µs sync (CUDA-par — NOT 100 µs as some docs claim)
• forceHostBuffer=false: crashes in gko::mpi::all_to_all_v (no GPU-aware MPI for xe)

Full documentation:
https://github.com/heikogleu-dev/Openfoam13---GPU-Offloading-Intel-B70-Pro

Summary Question

With BJ(1) as only working SYCL preconditioner (always hits maxIter=200),
GPU CFD on Battlemage cannot compete with CPU GAMG (5-10 iter).
What is the realistic timeline for SYCL IC or Multigrid becoming
production-ready for distributed CFD workloads?

[nbeams]

Note the quote from the paper says

ParILU, ParIC, ParILUT, ParICT

(bolding mine). The Par* algorithms are not the same thing as ILU/IC. There is no discrepancy between the paper and the documentation.

There is a branch where @yhmtsai is working on getting direct factorizations working for SYCL, but we have run into issues on some Intel hardware (the A770, for example) that we haven't yet figured out.

How are you setting up your ISAI?

[heikogleu-dev]

Thank you for the clarification — you're absolutely right, I conflated
ParIC/ParILU with IC/ILU. We will correct this in our documentation.

That's also great news about @yhmtsai's branch — we'll keep watching.

To answer your question about ISAI setup, here is our exact OGL configuration
in fvSolution:
solvers
{
p
{
solver GKOCG;
executor sycl;
matrixFormat Csr;
preconditioner
{
preconditioner ISAI;
sparsityPower 1; // also tested 3 → int32 overflow
scaling -1.0; // for SPD pressure equation
}
forceHostBuffer true;
ranksPerGPU 8;
tolerance 1e-5;
relTol 0.01;
maxIter 200;
}
}

System: Distributed pressure equation from incompressible flow CFD

• 34M cells (4.25M per rank), np=8 ranks on single B70 Pro
• Symmetric Positive Definite (Laplacian, scaled by -1.0 per OGL docs)
• diagonally dominant

Observed: Initial residual = 1.0, Final residual = 8.68 after first solve.
Solver appears to diverge immediately rather than converge.

Question: Should we be testing ParIC instead of ISAI for this case?
We see in the documentation that ParIC is available for SYCL — we have
not tried it yet.

Also: the int32 overflow with ISAI sparsityPower=3 — is there a 64-bit
index path planned for SYCL ISAI sparsity pattern generation?

Full reproduction context:
https://github.com/heikogleu-dev/Openfoam13---GPU-Offloading-Intel-B70-Pro

[heikogleu-dev]

Quick follow-up — re-tested with the correct OGL→Ginkgo mapping
(OGL ILU = ParIlu, OGL ICT = ParIct). Both reach factorization
but fail later:

ILU (ParIlu factorization OK, apply fails):
gko::NotImplemented at:
dpcpp/solver/lower_trs_kernels.dp.cpp:43: generate is not implemented

ICT (ParIct factorization itself crashes):
SIGABRT in gko::kernels::dpcpp::par_ict_factorization::add_candidates<double, int>

Tested with and without scaling=-1.0 — same outcome.

So the real gap on Battlemage seems to be on the apply side
(lower_trs/upper_trs SYCL kernels missing) plus the
par_ict_factorization::add_candidates crash. ISAI question from my
previous comment still stands.

Updated documentation:
https://github.com/heikogleu-dev/Openfoam13---GPU-Offloading-Intel-B70-Pro/blob/main/findings/05_sycl_preconditioner_status.md

Happy to test @yhmtsai's branch on B70 Pro if that would be useful —
the Battlemage bug profile looks different from the A770 issues you
mentioned.

[nbeams]

The Isai class can have an IndexType of int64 as far as I know, but I'm not sure if that's currently an option within OGL. I will have to defer to @greole on that.

Do you have any more information about the ParICT crash? (Backtrace, etc.)

Since you have a diagonally dominant matrix, we can hope that ParIC should work pretty well. One way that the ISAI can be used is to approximate the triangular solves for applying the (Par)IC factors, rather than as a stand-alone preconditioner. You could try that to get around the missing lower triangular solve (with ISAI sparsity power of 1). But again, I'm not sure how to do that in OGL.

[heikogleu-dev]

Thank you very much — that's a really helpful suggestion. The
ISAI-as-trs-approximation approach is exactly the kind of workaround
we were hoping for. We'll investigate whether OGL exposes this
configuration or whether it would require an OGL patch (cc @greole).

Regarding the ParICT crash backtrace, here is the full stack from
our run with OGL preconditioner=ICT, scaling=-1.0, np=8 on B70 Pro
(BMG-G31, Battlemage Xe2), 34M-cell distributed CFD pressure system:
SIGABRT in:
gko::kernels::dpcpp::par_ict_factorization::add_candidates<double, int>(
std::shared_ptr,
matrix::Csr<double, int> const*,
matrix::Csr<double, int> const*,
matrix::Csr<double, int> const*,
matrix::Csr<double, int>*)
in libginkgo_dpcpp.so.1.10.0
Called from:
gko::factorization::ParIct<double, int>::generate_l_lt(...) const
→ ExecutorBase::run(Operation const&) const
→ EnableDefaultFactory<ParIct<double,int>::Factory, ...>::generate_impl(...)
→ LinOpFactory::generate(...)
All 8 MPI ranks abort identically. No prior gko::NotImplemented
exception — direct SIGABRT inside the kernel.

Stack: Compute Runtime 26.05.37020.3-1, IGC 2.32.7, oneAPI 2026.0,
OGL master built against Ginkgo 1.10 with GINKGO_JACOBI_FULL_OPTIMIZATIONS=ON
and -fp-model=precise. Linux 7.0, xe driver, Ubuntu 26.04.

Full crash log preserved at the repo if useful for debugging:
https://github.com/heikogleu-dev/Openfoam13---GPU-Offloading-Intel-B70-Pro/blob/main/findings/05_sycl_preconditioner_status.md

The crash occurs deterministically in add_candidates — no other
output between the factorization start log and the SIGABRT. This
looks like a real SYCL kernel issue rather than a stack/setup problem,
since ParIc factorization (via OGL ILU) succeeds on the same setup
— only the apply side fails there with lower_trs_kernels.dp.cpp:43: generate is not implemented.

We'll watch for @greole's input on whether ISAI-as-trs is reachable
through OGL today, and we'd be glad to test @yhmtsai's branch on
Battlemage — the bug profile looks distinctly different from A770
(direct SIGABRT in add_candidates vs. whatever you've seen there).

[nbeams]

Full crash log preserved at the repo if useful for debugging

I didn't see the full crash log there. Can you point me directly to it?

[heikogleu-dev]

Apologies — the full raw log was on my system but not yet in the repo.
It's there now:

https://github.com/heikogleu-dev/Openfoam13---GPU-Offloading-Intel-B70-Pro/blob/main/logs/paricT_crash_B70Pro_CR26.05_IGC2.32.7.log

Stack: B70 Pro (BMG-G31), CR 26.05.37020.3-1, IGC 2.32.7,
oneAPI 2026.0, xe driver, Ubuntu 26.04 LTS, Kernel 7.0.0-15.
np=8 ranks, 34M cells, OGL preconditioner=ICT, scaling=-1.0.

The SIGABRT occurs deterministically at the same point in all 8 ranks,
immediately after OGL's verbose "scale_RHS" log entry — before any
factorization output. Let me know if you need any additional context
or a different configuration.

[nbeams]

Just a note, #2023 has been merged, so now triangular solvers are available for dpcpp.