[petsc-users] How to efficiently fill in, in parallel, a PETSc matrix from a COO sparse matrix?
Matthew Knepley
knepley at gmail.com
Tue Jun 20 12:50:18 CDT 2023
On Tue, Jun 20, 2023 at 1:43 PM Diego Magela Lemos <diegomagela at usp.br>
wrote:
> 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*
>
I can answer this one.
> 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.
>
No. This was mostly intended for GPUs, where there is 1 process. If you
want to use multiple MPI processes, then each process can only introduce
some disjoint subset of the values. This is also how MatSetValues() works,
but it might not be as obvious.
Thanks,
Matt
> *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/>
>>
>>
>>
--
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/4e89ac5f/attachment-0001.html>
More information about the petsc-users
mailing list