[GPU][LCAO] H(R) CSR output keeps extra near-zero entries on GPU vs CPU

[GenZeal-Lin]

Describe the bug

While validating existing 03_NAO_multik integration tests on GPU, I found that two LCAO cases produce different H(R) CSR sparsity patterns on CPU and GPU:

• tests/03_NAO_multik/scf_out_hsr_spin4
• tests/03_NAO_multik/nscf_out_hsr_tr_rr

Both CPU and GPU calculations finish successfully.

In these comparisons, hrs1_nao.csr contains additional near-zero entries on GPU, while the other generated CSR outputs compared below are byte-identical between CPU and GPU in file-level cmp comparisons.

For scf_out_hsr_spin4, the R = (0, 0, 0) block has CPU NNZ = 42 and GPU NNZ = 50. The eight GPU-only H(R) entries are around 1.1e-10 to 1.5e-10.

For nscf_out_hsr_tr_rr, the same block has CPU NNZ = 22 and GPU NNZ = 24. The two GPU-only entries are about 4.83e-10.

No CPU-only entries were found. At the current output precision, no differences were found among H(R) entries retained by both CPU and GPU.

The issue is reproducible on:

• 12ecc795f Refactor DFT+U codes (Refactor DFT+U codes #7852)
• 0dca9dd78 (earlier validation baseline)

These two cases are currently commented out of the CUDA CI case lists in #7856.

Expected behavior

What is the expected behavior for H(R) CSR output when CPU/GPU values differ only near the sparse-output threshold?

In particular, should this be addressed in:

1. the upstream H(R) numerical path;
2. CSR sparse-threshold handling; or
3. sparse-matrix integration-test comparison logic, subject to the project's test-governance requirements?

Numerically equivalent CPU/GPU calculations should ideally not produce materially different CSR sparsity structures solely because of near-zero values around the output threshold, unless that difference is expected by design.

To Reproduce

1. Checkout a tested revision, for example:
   12ecc795f

2. Build a CPU version with:
   USE_CUDA=OFF

3. Build a GPU version with:
   USE_CUDA=ON

4. Run the following cases with 2 MPI processes:

   • tests/03_NAO_multik/scf_out_hsr_spin4
   • tests/03_NAO_multik/nscf_out_hsr_tr_rr

5. For the GPU run, add:
   device gpu
   to the INPUT file.

6. Compare the generated CSR files with file-level cmp.

Observed results:

Case 1: scf_out_hsr_spin4

• CPU exit = 0
• GPU exit = 0

CSR comparison:

• hrs1_nao.csr DIFFERENT
• rr.csr SAME
• sr_nao.csr SAME

For R = (0, 0, 0):

• CPU NNZ = 42
• GPU NNZ = 50

GPU-only entries:

• (0, 4) ≈ -1.1448e-10
• (0, 10) ≈ -1.4810e-10
• (1, 5) ≈ -1.1448e-10
• (1, 11) ≈ -1.4810e-10
• (4, 0) ≈ -1.1448e-10
• (5, 1) ≈ -1.1448e-10
• (10, 0) ≈ -1.4810e-10
• (11, 1) ≈ -1.4810e-10

Case 2: nscf_out_hsr_tr_rr

• CPU exit = 0
• GPU exit = 0

CSR comparison:

• hrs1_nao.csr DIFFERENT
• rr.csr SAME
• sr_nao.csr SAME
• trs1_nao.csr SAME

For R = (0, 0, 0):

• CPU NNZ = 22
• GPU NNZ = 24

GPU-only entries:

• (2, 3) ≈ -4.82706319e-10
• (3, 2) ≈ -4.82706763e-10

Environment

• OS: Ubuntu 24.04
• GPU: NVIDIA GeForce RTX 4090
• MPI processes: 2
• Basis: LCAO
• ABACUS CPU binary: v3.11.0-beta8
• ABACUS GPU binary: v3.11.0-beta8
• Latest reproduced commit: 12ecc79
• Earlier reproduced commit: 0dca9dd
• GCC: gcc (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0
• NVCC / CUDA toolkit: Cuda compilation tools, release 12.6, V12.6.85
• CMake: 3.28.3
• CPU build: USE_CUDA=OFF
• GPU build: USE_CUDA=ON
• ENABLE_MPI=ON
• ENABLE_LCAO=ON
• ENABLE_ELPA=ON
• ENABLE_LIBXC=OFF
• ENABLE_DFTD4=OFF

Additional Context

The text H(R) CSR writer in source/source_io/module_hs/write_hs_r.cpp uses a fixed sparse threshold of 1e-10:

const double sparse_threshold = 1e-10;
hamilt::Output_HContainer out(
mat_serial,
ofs,
sparse_threshold,
precision
);

Output_HContainer in source/source_hamilt/module_hcontainer/output_hcontainer.cpp applies the same threshold when converting the already-computed HContainer values into CSR form:

sparse_matrix.setSparseThreshold(this->_sparse_threshold);

The matrix values are obtained from the existing HContainer and then inserted into the sparse matrix. The native binary writer in write_hs_r.cpp also uses a fixed 1e-10 threshold.

Therefore, CPU and GPU use the same CSR output implementation and the same sparse threshold.

The observed GPU-only H(R) values are close to the fixed 1e-10 sparse-output threshold. This suggests that small CPU/GPU numerical differences near zero may be interacting with the threshold and changing the resulting CSR sparsity pattern.

The upstream source of the GPU near-zero residual values has not yet been identified.

Related:

• Refactor(io): modularize sparse matrix writer interface #7468 Refactor(io): modularize sparse matrix writer interface
• Fix module_hs sparse output controls #7510 Fix module_hs sparse output controls
• test: add validated NAO cases to CUDA CI #7856 test: add validated NAO cases to CUDA CI

Task list for Issue attackers (only for developers)

• Verify the issue is not a duplicate.
• Describe the bug.
• Steps to reproduce.
• Expected behavior.
• Error message.
• Environment details.
• Additional context.
• Assign a priority level (low, medium, high, urgent).
• Assign the issue to a team member.
• Label the issue with relevant tags.
• Identify possible related issues.
• Create a unit test or automated test to reproduce the bug (if applicable).
• Fix the bug.
• Test the fix.
• Update documentation (if necessary).
• Close the issue and inform the reporter (if applicable).