[petsc-users] How to efficiently fill in, in parallel, a PETSc matrix from a COO sparse matrix?

[Diego Magela Lemos]

Using all recommended approaches it worked!
Thank you all in advance.

Now, I'm facing problems when solving a linear system using each approach.

*COO approach*

Using MatSetPreallocationCOO and then MatSetValuesCOO, I'm able to fill in
the matrix when running with 1 MPI process.
But, if I run with more than one MPI process, the values entries are
multiplied by the number of MPI processes being used.
Is this behavior correct?

Consider the following code:

// fill_in_matrix.cc

static char help[] = "Fill in a parallel COO format sparse matrix.";

#include <petsc.h>
#include <vector>

int main(int argc, char **args)
{
    std::vector<PetscInt> coo_i{0, 0, 1, 2, 3, 4};
    std::vector<PetscInt> coo_j{0, 0, 1, 2, 3, 4};
    std::vector<PetscScalar> coo_v{2, -1, 2, 3, 4, 5};

    Mat A;

    PetscInt size = 5;

    PetscCall(PetscInitialize(&argc, &args, NULL, help));

    // Create matrix
    PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, size, size));
    PetscCall(MatSetFromOptions(A));

    // Populate matrix
    PetscCall(MatSetPreallocationCOO(A, coo_v.size(), coo_i.data(),
coo_j.data()));
    PetscCall(MatSetValuesCOO(A, coo_v.data(), ADD_VALUES));

    // View matrix
    PetscCall(MatView(A, PETSC_VIEWER_STDERR_WORLD));

    PetscCall(MatDestroy(&A));

    PetscCall(PetscFinalize());
    return 0;
}

When running with petscmpiexec -n 1 ./fill_in_matrix, I got

Mat Object: 1 MPI process
>   type: seqaij
> row 0: (0, 1.)
> row 1: (1, 2.)
> row 2: (2, 3.)
> row 3: (3, 4.)
> row 4: (4, 5.)


Which is a correct result. But, when running it with petscmpiexec -n 2
./fill_in_matrix, I get

Mat Object: 2 MPI process
>   type: mpiaij
> row 0: (0, 2.)
> row 1: (1, 4.)
> row 2: (2, 6.)
> row 3: (3, 8.)
> row 4: (4, 10.)


The matrix entries are multiplied by 2, that is, the number of processes
used to execute the code.

*MatSetValues approach*

I obtain the same behavior when filling in the matrix by using MatSetValues

static char help[] = "Fill in a parallel COO format sparse matrix.";

// fill_in_matrix.cc

#include <petsc.h>
#include <vector>

int main(int argc, char **args)
{
    std::vector<PetscInt> coo_i{0, 0, 1, 2, 3, 4};
    std::vector<PetscInt> coo_j{0, 0, 1, 2, 3, 4};
    std::vector<PetscScalar> coo_v{2, -1, 2, 3, 4, 5};

    Mat A;
    PetscInt size = 5;

    PetscCall(PetscInitialize(&argc, &args, NULL, help));

    // Create matrix
    PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, size, size));
    PetscCall(MatSetFromOptions(A));

    // Populate matrix
    for (size_t i = 0; i < coo_v.size(); i++)
        PetscCall(MatSetValues(A, 1, &coo_i.at(i), 1, &coo_j.at(i), &
coo_v.at(i), ADD_VALUES));

    PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    // View matrix
    PetscCall(MatView(A, PETSC_VIEWER_STDERR_WORLD));

    PetscCall(MatDestroy(&A));

    PetscCall(PetscFinalize());
    return 0;
}

When solving a linear system, I get the correct answer no matter the number
of MPI processes when using MatSetValues approach.
The same is not true when using COO approach, whose result is only correct
when using 1 MPI process.

static char help[] = "Fill in a parallel COO format sparse matrix and solve
a linear system.";

#include <petsc.h>
#include <vector>

int main(int argc, char **args)
{
    std::vector<PetscInt> coo_i{0, 0, 1, 2, 3, 4};
    std::vector<PetscInt> coo_j{0, 0, 1, 2, 3, 4};
    std::vector<PetscScalar> coo_v{2, -1, 2, 3, 4, 5};

    Mat A;
    Vec B, X, U;
    KSP ksp;
    PC pc;
    PetscReal norm; // norm of solution error
    PetscInt its;

    PetscInt size = 5;

    PetscCall(PetscInitialize(&argc, &args, NULL, help));

    // Create matrix

    PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
    PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, size, size));
    PetscCall(MatSetFromOptions(A));


    // Populate matrix

    // COO
    PetscCall(MatSetPreallocationCOO(A, coo_v.size(), coo_i.data(),
coo_j.data()));
    PetscCall(MatSetValuesCOO(A, coo_v.data(), ADD_VALUES));

    // MatSetValues for-loop
    // for (size_t i = 0; i < coo_v.size(); i++)
    //     PetscCall(MatSetValues(A, 1, &coo_i.at(i), 1, &coo_j.at(i), &
coo_v.at(i), ADD_VALUES));

    // PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    // PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

    // View matrix
    PetscCall(MatView(A, PETSC_VIEWER_STDERR_WORLD));

    // Create vector B
    PetscCall(VecCreate(PETSC_COMM_WORLD, &B));
    PetscCall(VecSetSizes(B, PETSC_DECIDE, size));
    PetscCall(VecSetFromOptions(B));
    PetscCall(VecSetUp(B));

    // Populate vector
    PetscCall(VecSetValues(B, coo_i.size(), coo_i.data(), coo_v.data(),
ADD_VALUES));
    PetscCall(VecAssemblyBegin(B));
    PetscCall(VecAssemblyEnd(B));

    // View vector
    PetscCall(VecView(B, PETSC_VIEWER_STDERR_WORLD));

    // Define solution and auxiliary vector
    PetscCall(VecDuplicate(B, &X));
    PetscCall(VecDuplicate(B, &U));
    PetscCall(VecSet(U, 1.0));

    // Create solver
    PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
    PetscCall(KSPSetOperators(ksp, A, A));
    PetscCall(KSPGetPC(ksp, &pc));
    PetscCall(PCSetType(pc, PCJACOBI));
    PetscCall(KSPSetFromOptions(ksp));
    PetscCall(KSPSetTolerances(ksp, 1.e-5, PETSC_DEFAULT, PETSC_DEFAULT,
PETSC_DEFAULT));

    // Solve
    PetscCall(KSPSolve(ksp, B, X));

    // View solution vector
    PetscCall(VecView(X, PETSC_VIEWER_STDERR_WORLD));

    // Verify the solution
    PetscCall(VecAXPY(X, -1.0, U));
    PetscCall(VecNorm(X, NORM_2, &norm));
    PetscCall(KSPGetIterationNumber(ksp, &its));
    PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations
%" PetscInt_FMT "\n", (double)norm, its));

    PetscCall(MatDestroy(&A));
    PetscCall(VecDestroy(&B));
    PetscCall(VecDestroy(&X));
    PetscCall(VecDestroy(&U));
    PetscCall(KSPDestroy(&ksp));

    PetscCall(PetscFinalize());
    return 0;
}

Why am I getting wrong results using the COO approach with more than one
MPI process?

Em ter., 20 de jun. de 2023 às 13:13, Barry Smith <bsmith at petsc.dev>
escreveu:

>
>   Since you have 6 entries that needed to be added to the matrix you will
> need to call MatSetValues() six time for the six entries.
>
> On Jun 20, 2023, at 11:06 AM, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Tue, Jun 20, 2023 at 10:55 AM Diego Magela Lemos via petsc-users <
> petsc-users at mcs.anl.gov> wrote:
>
>> Considering, for instance, the following COO sparse matrix format, with
>> repeated indices:
>>
>> std::vector<size_t> rows{0, 0, 1, 2, 3, 4};
>> std::vector<size_t> cols{0, 0, 1, 2, 3, 4};
>> std::vector<double> values{2, -1, 2, 3, 4, 5};
>>
>> that represents a 5x5 diagonal matrix A.
>>
>> So far, the code that I have is:
>>
>> // fill_in_matrix.cc
>> static char help[] = "Fill in a parallel COO format sparse matrix.";
>> #include <petsc.h>#include <vector>
>> int main(int argc, char **args){
>>     Mat A;
>>     PetscInt m = 5, i, Istart, Iend;
>>
>>     PetscCall(PetscInitialize(&argc, &args, NULL, help));
>>
>>     PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
>>     PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, m));
>>     PetscCall(MatSetFromOptions(A));
>>     PetscCall(MatSetUp(A));
>>     PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
>>
>>     std::vector<PetscInt> II{0, 0, 1, 2, 3, 4};
>>     std::vector<PetscInt> JJ{0, 0, 1, 2, 3, 4};
>>     std::vector<PetscScalar> XX{2, -1, 2, 3, 4, 5};
>>
>>     for (i = Istart; i < Iend; i++)
>>         PetscCall(MatSetValues(A, 1, &II.at(i), 1, &JJ.at(i), &XX.at(i), ADD_VALUES));
>>
>>     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
>>     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
>>     PetscCall(MatView(A, PETSC_VIEWER_STDERR_WORLD));
>>
>>     PetscCall(MatDestroy(&A));
>>     PetscCall(PetscFinalize());
>>     return 0;
>> }
>>
>> When running it with
>>
>> petscmpiexec -n 4 ./fill_in_matrix
>>
>>
>> I get
>>
>>
>>  Mat Object: 4 MPI processes
>>
>>   type: mpiaij
>> row 0: (0, 1.)
>> row 1: (1, 2.)
>> row 2: (2, 3.)
>> row 3: (3, 4.)
>> row 4:
>>
>>
>> Which is missing the entry of the last row.
>> What am I missing? Even better, which would be the best way to fill in this matrix?
>>
>> We have a new interface for this:
>
>   https://petsc.org/main/manualpages/Mat/MatSetValuesCOO/
>
>   Thanks,
>
>      Matt
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230620/c8812220/attachment-0001.html>