Invalid memory acceses when executing ParICT preconditioner

[codecircuit]

I try to reproduce the results published in “PARILUT—A NEW PARALLEL THRESHOLD ILU FACTORIZATION” 2018. There, the ParILUT/ParICT preconditioners are proposed. I want to use the ParICT preconditioner with a PCG as an outer solver and matrix thermal1, which is available in the Sparse Matrix Collection. I checked that I read in the matrix correctly and the solver execution with the Jacobi preconditioner works fine. If the ParICT preconditioner is used cuda-memcheck reports invalid memory accesses, valgrind reports that conditional jumps depend on reads of uninitialized variables and the preconditioner terminates in a segmentation fault. Below you find my snippet which I used. I compiled with CUDA 11.0.2, GCC 9.2.0 and MAGMA with commit hash 23f549e9a238ca51710eadae1ff6a9d840f09aed, which is up to this moment the newest commit on the master branch. Maybe you can provide me with the code, which generates the ANI5, ANI6, or ANI7 matrices, which are mentioned in the paper? Then I can check if the errors still occur and provide a stand-alone snippet.

‌

‌

int main() {
Eigen::SparseMatrix<double> mtx = matrixmarket::legacy::readEigenSparseFromFile<double>(
    "thermal1.mtx");

auto csr = eigen2csr(mtx);

int i, m = mtx.cols(), n = 1;
double *rhs, *sol;

cudaMallocHost(&rhs, m * sizeof(double));
cudaMallocHost(&sol, m * sizeof(double));

std::default_random_engine rng(0);
std::uniform_real_distribution<double> dist(-0.5, 0.5);

for (int i = 0; i < m; ++i) {
    rhs[i] = dist(rng);
    sol[i] = 0.0;
}

magma_init();
magma_dopts opts;
magma_queue_t queue;
magma_queue_create(0, &queue);

memset(&opts, 0, sizeof(magma_dopts));

magma_d_matrix A = {Magma_CSR}, dA = {Magma_CSR};
magma_d_matrix b = {Magma_CSR}, db = {Magma_CSR};
magma_d_matrix x = {Magma_CSR}, dx = {Magma_CSR};

magma_dcsrset(m,
              m,
              csr.row_ptrs.data(),
              csr.col_indices.data(),
              csr.values.data(),
              &A,
              queue);
magma_dvset(m, 1, rhs, &b, queue);

opts.solver_par.solver = Magma_PCG;
opts.solver_par.restart = 0;
opts.solver_par.maxiter = 1000;
opts.solver_par.rtol = 1e-10;
opts.precond_par.solver = Magma_PARICT;
opts.precond_par.sweeps = 1;
opts.precond_par.levels = 0;
opts.precond_par.atol = 1.0;
opts.precond_par.trisolver = Magma_CUSOLVE;

// Initialize the solver
magma_int_t status_solver =
    magma_dsolverinfo_init(&opts.solver_par, &opts.precond_par, queue);
if (status_solver != MAGMA_SUCCESS) {
    std::cout << "Magma error while settung up the solver!" << std::endl;
    std::terminate();
}

// Copy the system to the device
magma_int_t status_dmtransfer0 =
    magma_dmtransfer(A, &dA, Magma_CPU, Magma_DEV, queue);
if (status_dmtransfer0 != MAGMA_SUCCESS) {
    std::cout << "Magma error while dmtransfer0" << std::endl;
    std::terminate();
}
magma_int_t status_dmtransfer1 =
    magma_dmtransfer(b, &db, Magma_CPU, Magma_DEV, queue);
if (status_dmtransfer1 != MAGMA_SUCCESS) {
    std::cout << "Magma error while dmtransfer1" << std::endl;
    std::terminate();
}
// initialize an initial guess for the iteration vector
magma_int_t status_dvinit =
    magma_dvinit(&dx, Magma_DEV, b.num_rows, b.num_cols, 0.0, queue);
if (status_dvinit != MAGMA_SUCCESS) {
    std::cout << "Magma error while dvinit" << std::endl;
    std::terminate();
}

// Generate the preconditioner
magma_int_t status = magma_d_precondsetup(
    dA, db, &opts.solver_par, &opts.precond_par, queue);
if (status != MAGMA_SUCCESS) {
    printf("Magma error while setting up the preconditioner: status = %d\n",
           (int) status);
    exit(EXIT_FAILURE);
}

magma_d_solver(dA, db, &dx, &opts, queue);

magma_dmtransfer(dx, &x, Magma_DEV, Magma_CPU, queue);
magma_dvcopy(x, &m, &n, sol, queue);

cudaDeviceSynchronize();

std::cout << "magma_numiter = " << opts.solver_par.numiter << std::endl;
std::cout << "init_res = " << opts.solver_par.init_res << std::endl;
std::cout << "final_res = " << opts.solver_par.final_res << std::endl;
std::cout << "relative_res = "
          << opts.solver_par.final_res / opts.solver_par.init_res
          << std::endl;

magma_dmfree(&dx, queue);
magma_dmfree(&db, queue);
magma_dmfree(&dA, queue);
magma_dmfree(&b,
             queue); // won't do anything as MAGMA does not own the data.
magma_dmfree(&A,
             queue); // won't do anything as MAGMA does not own the data.
magma_queue_destroy(queue);
magma_finalize();

for (i = 0; i < 20; ++i) {
    printf("%.4f\n", sol[i]);
}
cudaFreeHost(sol);
cudaFreeHost(rhs);
return 0;
}

‌