Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

HYPRE_ParCSRPCGSolve returns unexpected result #1187

Open
hhrecode001 opened this issue Dec 3, 2024 · 0 comments
Open

HYPRE_ParCSRPCGSolve returns unexpected result #1187

hhrecode001 opened this issue Dec 3, 2024 · 0 comments

Comments

@hhrecode001
Copy link

I'm trying to use Hypre to solve a system of linear equations:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "HYPRE_krylov.h"
#include "HYPRE.h"
#include "HYPRE_parcsr_ls.h"

int main(int argc, char* argv[])
{
    int i;
    int myid, num_procs;
    int N, n;

    int ilower, iupper;
    int local_size;

    HYPRE_IJMatrix A;
    HYPRE_ParCSRMatrix parcsr_A;
    HYPRE_IJVector b;
    HYPRE_ParVector par_b;
    HYPRE_IJVector x;
    HYPRE_ParVector par_x;

    HYPRE_Solver solver, precond;

    // Initialize
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    HYPRE_Initialize();

    N = 6;
    local_size = N / num_procs;

    ilower = local_size * myid;
    iupper = local_size * (myid + 1) - 1;
    local_size = iupper - ilower + 1;

    HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);
    HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);
    HYPRE_IJMatrixInitialize(A);

    // Fill the matrix anti-diagonal
    {
        int nnz;
        double values[1];
        int cols[1];
        int tmp[2];

        for (i = ilower; i <= iupper; i++)
        {
            nnz = 0;
            cols[nnz] = N - i - 1;
            values[nnz] = 1;
            nnz++;

            tmp[0] = nnz;
            tmp[1] = i;
            HYPRE_IJMatrixSetValues(A, 1, &tmp[0], &tmp[1], cols, values);
        }
    }

    HYPRE_IJMatrixAssemble(A);
    HYPRE_IJMatrixGetObject(A, (void**)&parcsr_A);

    // Create the rhs and solution
    HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &b);
    HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
    HYPRE_IJVectorInitialize(b);

    HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &x);
    HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR);
    HYPRE_IJVectorInitialize(x);

    {
        double* rhs_values, * x_values;
        int* rows;

        rhs_values = (double*)calloc(local_size, sizeof(double));
        x_values = (double*)calloc(local_size, sizeof(double));
        rows = (int*)calloc(local_size, sizeof(int));

        for (i = 0; i < local_size; i++)
        {
            rhs_values[i] = i + 1;
            x_values[i] = 0.0;
            rows[i] = ilower + i;
        }

        HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values);
        HYPRE_IJVectorSetValues(x, local_size, rows, x_values);

        free(x_values);
        free(rhs_values);
        free(rows);
    }

    HYPRE_IJVectorAssemble(b);
    HYPRE_IJVectorGetObject(b, (void**)&par_b);

    HYPRE_IJVectorAssemble(x);
    HYPRE_IJVectorGetObject(x, (void**)&par_x);

    HYPRE_IJMatrixPrint(A, "IJ.out.A");
    HYPRE_IJVectorPrint(b, "IJ.out.b");

    // Solve the system using PCG
    HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
    HYPRE_PCGSetMaxIter(solver, 100);
    HYPRE_PCGSetTol(solver, 1e-6);
    HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
    HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);

    // print the local solution
    int nvalues = local_size;
    int* rows = (int*)calloc(nvalues, sizeof(int));
    double* values = (double*)calloc(nvalues, sizeof(double));
    for (i = 0; i < nvalues; i++)
    {
        rows[i] = ilower + i;
    }
    HYPRE_IJVectorGetValues(x, nvalues, rows, values);

    for (int i = 0; i < nvalues; i++) {
        printf("Rank %d, x[%d] = %lf\n", myid, rows[i], values[i]);
    }

    // Finalize HYPRE
    HYPRE_Finalize();
    MPI_Finalize();
    return (0);
}

With N=6, the matrix A will be

0 0 0 0 0 1
0 0 0 0 1 0
0 0 0 1 0 0
0 0 1 0 0 0
0 1 0 0 0 0
1 0 0 0 0 0

And vector B = 1 2 3 4 5 6
The expected result X will be 6 5 4 3 2 1

But the result I got is very strange
With mpi world size = 1, X = 1.625 3.250 4.875 6.500 8.125 9.750
With mpi world size = 2, X = 1.4 2.8 4.2 1.4 2.8 4.2

I tried to modify this line from
HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);
to
HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, 0, N - 1, &A);
but the results are also strange.

I'm using Visual Studio 2022 on Windows with Intel MPI library.
The hypre install cmake config is: cmake -S ".." -B . -G "Visual Studio 17 2022" -DHYPRE_ENABLE_ONEMKLBLAS=ON -DHYPRE_ENABLE_ONEMKLRAND=ON -DHYPRE_ENABLE_ONEMKLSPARSE=ON -DHYPRE_HAVE_MPI=ON -DHYPRE_WITH_MPI=ON -DHYPRE_USING_HOST_MEMORY=ON -DHYPRE_USING_HYPRE_BLAS=ON -DHYPRE_USING_HYPRE_LAPACK=ON

@hhrecode001 hhrecode001 changed the title Solver returns unexpected result HYPRE_ParCSRPCGSolve returns unexpected result Dec 3, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant