[Issue] Solving the same system on multi GPUs repeatedly gives different results

For some linear systems, repeatedly solving the same system on two or more GPUs with classical AMG as the preconditioner may not produce consistent results. In the reproducible example illustrated below, the first solve on two GPUs reaches the iteration limit without converging, whereas every subsequent solve of the exact same system converges in 13 iterations. The specific behavior may vary across machines.

Matrix Data

G3_circuit from SuiteSparse (https://sparse.tamu.edu/AMD/G3_circuit)

Reproduction steps

The following program (amgx_bug.c) reproduces the issue:

#include <stdio.h>
#include <mpi.h>
#include <cuda_runtime.h>
#include "amgx_c.h"

int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);
    MPI_Comm comm = MPI_COMM_WORLD;
    int rank, nranks, dev;
    MPI_Comm_rank(comm, &rank);
    MPI_Comm_size(comm, &nranks);
    cudaGetDeviceCount(&dev);
    dev = rank % dev;
    cudaSetDevice(dev);

    AMGX_initialize();
    AMGX_config_handle cfg;
    AMGX_config_create(&cfg,
        "config_version=2, solver(s)=PBICGSTAB, s:tolerance=1e-6, "
        "s:convergence=RELATIVE_INI, s:monitor_residual=1, "
        "s:preconditioner(p)=AMG, p:algorithm=CLASSICAL, p:interpolator=D2, "
        "p:strength_threshold=0.5, p:max_iters=1");
    AMGX_resources_handle rsrc;
    AMGX_resources_create(&rsrc, cfg, &comm, 1, &dev);
    AMGX_matrix_handle A;
    AMGX_vector_handle x, b;
    AMGX_solver_handle solver;
    AMGX_matrix_create(&A, rsrc, AMGX_mode_dDDI);
    AMGX_vector_create(&x, rsrc, AMGX_mode_dDDI);
    AMGX_vector_create(&b, rsrc, AMGX_mode_dDDI);
    AMGX_solver_create(&solver, rsrc, AMGX_mode_dDDI, cfg);

    // read the matrix and split its rows evenly across the ranks (b = 1 by default)
    AMGX_read_system_distributed(A, b, x, argv[1], 2, nranks, NULL, 0, NULL);

    for (int k = 0; k < 8; k++) {   // solve the same system 8 times
        if (k == 0) AMGX_solver_setup(solver, A);
        else        AMGX_solver_resetup(solver, A);
        AMGX_solver_solve_with_0_initial_guess(solver, b, x);
        int iters;
        AMGX_solver_get_iterations_number(solver, &iters);
        // all 8 counts should be equal; with the bug they are not (or the run crashes)
        if (rank == 0) printf("solve %d: %d iterations\n", k, iters);
    }

    AMGX_solver_destroy(solver);
    AMGX_vector_destroy(x);
    AMGX_vector_destroy(b);
    AMGX_matrix_destroy(A);
    AMGX_resources_destroy(rsrc);
    AMGX_config_destroy(cfg);
    AMGX_finalize();
    MPI_Finalize();
    return 0;
}

To reproduce the issue:

wget https://suitesparse-collection-website.herokuapp.com/MM/AMD/G3_circuit.tar.gz
tar xf G3_circuit.tar.gz
export AMGX=/path/to/AMGX CUDA_HOME=/usr/local/cuda
mpicc -I$AMGX/include -I$CUDA_HOME/include -L$AMGX/build -L$CUDA_HOME/lib64 \
      amgx_bug.c -lamgxsh -lcudart
export LD_LIBRARY_PATH=$AMGX/build:$CUDA_HOME/lib64:$LD_LIBRARY_PATH
mpirun -n 2 ./a.out G3_circuit/G3_circuit.mtx

On the current main branch, a typical output is shown below, although the exact behavior may vary across machines:

solve 0: 100 iterations      <- hit the iteration limit, did not converge
solve 1: 13 iterations
solve 2: 13 iterations
...
solve 7: 13 iterations

With 4 GPUs, it is the second solve that fails: 17, 100, 17, 17, 17, 17, 17, 17.

Environment information:

• OS: Ubuntu 22.04.5
• CUDA runtime: CUDA 13.1
• MPI version: OpenMPI 5.0.8
• AMGX version or commit hash: 91a8413e (current main branch)
• NVIDIA driver: 580.159.03
• NVIDIA GPU: NVIDIA L40 GPUs

Start with the amgx_bug.c reproduction and run it with the provided MPI, CUDA, and G3_circuit commands. Compare the AMGX_solver_setup and AMGX_solver_resetup paths across repeated solves and inspect the iteration counts. Done means repeated solves produce consistent convergence behavior without reaching the limit or crashing.