September 2008 - petsc-dev - mailman3.cels.anl.gov

Re: User-difined PC
by Lisandro Dalcin 22 Sep '08

22 Sep '08

BTW, I'm pretty sure that long, long time ago I asked for this, and some of you rejected the idea, but I cannot remember the rationale right now. The current state of PCShell somewhat complicates the access to internal PC stuff (eg. linear operators). On Mon, Sep 22, 2008 at 1:08 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > On Sep 22, 2008, at 7:59 AM, Dave May wrote: > >> >> >> Why is the first pointer in PCSHELL operations the user context >> whereas in MATSHELL the first pointer is of type Mat ?? Seems slightly >> inconsistent. >> > Good question, these were done at different times and only later > did we realize the inconsistency. This is something that should be fixed in > petsc-dev. > > Barry > > I guess the first argument for the PCSHELL should become the pc. > >> >> >> On Mon, Sep 22, 2008 at 8:31 PM, Michel Cancelliere >> <fernandez858(a)gmail.com> wrote: >>> >>> Hello, >>> >>> I have problems with the implementation of a user-defined pc, basically >>> my >>> program is in a cycle, for with which he is called by matlab each time it >>> seeks to solve a linear system, the problem is that in the first >>> iteration >>> of the for-cycle the preconditioner works very well, then fails to >>> convergence. It may be some parameters that I setting wrong? >>> I'am attaching my program code. >>> >>> Thank you, >>> >>> Michel >>> >> > > -- Lisandro Dalcï¿½n --------------- Centro Internacional de Mï¿½todos Computacionales en Ingenierï¿½a (CIMEC) Instituto de Desarrollo Tecnolï¿½gico para la Industria Quï¿½mica (INTEC) Consejo Nacional de Investigaciones Cientï¿½ficas y Tï¿½cnicas (CONICET) PTLC - Gï¿½emes 3450, (3000) Santa Fe, Argentina Tel/Fax: +54-(0)342-451.1594

1 0

Re: Modified solver routines cr.c, bcgs.c and bcgsl.c
by Barry Smith 22 Sep '08

22 Sep '08

Matt, See my email to Satish about the bcgsl.c Barry On Sep 22, 2008, at 11:03 AM, Matthew Knepley wrote: > On Mon, Sep 22, 2008 at 11:01 AM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: >> >> Stephen, >> >> Thank you. We'll update PETSc with some of the changes. > > I am doing it. However, there is a strange bug in the cr.c you sent. > If you look at > line 111, you are using ai instead of -ai. I went back through the > log, and Barry > changed this, but also removed the line > > ai = -ai; > > beforehand. This is missing in your version, even though you have > not made the > sign change. Not sure how this could have run. > > Matt > >> Barry >> >> On Sep 22, 2008, at 9:54 AM, Stephen R Ball wrote: >> >>> >>> Hi >>> >>> I thought it best to send this to petsc-dev rather than petsc-users. >>> >>> One of our PETSc users has noticed that they obtain better >>> performance >>> when they use the modified versions of solver routines cr.c, >>> bcgs.c and >>> bcgsl.c attached. >>> In places they are using routines VecDotBegin()/VecDotEnd() >>> instead of >>> VecDot() and VecNormBegin()/VecNormEnd() instead of VecNorm(). >>> >>> I am sending them to you so that you may test them with a view to >>> you >>> possibly implementing them within PETSc if you find they improve >>> performance. >>> >>> Regards >>> >>> Stephen R. Ball >>> Advanced Technologies >>> HPC >>> DRAS >>> Rm: G17 >>> Bldg: E1.1 >>> AWE(A) >>> Aldermaston >>> Reading >>> Berkshire >>> ENGLAND >>> RG7 4PR >>> Tel: +44 (0)118 982 4528 >>> e-mail: stephen.r.ball(a)awe.co.uk >>> >>> >>> <<cr.c>> <<bcgs.c>> <<bcgsl.c>> >>> <cr.c><bcgs.c><bcgsl.c> >> >> > > > > -- > 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 >

1 0

Re: Modified solver routines cr.c, bcgs.c and bcgsl.c
by Satish Balay 22 Sep '08

22 Sep '08

On Mon, 22 Sep 2008, Matthew Knepley wrote: > On Mon, Sep 22, 2008 at 11:01 AM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > > > Stephen, > > > > Thank you. We'll update PETSc with some of the changes. > > I am doing it. However, there is a strange bug in the cr.c you sent. > If you look at > line 111, you are using ai instead of -ai. I went back through the > log, and Barry > changed this, but also removed the line > > ai = -ai; > > beforehand. This is missing in your version, even though you have not made the > sign change. Not sure how this could have run. Matt, You should do a diff with petsc-2.3.3 to see the actual changes.. Satish

1 0

Re: Modified solver routines cr.c, bcgs.c and bcgsl.c
by Matthew Knepley 22 Sep '08

22 Sep '08

On Mon, Sep 22, 2008 at 11:01 AM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Stephen, > > Thank you. We'll update PETSc with some of the changes. I am doing it. However, there is a strange bug in the cr.c you sent. If you look at line 111, you are using ai instead of -ai. I went back through the log, and Barry changed this, but also removed the line ai = -ai; beforehand. This is missing in your version, even though you have not made the sign change. Not sure how this could have run. Matt > Barry > > On Sep 22, 2008, at 9:54 AM, Stephen R Ball wrote: > >> >> Hi >> >> I thought it best to send this to petsc-dev rather than petsc-users. >> >> One of our PETSc users has noticed that they obtain better performance >> when they use the modified versions of solver routines cr.c, bcgs.c and >> bcgsl.c attached. >> In places they are using routines VecDotBegin()/VecDotEnd() instead of >> VecDot() and VecNormBegin()/VecNormEnd() instead of VecNorm(). >> >> I am sending them to you so that you may test them with a view to you >> possibly implementing them within PETSc if you find they improve >> performance. >> >> Regards >> >> Stephen R. Ball >> Advanced Technologies >> HPC >> DRAS >> Rm: G17 >> Bldg: E1.1 >> AWE(A) >> Aldermaston >> Reading >> Berkshire >> ENGLAND >> RG7 4PR >> Tel: +44 (0)118 982 4528 >> e-mail: stephen.r.ball(a)awe.co.uk >> >> >> <<cr.c>> <<bcgs.c>> <<bcgsl.c>> >> <cr.c><bcgs.c><bcgsl.c> > > -- 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

1 0

Re: Modified solver routines cr.c, bcgs.c and bcgsl.c
by Barry Smith 22 Sep '08

22 Sep '08

Stephen, Thank you. We'll update PETSc with some of the changes. Barry On Sep 22, 2008, at 9:54 AM, Stephen R Ball wrote: > > Hi > > I thought it best to send this to petsc-dev rather than petsc-users. > > One of our PETSc users has noticed that they obtain better performance > when they use the modified versions of solver routines cr.c, bcgs.c > and > bcgsl.c attached. > In places they are using routines VecDotBegin()/VecDotEnd() instead of > VecDot() and VecNormBegin()/VecNormEnd() instead of VecNorm(). > > I am sending them to you so that you may test them with a view to you > possibly implementing them within PETSc if you find they improve > performance. > > Regards > > Stephen R. Ball > Advanced Technologies > HPC > DRAS > Rm: G17 > Bldg: E1.1 > AWE(A) > Aldermaston > Reading > Berkshire > ENGLAND > RG7 4PR > Tel: +44 (0)118 982 4528 > e-mail: stephen.r.ball(a)awe.co.uk > > > <<cr.c>> <<bcgs.c>> <<bcgsl.c>> > <cr.c><bcgs.c><bcgsl.c>

1 0

Modified solver routines cr.c, bcgs.c and bcgsl.c
by Stephen R Ball 22 Sep '08

22 Sep '08

Hi I thought it best to send this to petsc-dev rather than petsc-users. One of our PETSc users has noticed that they obtain better performance when they use the modified versions of solver routines cr.c, bcgs.c and bcgsl.c attached. In places they are using routines VecDotBegin()/VecDotEnd() instead of VecDot() and VecNormBegin()/VecNormEnd() instead of VecNorm(). I am sending them to you so that you may test them with a view to you possibly implementing them within PETSc if you find they improve performance. Regards Stephen R. Ball Advanced Technologies HPC DRAS Rm: G17 Bldg: E1.1 AWE(A) Aldermaston Reading Berkshire ENGLAND RG7 4PR Tel: +44 (0)118 982 4528 e-mail: stephen.r.ball(a)awe.co.uk <<cr.c>> <<bcgs.c>> <<bcgsl.c>>

1 0

Re: PETSc on GPUs
by Harald Pfeiffer 19 Sep '08

19 Sep '08

Hi All, thanks for the many responses - people are certainly excited about GPUs and their possibilities. Our group at Caltech & Cornell uses PETSc to compute initial data for collisions between two black holes. The black hole evolution code is currently fully explicit, but we are working on implicit time-stepping (which heavily uses PETSc). The trick Barry mentions below (isolate heavy calculations into a few computational kernels; rewrite those kernels for GPUs) will apply to our code, too: There is one "data-class" which holds essentially all data, and it could easily be changed to keep the data in GPU memory. All operations on the data are performed via expression templates, which end up in only a handful routines with explicit loops over grid-points. Making these few routines GPU aware would then speed up all calculations (or so I hope). My original email was driven by curiosity. But when this gets more serious for us, I'll certainly write back with an update. Cheers, Harald On 09/18/2008 02:09 PM, Barry Smith wrote: > > In implementing Vec in the GPU a VecScatter would take the entries > directly from GPU copy of the > vector and put them appropriately into the GPU copy of the other vector, > thus it would be completely > transparent from the MatMult_MPIAIJ() code. In other words, if you have > to modify the MatMult_MPIAIJ > IN ANYWAY for the GPU then there is something wrong. > > Barry > > Note the VecScatter would need to marshal the values up from the GPU > into the regular memory, > do the appropriate message passing to the resulting processes and then > marshal the values back > down to the destination GPU. Since the current VecScatter handles all > the message passing I think > the GPU VecScatter() could be written to 1) move the needed GPU values > up to regular memory > (inside work MPI Vecs) then 2) use a regular VecScatter on those work > vectors then 3) move the results > down from the regular memory to the GPU. > > In other words PETSc has been designed since 1991 to be able to reuse > almost all the infrastructure > and only need new low level kernels written to run on GPUs. :-) > > Fun stuff, sure beats working :-) > > > > On Sep 18, 2008, at 3:18 PM, Richard Tran Mills wrote: > >> Ahmed, >> >> If you take a look in >> >> src/mat/impls/aij/mpi/ >> >> at mpiaij.h and mpiaij.c, you can see that an MPIAIJ matrix is stored >> locally as two "sequential" matrics, A and B. A contains the >> "diagonal" portion of the matrix, and B contains the "off-diagonal" >> portion. See page 59 of the PETSc User manual for a quick >> illustration of how this partitioning works. Because vectors are >> partitioned row-wise among processors in the same manner as the >> matrix, each process already has the vector entries required to form >> the portion of the matrix-vector product ("mat-vec") involving the >> diagonal portion 'A' can be done with data that reside locally. The >> off-diagonal portion 'B', however, requires gathering vector entries >> that reside off-processor. If you look in mpiaij.c at the >> MatMult_MPIAIJ, you will see the following lines: >> >> ierr = >> VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); >> >> ierr = (*a->A->ops->mult)(a->A,xx,yy);CHKERRQ(ierr); >> ierr = >> VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); >> >> ierr = (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);CHKERRQ(ierr); >> >> The mat-vec involving a->A can be completed without any VecScatter, >> but the one involving a->B needs the VecScatter to complete before it >> can be done. (Note that the VecScatterBegin() occurs before the >> mat-vec for a->A so that there might be some overlap of communication >> and computation). >> >> Hopefully this helps elucidate what Matt meant, and didn't just >> confuse you. >> >> --Richard >> >> Ahmed El Zein wrote: >>> On Thu, 2008-09-18 at 10:20 -0500, Matthew Knepley wrote: >>>>> A question that I had regarding the PETSc code when I was thinking >>>> about >>>>> this was: >>>>> You have the SeqAIJ matrix type and the the MPIAIJ type built around >>>> it >>>>> (or that is what I understand from the code). So basically you >>>> implement >>>>> the SeqAIJ type for the GPU and you get the MPI type for free? >>>> Yes, that is true. However, note that in the MPI step, you will need a >>>> gather >>>> operation to get the second matrix multiply to work. >>> Matt, >>> Could you explain a bit more? >>> Ahmed >> >> -- >> Richard Tran Mills, Ph.D. | E-mail: rmills(a)climate.ornl.gov >> Computational Scientist | Phone: (865) 241-3198 >> Computational Earth Sciences Group | Fax: (865) 574-0405 >> Oak Ridge National Laboratory | http://climate.ornl.gov/~rmills >> -- Harald P. Pfeiffer harald(a)tapir.caltech.edu Theoretical Astrophysics Phone (626) 395-8413 Caltech 130-33 Fax (626) 796-5675 Pasadena, CA 91125, USA

1 0

Re: PETSc on GPUs
by Barry Smith 18 Sep '08

18 Sep '08

In implementing Vec in the GPU a VecScatter would take the entries directly from GPU copy of the vector and put them appropriately into the GPU copy of the other vector, thus it would be completely transparent from the MatMult_MPIAIJ() code. In other words, if you have to modify the MatMult_MPIAIJ IN ANYWAY for the GPU then there is something wrong. Barry Note the VecScatter would need to marshal the values up from the GPU into the regular memory, do the appropriate message passing to the resulting processes and then marshal the values back down to the destination GPU. Since the current VecScatter handles all the message passing I think the GPU VecScatter() could be written to 1) move the needed GPU values up to regular memory (inside work MPI Vecs) then 2) use a regular VecScatter on those work vectors then 3) move the results down from the regular memory to the GPU. In other words PETSc has been designed since 1991 to be able to reuse almost all the infrastructure and only need new low level kernels written to run on GPUs. :-) Fun stuff, sure beats working :-) On Sep 18, 2008, at 3:18 PM, Richard Tran Mills wrote: > Ahmed, > > If you take a look in > > src/mat/impls/aij/mpi/ > > at mpiaij.h and mpiaij.c, you can see that an MPIAIJ matrix is > stored locally as two "sequential" matrics, A and B. A contains the > "diagonal" portion of the matrix, and B contains the "off-diagonal" > portion. See page 59 of the PETSc User manual for a quick > illustration of how this partitioning works. Because vectors are > partitioned row-wise among processors in the same manner as the > matrix, each process already has the vector entries required to form > the portion of the matrix-vector product ("mat-vec") involving the > diagonal portion 'A' can be done with data that reside locally. The > off-diagonal portion 'B', however, requires gathering vector entries > that reside off-processor. If you look in mpiaij.c at the > MatMult_MPIAIJ, you will see the following lines: > > ierr = VecScatterBegin(a->Mvctx,xx,a- > >lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); > ierr = (*a->A->ops->mult)(a->A,xx,yy);CHKERRQ(ierr); > ierr = VecScatterEnd(a->Mvctx,xx,a- > >lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); > ierr = (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);CHKERRQ(ierr); > > The mat-vec involving a->A can be completed without any VecScatter, > but the one involving a->B needs the VecScatter to complete before > it can be done. (Note that the VecScatterBegin() occurs before the > mat-vec for a->A so that there might be some overlap of > communication and computation). > > Hopefully this helps elucidate what Matt meant, and didn't just > confuse you. > > --Richard > > Ahmed El Zein wrote: >> On Thu, 2008-09-18 at 10:20 -0500, Matthew Knepley wrote: >>>> A question that I had regarding the PETSc code when I was thinking >>> about >>>> this was: >>>> You have the SeqAIJ matrix type and the the MPIAIJ type built >>>> around >>> it >>>> (or that is what I understand from the code). So basically you >>> implement >>>> the SeqAIJ type for the GPU and you get the MPI type for free? >>> Yes, that is true. However, note that in the MPI step, you will >>> need a >>> gather >>> operation to get the second matrix multiply to work. >> Matt, >> Could you explain a bit more? >> Ahmed > > -- > Richard Tran Mills, Ph.D. | E-mail: rmills(a)climate.ornl.gov > Computational Scientist | Phone: (865) 241-3198 > Computational Earth Sciences Group | Fax: (865) 574-0405 > Oak Ridge National Laboratory | http://climate.ornl.gov/~rmills >

1 0

Re: PETSc on GPUs
by Richard Tran Mills 18 Sep '08

18 Sep '08

Ahmed, If you take a look in src/mat/impls/aij/mpi/ at mpiaij.h and mpiaij.c, you can see that an MPIAIJ matrix is stored locally as two "sequential" matrics, A and B. A contains the "diagonal" portion of the matrix, and B contains the "off-diagonal" portion. See page 59 of the PETSc User manual for a quick illustration of how this partitioning works. Because vectors are partitioned row-wise among processors in the same manner as the matrix, each process already has the vector entries required to form the portion of the matrix-vector product ("mat-vec") involving the diagonal portion 'A' can be done with data that reside locally. The off-diagonal portion 'B', however, requires gathering vector entries that reside off-processor. If you look in mpiaij.c at the MatMult_MPIAIJ, you will see the following lines: ierr = VecScatterBegin(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); ierr = (*a->A->ops->mult)(a->A,xx,yy);CHKERRQ(ierr); ierr = VecScatterEnd(a->Mvctx,xx,a->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); ierr = (*a->B->ops->multadd)(a->B,a->lvec,yy,yy);CHKERRQ(ierr); The mat-vec involving a->A can be completed without any VecScatter, but the one involving a->B needs the VecScatter to complete before it can be done. (Note that the VecScatterBegin() occurs before the mat-vec for a->A so that there might be some overlap of communication and computation). Hopefully this helps elucidate what Matt meant, and didn't just confuse you. --Richard Ahmed El Zein wrote: > On Thu, 2008-09-18 at 10:20 -0500, Matthew Knepley wrote: >>> A question that I had regarding the PETSc code when I was thinking >> about >>> this was: >>> You have the SeqAIJ matrix type and the the MPIAIJ type built around >> it >>> (or that is what I understand from the code). So basically you >> implement >>> the SeqAIJ type for the GPU and you get the MPI type for free? >> Yes, that is true. However, note that in the MPI step, you will need a >> gather >> operation to get the second matrix multiply to work. > Matt, > Could you explain a bit more? > > Ahmed -- Richard Tran Mills, Ph.D. | E-mail: rmills(a)climate.ornl.gov Computational Scientist | Phone: (865) 241-3198 Computational Earth Sciences Group | Fax: (865) 574-0405 Oak Ridge National Laboratory | http://climate.ornl.gov/~rmills

1 0

Re: PETSc on GPUs
by Matthew Knepley 18 Sep '08

18 Sep '08

On Thu, Sep 18, 2008 at 1:34 PM, Ahmed El Zein <ahmed(a)azein.com> wrote: > On Thu, 2008-09-18 at 10:20 -0500, Matthew Knepley wrote: >> > A question that I had regarding the PETSc code when I was thinking >> about >> > this was: >> > You have the SeqAIJ matrix type and the the MPIAIJ type built around >> it >> > (or that is what I understand from the code). So basically you >> implement >> > the SeqAIJ type for the GPU and you get the MPI type for free? >> >> Yes, that is true. However, note that in the MPI step, you will need a >> gather >> operation to get the second matrix multiply to work. > Matt, > Could you explain a bit more? The diagonal block performs a simple, serial AIJ MatMult(). However, for the offdiagonal block, you must gather the Vec entries for nonzero columns in that block, fit them into a small Vec, do another serial AIJ MatMult(), and then scatter them back to the correct locations. Matt > Ahmed -- 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

1 0