Seeking Advice on Configuring AMGX for Solving a Poorly-Conditioned Matrix Equation

[fxuecnu]

Hello,

I'm working on solving a matrix equation that seems to be ill-conditioned, and I'm encountering convergence issues with all the configuration files provided in AMGX. I've attached a JPEG image of the matrix to give an idea of its structure.

Despite confirming that the matrix rank is appropriate, I'm struggling with the equation's principal axis, which contains many zero values. This characteristic might be contributing to the convergence problems.

I'm seeking advice on how to properly configure AMGX for this type of equation. To provide more context, I've also attached a snippet of my code that illustrates how I'm constructing and attempting to solve this system.

Any suggestions or guidance on the appropriate AMGX configuration settings for this scenario would be greatly appreciated.

Thank you for your time and help.

//// below is code

/* Copyright (c) 2011-2017, NVIDIA CORPORATION. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *  * Neither the name of NVIDIA CORPORATION nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include "cuda_runtime.h"

#include <stdint.h>

#define M_PI 3.14159265358979323846

/* CUDA error macro */
#define CUDA_SAFE_CALL(call) do {                                 \
  cudaError_t err = call;                                         \
  if(cudaSuccess != err) {                                        \
    fprintf(stderr, "Cuda error in file '%s' in line %i : %s.\n", \
            __FILE__, __LINE__, cudaGetErrorString( err) );       \
    exit(EXIT_FAILURE);                                           \
  } } while (0)

//#define AMGX_DYNAMIC_LOADING
//#undef AMGX_DYNAMIC_LOADING
#define MAX_MSG_LEN 4096

/* standard or dynamically load library */
#ifdef AMGX_DYNAMIC_LOADING
#include "amgx_capi.h"
#else
#include "amgx_c.h"
#endif

/* print error message and exit */
void errAndExit(const char *err)
{
    printf("%s\n", err);
    fflush(stdout);
    MPI_Abort(MPI_COMM_WORLD, 1);
    exit(1);
}

/* print callback (could be customized) */
void print_callback(const char *msg, int length)
{
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    if (rank == 0) { printf("%s", msg); }
}

/* print usage and exit */
void printUsageAndExit()
{
    char msg[MAX_MSG_LEN] = "Usage: mpirun [-n nranks] ./example [-mode [dDDI | dDFI | dFFI]] [-p nx ny] [-c config_file] [-amg \"variable1=value1 ... variable3=value3\"] [-gpu] [-it k]\n";
    strcat(msg, "     -mode:                select the solver mode\n");
    strcat(msg, "     -p nx ny:             select x- and y-dimensions of the 2D (5-points) local discretization of the Poisson operator (the global problem size will be nranks*nx*ny)\n");
    strcat(msg, "     -c:                   set the amg solver options from the config file\n");
    strcat(msg, "     -amg:                 set the amg solver options from the command line\n");
    print_callback(msg, MAX_MSG_LEN);
    MPI_Finalize();
    exit(0);
}

/* parse parameters */
int findParamIndex(char **argv, int argc, const char *parm)
{
    int count = 0;
    int index = -1;

    for (int i = 0; i < argc; i++)
    {
        if (strncmp(argv[i], parm, 100) == 0)
        {
            index = i;
            count++;
        }
    }

    if (count == 0 || count == 1)
    {
        return index;
    }
    else
    {
        char msg[MAX_MSG_LEN];
        sprintf(msg, "ERROR: parameter %s has been specified more than once, exiting\n", parm);
        print_callback(msg, MAX_MSG_LEN);
        exit(1);
    }

    return -1;
}

void continuityXdirKernel(
    int n,
    int nx,
    int ny,
    double dx,
    int* ndd,
    double* u,
    double* v,
    double* h,
    double* detadt,
    void* values,
    void* rhs,
    int* row_ptrs) {

    /* -------------------- continuity equation ---------------------------
    -----------------------detadt = dhu/dx + dhv/dy------------------------
    detadt is constant, so it is not included in the matrix. In other words,
    the pivot element is 0 for the continuity equation.
    -------------------------------------------------------------------- */

    for (int i = 0; i < n; i++) {

        // note that here n, nx, ny is different from those of momentum equation

        int row = i / nx;                               // row index in the 2D work array       
        int col = i % nx;                               // col index in the 1D work array    
        int idx = 2 * col + 1 + (2 * nx + 1) * row;     // pivot index in the row of sparse matrix
        int nnz = row_ptrs[idx];						// index of the first non-zero element in the row


        // open boundary condition for water level
        if (ndd[i] == 2) {

            ((double*)values)[nnz] = 1.0;
            ((double*)rhs)[idx] = 0.0;

        }

        // in other cases, just discrete the equation 
        else {

            // left u coefficient
            ((double*)values)[nnz] = -0.5 * (h[i] + h[i - 1]) / dx;

            // right u coefficient
            ((double*)values)[nnz + 1] = 0.5 * (h[i] + h[min(i + 1, nx - 1)]) / dx;

            //rhs values
            ((double*)rhs)[idx] = detadt[i];

        }

    }

}


void momentumXdirKernel(int n,
    int nx,
    int ny,
    double dx,
    double gra,
    double mann,
    double Du,
    int* ndd,
    double* u,
    double* v,
    double* eta,
    double* h,
    void* values,
    void* rhs,
    int* row_ptrs) {

    /* -------------------- momentum equation u ---------------------------
    * 0 = g * mann^2 / h^(4/3) * |uv| * u     // friction term
          + udu/dx + vdu/dy                   // advection term
          g * deta/dx                         // gravity term
          - d/dx (Du * du/dx)                 // diffusion term
    -------------------------------------------------------------------- */

    //int i = blockIdx.x * blockDim.x + threadIdx.x;
    for (int i = 0; i < n; i++) {

        // note that here n, nx is different from those of continuity equation

        int row = i / nx;                               // row index in the 2D work array       
        int col = i % nx;                               // col index in the 1D work array    
        int idx = 2 * col + (2 * nx - 1) * row;         // pivot index in the row of sparse matrix
        int nnz = row_ptrs[idx];						// index of the first non-zero element in the row


        // inner points
        if (col > 0 && col < nx - 1) {

            int i0 = col + row * (nx - 1);

            double h_u_c = 0.5 * (h[i0] + h[i0 - 1]);   // h_u_c = h_u_center
            double D_h_l = Du;
            double D_h_r = Du;

            // left U
            ((double*)values)[nnz] = -fmax(0., u[i]) / dx      // upwind advection term
                - D_h_l / dx / dx;         // diffusion term

            // left eta
            ((double*)values)[nnz + 1] = -gra / dx;               // gravity term

            // center U
            ((double*)values)[nnz + 2] = gra * pow(mann, 2) / pow(h_u_c, 4.0 / 3.0) * fabs(u[i])       // friction term
                + (fmax(0., u[i]) - fmin(0., u[i])) / dx                        // upwind advection term
                + (D_h_l + D_h_r) / dx / dx;                                    // diffusion term

            // right eta
            ((double*)values)[nnz + 3] = gra / dx;                 // gravity term

            // right U
            ((double*)values)[nnz + 4] = fmin(0., u[i]) / dx       // upwind advection term
                - D_h_r / dx / dx;          // diffusion term

            //rhs values
            ((double*)rhs)[idx] = 0.;

        }
        // left boundary
        else if (col == 0) {

            int i0 = col + row * (nx - 1);

            // open boundary condition
            if (ndd[i0] == 2) {

                // center U
                ((double*)values)[nnz] = 1.0;

                // right U
                ((double*)values)[nnz + 1] = -2.0;

                // right right U
                ((double*)values)[nnz + 2] = 1.0;

                //rhs values
                ((double*)rhs)[idx] = 0.;

            }

            // closed boundary condition
            else {

                ((double*)values)[nnz] = 1.0;
                ((double*)rhs)[idx] = 0.;

            }

        }
        // right boundary
        else if (col == nx - 1) {

            int i0 = col - 1 + row * (nx - 1);

            // open boundary condition
            if (ndd[i0] == 2) {

                // left left U
                ((double*)values)[nnz] = 1.0;

                // left U
                ((double*)values)[nnz + 1] = -2.0;

                // center U
                ((double*)values)[nnz + 2] = 1.0;

                //rhs values
                ((double*)rhs)[idx] = 0.;

            }

            // closed boundary condition
            else {

                ((double*)values)[nnz] = 1.0;
                ((double*)rhs)[idx] = 0.;

            }

        }
    }

}


int main(int argc, char **argv)
{
    //parameter parsing
    int pidx = 0;
    int pidy = 0;
    //MPI (with CUDA GPUs)
    int rank = 0;
    int lrank = 0;
    int nranks = 0;
    int n;
    int nx, ny;
    int gpu_count = 0;
    MPI_Comm amgx_mpi_comm = MPI_COMM_WORLD;
    //versions
    int major, minor;
    char *ver, *date, *time;
    //input matrix and rhs/solution
    int *partition_sizes = NULL;
    int *partition_vector = NULL;
    int partition_vector_size = 0;
    //library handles
    AMGX_Mode mode;
    AMGX_config_handle cfg;
    AMGX_resources_handle rsrc;
    AMGX_matrix_handle A;
    AMGX_vector_handle b, x;
    AMGX_solver_handle solver;
    //status handling
    AMGX_SOLVE_STATUS status;
    /* MPI init (with CUDA GPUs) */
    //MPI
    MPI_Init(&argc, &argv);
    MPI_Comm_size(amgx_mpi_comm, &nranks);
    MPI_Comm_rank(amgx_mpi_comm, &rank);
    //CUDA GPUs
    CUDA_SAFE_CALL(cudaGetDeviceCount(&gpu_count));
    lrank = rank % gpu_count;
    CUDA_SAFE_CALL(cudaSetDevice(lrank));
    printf("Process %d selecting device %d\n", rank, lrank);

    /* check arguments */
    if (argc == 1)
    {
        printUsageAndExit();
    }

    /* load the library (if it was dynamically loaded) */
#ifdef AMGX_DYNAMIC_LOADING
    void *lib_handle = NULL;
#ifdef _WIN32
    lib_handle = amgx_libopen("amgxsh.dll");
#else
    lib_handle = amgx_libopen("libamgxsh.so");
#endif

    if (lib_handle == NULL)
    {
        errAndExit("ERROR: can not load the library");
    }

    //load all the routines
    if (amgx_liblink_all(lib_handle) == 0)
    {
        amgx_libclose(lib_handle);
        errAndExit("ERROR: corrupted library loaded\n");
    }

#endif
    /* init */
    AMGX_SAFE_CALL(AMGX_initialize());
    AMGX_SAFE_CALL(AMGX_initialize_plugins());
    /* system */
    AMGX_SAFE_CALL(AMGX_register_print_callback(&print_callback));
    AMGX_SAFE_CALL(AMGX_install_signal_handler());

    /* get api and build info */
    if ((pidx = findParamIndex(argv, argc, "--version")) != -1)
    {
        AMGX_get_api_version(&major, &minor);
        printf("amgx api version: %d.%d\n", major, minor);
        AMGX_get_build_info_strings(&ver, &date, &time);
        printf("amgx build version: %s\nBuild date and time: %s %s\n", ver, date, time);
        AMGX_SAFE_CALL(AMGX_finalize_plugins());
        AMGX_SAFE_CALL(AMGX_finalize());
        /* close the library (if it was dynamically loaded) */
#ifdef AMGX_DYNAMIC_LOADING
        amgx_libclose(lib_handle);
#endif
        MPI_Finalize();
        exit(0);
    }

    /* get mode */
    if ((pidx = findParamIndex(argv, argc, "-mode")) != -1)
    {
        if (strncmp(argv[pidx + 1], "dDDI", 100) == 0)
        {
            mode = AMGX_mode_dDDI;
        }
        else if (strncmp(argv[pidx + 1], "dDFI", 100) == 0)
        {
            mode = AMGX_mode_dDFI;
        }
        else if (strncmp(argv[pidx + 1], "dFFI", 100) == 0)
        {
            mode = AMGX_mode_dFFI;
        }
        else
        {
            errAndExit("ERROR: invalid mode");
        }
    }
    else
    {
        printf("Warning: No mode specified, using dDDI by default.\n");
        mode = AMGX_mode_dDDI;
    }

    int sizeof_m_val = ((AMGX_GET_MODE_VAL(AMGX_MatPrecision, mode) == AMGX_matDouble)) ? sizeof(double) : sizeof(float);
    int sizeof_v_val = ((AMGX_GET_MODE_VAL(AMGX_VecPrecision, mode) == AMGX_vecDouble)) ? sizeof(double) : sizeof(float);
    /* create config */
    pidx = findParamIndex(argv, argc, "-amg");
    pidy = findParamIndex(argv, argc, "-c");

    if ((pidx != -1) && (pidy != -1))
    {
        printf("%s\n", argv[pidx + 1]);
        AMGX_SAFE_CALL(AMGX_config_create_from_file_and_string(&cfg, argv[pidy + 1], argv[pidx + 1]));
    }
    else if (pidy != -1)
    {
        AMGX_SAFE_CALL(AMGX_config_create_from_file(&cfg, argv[pidy + 1]));
    }
    else if (pidx != -1)
    {
        printf("%s\n", argv[pidx + 1]);
        AMGX_SAFE_CALL(AMGX_config_create(&cfg, argv[pidx + 1]));
    }
    else
    {
        errAndExit("ERROR: no config was specified");
    }

    /* example of how to handle errors */
    //char msg[MAX_MSG_LEN];
    //AMGX_RC err_code = AMGX_resources_create(NULL, cfg, &amgx_mpi_comm, 1, &lrank);
    //AMGX_SAFE_CALL(AMGX_get_error_string(err_code, msg, MAX_MSG_LEN));
    //printf("ERROR: %s\n",msg);
    /* switch on internal error handling (no need to use AMGX_SAFE_CALL after this point) */
    AMGX_SAFE_CALL(AMGX_config_add_parameters(&cfg, "exception_handling=1"));
    /* create resources, matrix, vector and solver */
    AMGX_resources_create(&rsrc, cfg, &amgx_mpi_comm, 1, &lrank);
    AMGX_matrix_create(&A, rsrc, mode);
    AMGX_vector_create(&x, rsrc, mode);
    AMGX_vector_create(&b, rsrc, mode);
    AMGX_solver_create(&solver, rsrc, mode, cfg);
    //generate 3D Poisson matrix, [and rhs & solution]
    //WARNING: use 1 ring for aggregation and 2 rings for classical path
    int nrings; //=1; //=2;
    AMGX_config_get_default_number_of_rings(cfg, &nrings);
    //printf("nrings=%d\n",nrings);
    int nglobal = 0;

    if  ((pidx = findParamIndex(argv, argc, "-p")) != -1)
    {
        nx = atoi(argv[++pidx]);
        ny = atoi(argv[++pidx]);
        n  = nx * ny; // each rank has the strip of size nx*ny
        nglobal = n * nranks; // global domain is just those strips stacked one on each other
    }
    else
    {
        printf("Please, use '-p nx ny' parameter for this example\n");
        exit(1);
    }

    nx = 5000;
    ny = 1;
    n = nx * ny;
    nglobal = n * nranks;

    /* generate the matrix
       In more detail, this routine will create 2D (5 point) discretization of the
       Poisson operator. The discretization is performed on a the 2D domain consisting
       of nx and ny points in x- and y-dimension respectively. Each rank processes it's
       own part of discretization points. Finally, the rhs and solution will be set to
       a vector of ones and zeros, respectively. */

    int *row_ptrs = (int*)malloc((2 * n + ny + 1) * sizeof(int));
    int64_t *col_idxs = (int64_t*)malloc(6 * (2 * n + ny) * sizeof(int64_t));

    void *values = malloc(6 * (2 * n + ny) * sizeof(double)); // maximum nnz
    int nnz = 0;
    int64_t count = 0;
    int64_t start_idx = rank * n;

    //for (int i = 0; i < n; i ++)
    //{
    //    row_ptrs[i] = nnz;

    //    if (rank > 0 || i > ny-1)
    //    {
    //        col_idxs[nnz] = (i + start_idx - ny);

    //        if (sizeof_m_val == 4)
    //        {
    //            ((float *)values)[nnz] = -1.f;
    //        }
    //        else if (sizeof_m_val == 8)
    //        {
    //            ((double *)values)[nnz] = -1.;
    //        }

    //        nnz++;
    //    }

    //    if (i % ny != 0)
    //    {
    //        col_idxs[nnz] = (i + start_idx - 1);

    //        if (sizeof_m_val == 4)
    //        {
    //            ((float *)values)[nnz] = -1.f;
    //        }
    //        else if (sizeof_m_val == 8)
    //        {
    //            ((double *)values)[nnz] = -1.;
    //        }

    //        nnz++;
    //    }

    //    {
    //        col_idxs[nnz] = (i + start_idx);

    //        if (sizeof_m_val == 4)
    //        {
    //            ((float *)values)[nnz] = 4.f;
    //        }
    //        else if (sizeof_m_val == 8)
    //        {
    //            ((double *)values)[nnz] = 4.;
    //        }

    //        nnz++;
    //    }

    //    if ((i + 1) % ny != 0)
    //    {
    //        col_idxs[nnz] = (i + start_idx + 1);

    //        if (sizeof_m_val == 4)
    //        {
    //            ((float *)values)[nnz] = -1.f;
    //        }
    //        else if (sizeof_m_val == 8)
    //        {
    //            ((double *)values)[nnz] = -1.;
    //        }

    //        nnz++;
    //    }

    //    if ( (rank != nranks - 1) || (i / ny != (nx - 1)) )
    //    {
    //        col_idxs[nnz] = (i + start_idx + ny);

    //        if (sizeof_m_val == 4)
    //        {
    //            ((float *)values)[nnz] = -1.f;
    //        }
    //        else if (sizeof_m_val == 8)
    //        {
    //            ((double *)values)[nnz] = -1.;
    //        }

    //        nnz++;
    //    }
    //}

    //row_ptrs[n] = nnz;

    ///////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////

    double dx = .5;
    double gra = 9.81;
    double mann = 0.012;
    double Du = 0.5;

    double Trange = 1.0;
    double Tperiod = 43200.0;

    int* ndd;
    double* u;
    double* v;
    double* eta;
    double* h;
    double* detadt;
    void* rhs;

    // Allocate memory on host
    ndd = (int*)malloc(n * sizeof(int));
    u = (double*)malloc((n + ny) * sizeof(double));
    v = (double*)malloc((n + nx) * sizeof(double));
    eta = (double*)malloc(n * sizeof(double));
    h = (double*)malloc(n * sizeof(double));
    detadt = (double*)malloc(n * sizeof(double));
    rhs = malloc((2 * n + ny) * sizeof(double));


    // Initialize host memory
    for (int i = 0; i < n; i++) {

        int row = i / nx;                               // row index in the 2D work array   
        int col = i % nx;                               // col index in the 1D work array

        if (col == 0) {
            ndd[i] = 2;
        }
        else {
            ndd[i] = 0;
        }

        eta[i] = 0.0;
        h[i] = 2.0;
        detadt[i] = M_PI * Trange / Tperiod;

    }

    for (int i = 0; i < n + ny; i++) {
        u[i] = 0.0;
    }

    for (int i = 0; i < n + nx; i++) {
        v[i] = 0.0;
    }

    // initialize sparse matrix (determine row_ptrs and col_idxs)


    nnz = 0;

    for (int row = 0; row < ny; row++) {

        for (int i = 0; i < 2 * nx + 1; i++) {

            // continuity equation
            if (i % 2 == 1) {

                int col = (i - 1) / 2;
                int idx = i + (2 * nx + 1) * row;
                int i0 = col + row * nx;
                row_ptrs[idx] = nnz;

                if (ndd[i0] == 2) {

                    col_idxs[nnz++] = idx;

                }
                else {

                    col_idxs[nnz++] = idx - 1;
                    col_idxs[nnz++] = idx + 1;

                }

            }

            // momentum equation
            if (i % 2 == 0) {

                int col = i / 2;
                int idx = i + (2 * nx + 1) * row;
                int i0 = col + row * nx;
                row_ptrs[idx] = nnz;

                if (col > 0 && col < nx) {

                    col_idxs[nnz++] = idx - 2;
                    col_idxs[nnz++] = idx - 1;
                    col_idxs[nnz++] = idx;
                    col_idxs[nnz++] = idx + 1;
                    col_idxs[nnz++] = idx + 2;

                }
                else if (col == 0) {

                    if (ndd[i0] == 2) {

                        printf("idx = %d, nnz = %d\n", idx, nnz);

                        col_idxs[nnz++] = idx;
                        col_idxs[nnz++] = idx + 2;
                        col_idxs[nnz++] = idx + 4;

                    }
                    else {

                        col_idxs[nnz++] = idx;

                    }
                }
                else if (col == nx) {

                    int i0 = col - 1 + row * nx;

                    if (ndd[i0] == 2) {

                        col_idxs[nnz++] = idx - 4;
                        col_idxs[nnz++] = idx - 2;
                        col_idxs[nnz++] = idx;

                    }
                    else {

                        col_idxs[nnz++] = idx;

                    }

                }

            }

        }

    }

    row_ptrs[2 * n + ny] = nnz;


    continuityXdirKernel(n, nx, ny, dx, ndd, u, v, h, detadt, (double*)values, rhs, row_ptrs);

    momentumXdirKernel(n + ny, nx + 1, ny, dx, gra, mann, Du, ndd, u, v, eta, h, (double*)values, rhs, row_ptrs);


    ///////////////////////////////////////////////////////
    ///////////////////////////////////////////////////////


        FILE* fp;
    fp = fopen("values.bin", "wb");
    fwrite(values, sizeof(double), 6 * (2 * n + ny), fp);
    fclose(fp);

    fp = fopen("rhs.bin", "wb");
    fwrite(rhs, sizeof(double), (2 * n + ny), fp);
    fclose(fp);

    fp = fopen("row_ptrs.bin", "wb");
    fwrite(row_ptrs, sizeof(int), (2 * n + ny + 1), fp);
    fclose(fp);

    fp = fopen("col_idxs.bin", "wb");
    fwrite(col_idxs, sizeof(int64_t), 6 * (2 * n + ny), fp);
    fclose(fp);


    AMGX_matrix_upload_all_global(A,
                                   2 * n + ny, 2 * n + ny, nnz, 1, 1,
                                  row_ptrs, col_idxs, values, NULL,
                                  nrings, nrings, NULL);
    free(values);
    free(row_ptrs);
    free(col_idxs);
    /* generate the rhs and solution */
    void *h_x = malloc( (2 * n + ny) * sizeof_v_val);
    //for (int i = 0; i < 2 * n + ny; i++) { 
    //    if (i % 2 == 0) {
    //        ((double*)h_x)[i] = 0.0;
    //    }
    //    else {
    //        ((double*)h_x)[i] = 0.0;
    //    }   
    //}


    memset(h_x, 50.0, (2 * n + ny) * sizeof_v_val);

    /* set the connectivity information (for the vector) */
    AMGX_vector_bind(x, A);
    AMGX_vector_bind(b, A);
    /* upload the vector (and the connectivity information) */
    AMGX_vector_upload(x, 2 * n + ny, 1, h_x);
    AMGX_vector_upload(b, 2 * n + ny, 1, rhs);
    /* solver setup */
    //MPI barrier for stability (should be removed in practice to maximize performance)
    MPI_Barrier(amgx_mpi_comm);
    AMGX_solver_setup(solver, A);
    /* solver solve */
    //MPI barrier for stability (should be removed in practice to maximize performance)
    MPI_Barrier(amgx_mpi_comm);
    AMGX_solver_solve(solver, b, x);
    /* example of how to change parameters between non-linear iterations */
    //AMGX_config_add_parameters(&cfg, "config_version=2, default:tolerance=1e-12");
    //AMGX_solver_solve(solver, b, x);
    /* example of how to replace coefficients between non-linear iterations */
    //AMGX_matrix_replace_coefficients(A, n, nnz, values, diag);
    //AMGX_solver_setup(solver, A);
    //AMGX_solver_solve(solver, b, x);
    AMGX_solver_get_status(solver, &status);
    /* example of how to get (the local part of) the solution */
    //int sizeof_v_val;
    //sizeof_v_val = ((NVAMG_GET_MODE_VAL(NVAMG_VecPrecision, mode) == NVAMG_vecDouble))? sizeof(double): sizeof(float);
    // 
    void* SLN = malloc((2 * n + ny) * sizeof(double));
    AMGX_vector_download(x, SLN);

    fp = fopen("SLN.bin", "wb");
    fwrite(SLN, sizeof(double), (2 * n + ny), fp);
    fclose(fp);

    //free(result_host);
    /* destroy resources, matrix, vector and solver */
    AMGX_solver_destroy(solver);
    AMGX_vector_destroy(x);
    AMGX_vector_destroy(b);
    AMGX_matrix_destroy(A);
    AMGX_resources_destroy(rsrc);
    /* destroy config (need to use AMGX_SAFE_CALL after this point) */
    AMGX_SAFE_CALL(AMGX_config_destroy(cfg))
    /* shutdown and exit */
    AMGX_SAFE_CALL(AMGX_finalize_plugins())
    AMGX_SAFE_CALL(AMGX_finalize())
    /* close the library (if it was dynamically loaded) */
#ifdef AMGX_DYNAMIC_LOADING
    amgx_libclose(lib_handle);
#endif
    MPI_Finalize();
    CUDA_SAFE_CALL(cudaDeviceReset());
    return status;
}

/// below is the message returned

C:\STUDY\testamgx\CudaRuntime1\x64\Release>CudaRuntime1.exe -p 3 3 -c C:\AMGX\AMGX-2.3.0\core\configs\FGMRES_NOPREC.json
Process 0 selecting device 0
AMGX version 2.2.0.132-opensource
Built on Nov 23 2023, 14:14:00
Compiled with CUDA Runtime 11.8, using CUDA driver 12.2
Warning: No mode specified, using dDDI by default.
Cannot read file as JSON object, trying as AMGX config
Converting config string to current config version
Parsing configuration string: exception_handling=1 ;
idx = 0, nnz = 0
Using Normal MPI (Hostbuffer) communicator...
iter Mem Usage (GB) residual rate
--------------------------------------------------------------
Ini 0 5.141711e-03
0 0 5.141711e-03 1.0000
1 0.0000 5.141369e-03 0.9999
2 0.0000 5.141338e-03 1.0000
3 0.0000 5.141076e-03 0.9999
4 0.0000 5.141015e-03 1.0000
5 0.0000 5.140766e-03 1.0000
6 0.0000 5.140733e-03 1.0000
7 0.0000 5.140463e-03 0.9999
8 0.0000 5.140457e-03 1.0000
9 0.0000 5.140179e-03 0.9999
10 0.0000 5.140179e-03 1.0000
11 0.0000 5.139921e-03 0.9999
12 0.0000 5.139921e-03 1.0000
13 0.0000 5.139696e-03 1.0000
14 0.0000 5.139696e-03 1.0000
15 0.0000 5.139506e-03 1.0000
16 0.0000 5.139504e-03 1.0000
17 0.0000 5.139340e-03 1.0000
18 0.0000 5.139340e-03 1.0000
19 0.0000 5.139190e-03 1.0000
20 0.0000 5.139190e-03 1.0000
21 0.0000 5.139173e-03 1.0000
22 0.0000 5.139173e-03 1.0000
23 0.0000 5.139165e-03 1.0000
24 0.0000 5.139165e-03 1.0000
25 0.0000 5.139152e-03 1.0000
26 0.0000 5.139150e-03 1.0000
27 0.0000 5.139132e-03 1.0000
28 0.0000 5.139132e-03 1.0000
29 0.0000 5.139103e-03 1.0000
30 0.0000 5.139097e-03 1.0000
31 0.0000 5.139063e-03 1.0000
32 0.0000 5.139048e-03 1.0000
33 0.0000 5.139012e-03 1.0000
34 0.0000 5.138989e-03 1.0000
35 0.0000 5.138955e-03 1.0000
36 0.0000 5.138920e-03 1.0000
37 0.0000 5.138907e-03 1.0000
38 0.0000 5.138854e-03 1.0000
39 0.0000 5.138854e-03 1.0000
40 0.0000 5.138854e-03 1.0000
41 0.0000 5.138852e-03 1.0000
42 0.0000 5.138851e-03 1.0000
43 0.0000 5.138846e-03 1.0000
44 0.0000 5.138846e-03 1.0000
45 0.0000 5.138834e-03 1.0000
46 0.0000 5.138834e-03 1.0000
47 0.0000 5.138827e-03 1.0000
48 0.0000 5.138827e-03 1.0000
49 0.0000 5.138819e-03 1.0000
50 0.0000 5.138819e-03 1.0000
51 0.0000 5.138807e-03 1.0000
52 0.0000 5.138807e-03 1.0000
53 0.0000 5.138803e-03 1.0000
54 0.0000 5.138802e-03 1.0000
55 0.0000 5.138796e-03 1.0000
56 0.0000 5.138795e-03 1.0000
57 0.0000 5.138795e-03 1.0000
58 0.0000 5.138795e-03 1.0000
59 0.0000 5.138793e-03 1.0000
60 0.0000 5.138793e-03 1.0000
61 0.0000 5.138793e-03 1.0000
62 0.0000 5.138793e-03 1.0000
63 0.0000 5.138793e-03 1.0000
64 0.0000 5.138793e-03 1.0000
65 0.0000 5.138793e-03 1.0000
66 0.0000 5.138793e-03 1.0000
67 0.0000 5.138793e-03 1.0000
68 0.0000 5.138792e-03 1.0000
69 0.0000 5.138792e-03 1.0000
70 0.0000 5.138792e-03 1.0000
71 0.0000 5.138791e-03 1.0000
72 0.0000 5.138789e-03 1.0000
73 0.0000 5.138789e-03 1.0000
74 0.0000 5.138787e-03 1.0000
75 0.0000 5.138785e-03 1.0000
76 0.0000 5.138783e-03 1.0000
77 0.0000 5.138783e-03 1.0000
78 0.0000 5.138777e-03 1.0000
79 0.0000 5.138777e-03 1.0000
80 0.0000 5.138777e-03 1.0000
81 0.0000 5.138777e-03 1.0000
82 0.0000 5.138777e-03 1.0000
83 0.0000 5.138776e-03 1.0000
84 0.0000 5.138776e-03 1.0000
85 0.0000 5.138775e-03 1.0000
86 0.0000 5.138775e-03 1.0000
87 0.0000 5.138774e-03 1.0000
88 0.0000 5.138774e-03 1.0000
89 0.0000 5.138774e-03 1.0000
90 0.0000 5.138773e-03 1.0000
91 0.0000 5.138772e-03 1.0000
92 0.0000 5.138772e-03 1.0000
93 0.0000 5.138772e-03 1.0000
94 0.0000 5.138772e-03 1.0000
95 0.0000 5.138772e-03 1.0000
96 0.0000 5.138772e-03 1.0000
97 0.0000 5.138771e-03 1.0000
98 0.0000 5.138771e-03 1.0000
99 0.0000 5.138771e-03 1.0000
--------------------------------------------------------------
Total Iterations: 100
Avg Convergence Rate: 1.0000
Final Residual: 5.138771e-03
Total Reduction in Residual: 9.994281e-01
Maximum Memory Usage: 0.000 GB
--------------------------------------------------------------
Total Time: 0.0932096
setup: 8.192e-05 s
solve: 0.0931277 s
solve(per iteration): 0.000931277 s

/// below is the shape of the sparse matrix

///

[evstigneevnm]

Hello!
Can you provide a matrix itself? I'm collecting data for extreme test cases for my solver and that would be a great help! Plus, i might be able to give some advice on the configuration of AMGX if i have the matrix.

[fxuecnu]

@evstigneevnm Hello! Thank you for getting back to me.

I've attached the matrix information as requested. If you encounter any issues or have further questions, please don't hesitate to reach out.

matrix_info.zip