Re: [petsc-users] Parallel matrix-vector product with Matrix Shell
Thanks again Jose for your prompt and very useful indication. By using
that, I could understand where the REAL problem was (and obviously it was
my fault).
Now everything works smoothly and produces the expected result.
All the best,
Mario
Il giorno sab 18 giu 2022 alle ore 08:27 Jose E. Roman
ha scritto:
> The initial vector of the Krylov method is by default a random vector,
> which is different when you change the number of processes. To avoid this,
> you can run with the undocumented option -bv_reproducible_random which will
> generate the same random initial vector irrespective of the number of
> processes.
>
> Alternatively, set an initial vector in your code with
> EPSSetInitialSpace(), see e.g.
> https://slepc.upv.es/documentation/current/src/eps/tutorials/ex5.c.html
>
> Jose
>
>
> > El 18 jun 2022, a las 8:13, Mario Rossi escribió:
> >
> > Dear Jose and Petsc users, I implemented the parallel matrix-vector
> product and it works meaning that it produces a result but it is different
> from the result produced with a single task.
> > Obviously, I could be wrong in my implementation but what puzzles me is
> that the input vector (x) to the product is different running with one and
> two tasks and this is from the very first iteration (so it can not be due
> to a previous error in the product).
> > I checked that X is different with one and two tasks with the following
> (naive) code
> > PetscErrorCode MatMult_TM(Mat A,Vec x,Vec y) {
> > void *ctx;
> > PetscInt nx /* ,lo,i,j*/;
> > const PetscScalar *px;
> > PetscScalar *py;
> > MPI_Comm comm;
> > PetscFunctionBeginUser;
> > PetscCall(MatShellGetContext(A,));
> > PetscCall(VecGetLocalSize(x,));
> > PetscCall(PetscObjectGetComm((PetscObject)A,));
> >
> > // nx = *(int*)ctx;
> > PetscCall(VecGetArrayRead(x,));
> > PetscCall(VecGetArray(y,));
> >
> > for(int i=0;i w[%d]=%f\n",myrank,i+offset,px[i],i+offset,w[i+offset]); }
> > PetscCall(MPI_Barrier(comm));
> > exit(0);
> > ..
> > }
> >
> > Then I reordered the output obtained with one and two tasks. The first
> part of the x vector is very similar (but not exactly the same) using one
> and two tasks but the second part (belonging to the second task) is pretty
> different
> > (here "offset" is offset=(n/size)*myrank;)
> > I create the matrix shell with
> >
>
> PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,,));
> > I am sure I am doing something wrong but I don't know what I need to
> look at.
> > Thanks in advance!
> > Mario
> >
> >
> > Il giorno ven 17 giu 2022 alle ore 17:56 Mario Rossi
> ha scritto:
> > Thanks a lot, Jose!
> > I looked at the eps folder (where I found the test8.c that has been my
> starting point) but I did not look at the nep folder (my fault!)
> > Thanks again,
> > Mario
> >
> > Il giorno ven 17 giu 2022 alle ore 17:34 Jose E. Roman <
> jro...@dsic.upv.es> ha scritto:
> > You can use VecGetOwnershipRange() to determine the range of global
> indices corresponding to the local portion of a vector, and VecGetArray()
> to access the values. In SLEPc, you can assume that X and Y will have the
> same parallel distribution.
> >
> > For an example of a shell matrix that implements the matrix-vector
> product in parallel, have a look at this:
> https://slepc.upv.es/documentation/current/src/nep/tutorials/ex21.c.html
> > It is a simple tridiagonal example, where neighborwise communication is
> done with two calls to MPI_Sendrecv().
> >
> > Jose
> >
> >
> > > El 17 jun 2022, a las 17:21, Mario Rossi escribió:
> > >
> > > I need to find the largest eigenvalues (say the first three) of a very
> large matrix and I am using
> > > a combination of PetSc and SLEPc. In particular, I am using a shell
> matrix. I wrote a "custom"
> > > matrix-vector product and everything works fine in serial (one task)
> mode for a "small" case.
> > > For the real case, I need multiple (at least 128) tasks for memory
> reasons so I need a parallel variant of the custom matrix-vector product. I
> know exactly how to write the parallel variant
> > > (in plain MPI) but I am, somehow, blocked because it is not clear to
> me what each task receives
> > > and what is expected to provide in the parallel matrix-vector product.
> > > More in detail, with a single task, the function receives the full X
> vector and is expected to provide the full Y vector resulting from Y=A*X.
> > > What does it happen with multiple tasks? If I understand correctly
> > > in the matrix shell definition, I can choose to split the matrix into
> blocks of rows so that the matrix-vector function should compute a block of
> elements of the vector Y but does it receive only the corresponding subset
> of the X (input vector)? (this is what I guess happens) and in output, does
> > > each task return its subset of elements of Y as if it were the whole
> array and then PetSc manages all the subsets? Is there anyone who has a
> working example of a
Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
Thank you for your reply. May I ask for some references on the order of the
dofs on PETSc's FE Space (especially high order elements)?
Thanks,
Zongze
Matthew Knepley 于2022年6月18日周六 20:02写道：
> On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang wrote:
>
>> In order to check if I made mistakes in the python code, I try to use c
>> code to show the issue on DMProjectCoordinates. The code and mesh file is
>> attached.
>> If the code is correct, there must be something wrong with
>> `DMProjectCoordinates` or `DMPlexCreateGmshFromFile` for high-order mesh.
>>
>
> Something is definitely wrong with high order, periodic simplices from
> Gmsh. We had not tested that case. I am at a conference and cannot look at
> it for a week.
> My suspicion is that the space we make when reading in the Gmsh
> coordinates does not match the values (wrong order).
>
> Thanks,
>
> Matt
>
>
>> The command and the output are listed below: (Obviously the bounding box
>> is changed.)
>> ```
>> $ ./test_gmsh_load_2rd -filename cube-p2.msh -old_fe_view -new_fe_view
>> Old Bounding Box:
>> 0: lo = 0. hi = 1.
>> 1: lo = 0. hi = 1.
>> 2: lo = 0. hi = 1.
>> PetscFE Object: OldCoordinatesFE 1 MPI processes
>> type: basic
>> Basic Finite Element in 3 dimensions with 3 components
>> PetscSpace Object: P2 1 MPI processes
>> type: sum
>> Space in 3 variables with 3 components, size 30
>> Sum space of 3 concatenated subspaces (all identical)
>> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
>> type: poly
>> Space in 3 variables with 1 components, size 10
>> Polynomial space of degree 2
>> PetscDualSpace Object: P2 1 MPI processes
>> type: lagrange
>> Dual space with 3 components, size 30
>> Discontinuous Lagrange dual space
>> Quadrature of order 5 on 27 points (dim 3)
>> PetscFE Object: NewCoordinatesFE 1 MPI processes
>> type: basic
>> Basic Finite Element in 3 dimensions with 3 components
>> PetscSpace Object: P2 1 MPI processes
>> type: sum
>> Space in 3 variables with 3 components, size 30
>> Sum space of 3 concatenated subspaces (all identical)
>> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
>> type: poly
>> Space in 3 variables with 1 components, size 10
>> Polynomial space of degree 2
>> PetscDualSpace Object: P2 1 MPI processes
>> type: lagrange
>> Dual space with 3 components, size 30
>> Continuous Lagrange dual space
>> Quadrature of order 5 on 27 points (dim 3)
>> New Bounding Box:
>> 0: lo = 2.5624e-17 hi = 8.
>> 1: lo = -9.23372e-17 hi = 7.
>> 2: lo = 2.72091e-17 hi = 8.5
>> ```
>>
>> Thanks,
>> Zongze
>>
>> Zongze Yang 于2022年6月17日周五 14:54写道：
>>
>>> I tried the projection operation. However, it seems that the projection
>>> gives the wrong solution. After projection, the bounding box is changed!
>>> See logs below.
>>>
>>> First, I patch the petsc4py by adding `DMProjectCoordinates`:
>>> ```
>>> diff --git a/src/binding/petsc4py/src/PETSc/DM.pyx
>>> b/src/binding/petsc4py/src/PETSc/DM.pyx
>>> index d8a58d183a..dbcdb280f1 100644
>>> --- a/src/binding/petsc4py/src/PETSc/DM.pyx
>>> +++ b/src/binding/petsc4py/src/PETSc/DM.pyx
>>> @@ -307,6 +307,12 @@ cdef class DM(Object):
>>> PetscINCREF(c.obj)
>>> return c
>>>
>>> +def projectCoordinates(self, FE fe=None):
>>> +if fe is None:
>>> +CHKERR( DMProjectCoordinates(self.dm, NULL) )
>>> +else:
>>> +CHKERR( DMProjectCoordinates(self.dm, fe.fe) )
>>> +
>>> def getBoundingBox(self):
>>> cdef PetscInt i,dim=0
>>> CHKERR( DMGetCoordinateDim(self.dm, ) )
>>> diff --git a/src/binding/petsc4py/src/PETSc/petscdm.pxi
>>> b/src/binding/petsc4py/src/PETSc/petscdm.pxi
>>> index 514b6fa472..c778e39884 100644
>>> --- a/src/binding/petsc4py/src/PETSc/petscdm.pxi
>>> +++ b/src/binding/petsc4py/src/PETSc/petscdm.pxi
>>> @@ -90,6 +90,7 @@ cdef extern from * nogil:
>>> int DMGetCoordinateDim(PetscDM,PetscInt*)
>>> int DMSetCoordinateDim(PetscDM,PetscInt)
>>> int DMLocalizeCoordinates(PetscDM)
>>> +int DMProjectCoordinates(PetscDM, PetscFE)
>>>
>>> int DMCreateInterpolation(PetscDM,PetscDM,PetscMat*,PetscVec*)
>>> int DMCreateInjection(PetscDM,PetscDM,PetscMat*)
>>> ```
>>>
>>> Then in python, I load a mesh and project the coordinates to P2:
>>> ```
>>> import firedrake as fd
>>> from firedrake.petsc import PETSc
>>>
>>> # plex = fd.mesh._from_gmsh('test-fd-load-p2.msh')
>>> plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')
>>> print('old bbox:', plex.getBoundingBox())
>>>
>>> dim = plex.getDimension()
>>> # (dim, nc, isSimplex, k,
>>> qorder, comm=None)
>>> fe_new = PETSc.FE().createLagrange(dim, dim, True, 2,
>>> PETSc.DETERMINE)
>>> plex.projectCoordinates(fe_new)
>>> fe_new.view()
>>>
>>> print('new bbox:', plex.getBoundingBox())
>>> ```
>>>
>>> The output is (The bounding box is
Re: [petsc-users] How to find the map between the high order coordinates of DMPlex and vertex numbering?
On Sat, Jun 18, 2022 at 2:16 AM Zongze Yang wrote:
> In order to check if I made mistakes in the python code, I try to use c
> code to show the issue on DMProjectCoordinates. The code and mesh file is
> attached.
> If the code is correct, there must be something wrong with
> `DMProjectCoordinates` or `DMPlexCreateGmshFromFile` for high-order mesh.
>
Something is definitely wrong with high order, periodic simplices from
Gmsh. We had not tested that case. I am at a conference and cannot look at
it for a week.
My suspicion is that the space we make when reading in the Gmsh coordinates
does not match the values (wrong order).
Thanks,
Matt
> The command and the output are listed below: (Obviously the bounding box
> is changed.)
> ```
> $ ./test_gmsh_load_2rd -filename cube-p2.msh -old_fe_view -new_fe_view
> Old Bounding Box:
> 0: lo = 0. hi = 1.
> 1: lo = 0. hi = 1.
> 2: lo = 0. hi = 1.
> PetscFE Object: OldCoordinatesFE 1 MPI processes
> type: basic
> Basic Finite Element in 3 dimensions with 3 components
> PetscSpace Object: P2 1 MPI processes
> type: sum
> Space in 3 variables with 3 components, size 30
> Sum space of 3 concatenated subspaces (all identical)
> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
> type: poly
> Space in 3 variables with 1 components, size 10
> Polynomial space of degree 2
> PetscDualSpace Object: P2 1 MPI processes
> type: lagrange
> Dual space with 3 components, size 30
> Discontinuous Lagrange dual space
> Quadrature of order 5 on 27 points (dim 3)
> PetscFE Object: NewCoordinatesFE 1 MPI processes
> type: basic
> Basic Finite Element in 3 dimensions with 3 components
> PetscSpace Object: P2 1 MPI processes
> type: sum
> Space in 3 variables with 3 components, size 30
> Sum space of 3 concatenated subspaces (all identical)
> PetscSpace Object: sum component (sumcomp_) 1 MPI processes
> type: poly
> Space in 3 variables with 1 components, size 10
> Polynomial space of degree 2
> PetscDualSpace Object: P2 1 MPI processes
> type: lagrange
> Dual space with 3 components, size 30
> Continuous Lagrange dual space
> Quadrature of order 5 on 27 points (dim 3)
> New Bounding Box:
> 0: lo = 2.5624e-17 hi = 8.
> 1: lo = -9.23372e-17 hi = 7.
> 2: lo = 2.72091e-17 hi = 8.5
> ```
>
> Thanks,
> Zongze
>
> Zongze Yang 于2022年6月17日周五 14:54写道：
>
>> I tried the projection operation. However, it seems that the projection
>> gives the wrong solution. After projection, the bounding box is changed!
>> See logs below.
>>
>> First, I patch the petsc4py by adding `DMProjectCoordinates`:
>> ```
>> diff --git a/src/binding/petsc4py/src/PETSc/DM.pyx
>> b/src/binding/petsc4py/src/PETSc/DM.pyx
>> index d8a58d183a..dbcdb280f1 100644
>> --- a/src/binding/petsc4py/src/PETSc/DM.pyx
>> +++ b/src/binding/petsc4py/src/PETSc/DM.pyx
>> @@ -307,6 +307,12 @@ cdef class DM(Object):
>> PetscINCREF(c.obj)
>> return c
>>
>> +def projectCoordinates(self, FE fe=None):
>> +if fe is None:
>> +CHKERR( DMProjectCoordinates(self.dm, NULL) )
>> +else:
>> +CHKERR( DMProjectCoordinates(self.dm, fe.fe) )
>> +
>> def getBoundingBox(self):
>> cdef PetscInt i,dim=0
>> CHKERR( DMGetCoordinateDim(self.dm, ) )
>> diff --git a/src/binding/petsc4py/src/PETSc/petscdm.pxi
>> b/src/binding/petsc4py/src/PETSc/petscdm.pxi
>> index 514b6fa472..c778e39884 100644
>> --- a/src/binding/petsc4py/src/PETSc/petscdm.pxi
>> +++ b/src/binding/petsc4py/src/PETSc/petscdm.pxi
>> @@ -90,6 +90,7 @@ cdef extern from * nogil:
>> int DMGetCoordinateDim(PetscDM,PetscInt*)
>> int DMSetCoordinateDim(PetscDM,PetscInt)
>> int DMLocalizeCoordinates(PetscDM)
>> +int DMProjectCoordinates(PetscDM, PetscFE)
>>
>> int DMCreateInterpolation(PetscDM,PetscDM,PetscMat*,PetscVec*)
>> int DMCreateInjection(PetscDM,PetscDM,PetscMat*)
>> ```
>>
>> Then in python, I load a mesh and project the coordinates to P2:
>> ```
>> import firedrake as fd
>> from firedrake.petsc import PETSc
>>
>> # plex = fd.mesh._from_gmsh('test-fd-load-p2.msh')
>> plex = fd.mesh._from_gmsh('test-fd-load-p2-rect.msh')
>> print('old bbox:', plex.getBoundingBox())
>>
>> dim = plex.getDimension()
>> # (dim, nc, isSimplex, k,
>> qorder, comm=None)
>> fe_new = PETSc.FE().createLagrange(dim, dim, True, 2,
>> PETSc.DETERMINE)
>> plex.projectCoordinates(fe_new)
>> fe_new.view()
>>
>> print('new bbox:', plex.getBoundingBox())
>> ```
>>
>> The output is (The bounding box is changed!)
>> ```
>>
>> old bbox: ((0.0, 1.0), (0.0, 1.0), (0.0, 1.0))
>> PetscFE Object: P2 1 MPI processes
>> type: basic
>> Basic Finite Element in 3 dimensions with 3 components
>> PetscSpace Object: P2 1 MPI processes
>> type: sum
>> Space in 3 variables with 3 components, size 30
>> Sum space of 3 concatenated
Re: [petsc-users] Parallel matrix-vector product with Matrix Shell
The initial vector of the Krylov method is by default a random vector, which is
different when you change the number of processes. To avoid this, you can run
with the undocumented option -bv_reproducible_random which will generate the
same random initial vector irrespective of the number of processes.
Alternatively, set an initial vector in your code with EPSSetInitialSpace(),
see e.g. https://slepc.upv.es/documentation/current/src/eps/tutorials/ex5.c.html
Jose
> El 18 jun 2022, a las 8:13, Mario Rossi escribió:
>
> Dear Jose and Petsc users, I implemented the parallel matrix-vector product
> and it works meaning that it produces a result but it is different from the
> result produced with a single task.
> Obviously, I could be wrong in my implementation but what puzzles me is that
> the input vector (x) to the product is different running with one and two
> tasks and this is from the very first iteration (so it can not be due to a
> previous error in the product).
> I checked that X is different with one and two tasks with the following
> (naive) code
> PetscErrorCode MatMult_TM(Mat A,Vec x,Vec y) {
> void *ctx;
> PetscInt nx /* ,lo,i,j*/;
> const PetscScalar *px;
> PetscScalar *py;
> MPI_Comm comm;
> PetscFunctionBeginUser;
> PetscCall(MatShellGetContext(A,));
> PetscCall(VecGetLocalSize(x,));
> PetscCall(PetscObjectGetComm((PetscObject)A,));
>
> // nx = *(int*)ctx;
> PetscCall(VecGetArrayRead(x,));
> PetscCall(VecGetArray(y,));
>
> for(int i=0;i w[%d]=%f\n",myrank,i+offset,px[i],i+offset,w[i+offset]); }
> PetscCall(MPI_Barrier(comm));
> exit(0);
> ..
> }
>
> Then I reordered the output obtained with one and two tasks. The first part
> of the x vector is very similar (but not exactly the same) using one and two
> tasks but the second part (belonging to the second task) is pretty different
> (here "offset" is offset=(n/size)*myrank;)
> I create the matrix shell with
>
> PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,,));
> I am sure I am doing something wrong but I don't know what I need to look at.
> Thanks in advance!
> Mario
>
>
> Il giorno ven 17 giu 2022 alle ore 17:56 Mario Rossi ha
> scritto:
> Thanks a lot, Jose!
> I looked at the eps folder (where I found the test8.c that has been my
> starting point) but I did not look at the nep folder (my fault!)
> Thanks again,
> Mario
>
> Il giorno ven 17 giu 2022 alle ore 17:34 Jose E. Roman
> ha scritto:
> You can use VecGetOwnershipRange() to determine the range of global indices
> corresponding to the local portion of a vector, and VecGetArray() to access
> the values. In SLEPc, you can assume that X and Y will have the same parallel
> distribution.
>
> For an example of a shell matrix that implements the matrix-vector product in
> parallel, have a look at this:
> https://slepc.upv.es/documentation/current/src/nep/tutorials/ex21.c.html
> It is a simple tridiagonal example, where neighborwise communication is done
> with two calls to MPI_Sendrecv().
>
> Jose
>
>
> > El 17 jun 2022, a las 17:21, Mario Rossi escribió:
> >
> > I need to find the largest eigenvalues (say the first three) of a very
> > large matrix and I am using
> > a combination of PetSc and SLEPc. In particular, I am using a shell matrix.
> > I wrote a "custom"
> > matrix-vector product and everything works fine in serial (one task) mode
> > for a "small" case.
> > For the real case, I need multiple (at least 128) tasks for memory reasons
> > so I need a parallel variant of the custom matrix-vector product. I know
> > exactly how to write the parallel variant
> > (in plain MPI) but I am, somehow, blocked because it is not clear to me
> > what each task receives
> > and what is expected to provide in the parallel matrix-vector product.
> > More in detail, with a single task, the function receives the full X vector
> > and is expected to provide the full Y vector resulting from Y=A*X.
> > What does it happen with multiple tasks? If I understand correctly
> > in the matrix shell definition, I can choose to split the matrix into
> > blocks of rows so that the matrix-vector function should compute a block of
> > elements of the vector Y but does it receive only the corresponding subset
> > of the X (input vector)? (this is what I guess happens) and in output, does
> > each task return its subset of elements of Y as if it were the whole array
> > and then PetSc manages all the subsets? Is there anyone who has a working
> > example of a parallel matrix-vector product for matrix shell?
> > Thanks in advance for any help you can provide!
> > Mario
> > i
> >
>
Re: [petsc-users] Parallel matrix-vector product with Matrix Shell
Dear Jose and Petsc users, I implemented the parallel matrix-vector product
and it works meaning that it produces a result but it is different from the
result produced with a single task.
Obviously, I could be wrong in my implementation but what puzzles me is
that the *input *vector (x) to the product is different running with one
and two tasks and this is from the very first iteration (so it can not be
due to a previous error in the product).
I checked that X is different with one and two tasks with the following
(naive) code
PetscErrorCode MatMult_TM(Mat A,Vec x,Vec y) {
void *ctx;
PetscInt nx /* ,lo,i,j*/;
const PetscScalar *px;
PetscScalar *py;
MPI_Comm comm;
PetscFunctionBeginUser;
PetscCall(MatShellGetContext(A,));
PetscCall(VecGetLocalSize(x,));
PetscCall(PetscObjectGetComm((PetscObject)A,));
// nx = *(int*)ctx;
PetscCall(VecGetArrayRead(x,));
PetscCall(VecGetArray(y,));
for(int i=0;i ha
scritto:
> Thanks a lot, Jose!
> I looked at the eps folder (where I found the test8.c that has been my
> starting point) but I did not look at the nep folder (my fault!)
> Thanks again,
> Mario
>
> Il giorno ven 17 giu 2022 alle ore 17:34 Jose E. Roman
> ha scritto:
>
>> You can use VecGetOwnershipRange() to determine the range of global
>> indices corresponding to the local portion of a vector, and VecGetArray()
>> to access the values. In SLEPc, you can assume that X and Y will have the
>> same parallel distribution.
>>
>> For an example of a shell matrix that implements the matrix-vector
>> product in parallel, have a look at this:
>> https://slepc.upv.es/documentation/current/src/nep/tutorials/ex21.c.html
>> It is a simple tridiagonal example, where neighborwise communication is
>> done with two calls to MPI_Sendrecv().
>>
>> Jose
>>
>>
>> > El 17 jun 2022, a las 17:21, Mario Rossi escribió:
>> >
>> > I need to find the largest eigenvalues (say the first three) of a very
>> large matrix and I am using
>> > a combination of PetSc and SLEPc. In particular, I am using a shell
>> matrix. I wrote a "custom"
>> > matrix-vector product and everything works fine in serial (one task)
>> mode for a "small" case.
>> > For the real case, I need multiple (at least 128) tasks for memory
>> reasons so I need a parallel variant of the custom matrix-vector product. I
>> know exactly how to write the parallel variant
>> > (in plain MPI) but I am, somehow, blocked because it is not clear to me
>> what each task receives
>> > and what is expected to provide in the parallel matrix-vector product.
>> > More in detail, with a single task, the function receives the full X
>> vector and is expected to provide the full Y vector resulting from Y=A*X.
>> > What does it happen with multiple tasks? If I understand correctly
>> > in the matrix shell definition, I can choose to split the matrix into
>> blocks of rows so that the matrix-vector function should compute a block of
>> elements of the vector Y but does it receive only the corresponding subset
>> of the X (input vector)? (this is what I guess happens) and in output, does
>> > each task return its subset of elements of Y as if it were the whole
>> array and then PetSc manages all the subsets? Is there anyone who has a
>> working example of a parallel matrix-vector product for matrix shell?
>> > Thanks in advance for any help you can provide!
>> > Mario
>> > i
>> >
>>
>>
