You should partition the entries so each entry is submitted by only one
process. Note that duplicate entries (on the same or different proceses) are
summed as you've seen. For example, in finite elements, it's typical to
partition the elements and each process submits entries from its elements.
Diego Magela Lemos via petsc-users <petsc-users@mcs.anl.gov> writes:
> 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 <bsm...@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 <knep...@gmail.com> wrote:
>>
>> On Tue, Jun 20, 2023 at 10:55 AM Diego Magela Lemos via petsc-users <
>> petsc-users@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/>
>>
>>
>>
