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

Re: [petsc4py] What should be the 'default' communicator?
by Lisandro Dalcin 25 Aug '08

25 Aug '08

On Mon, Aug 25, 2008 at 1:08 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > I would still maintain that PETSC_COMM_WORLD is the correct default. There > are better paradigms for embarassingly parallel operation, like Condor. PETSc > is intended for parallel, domain decomposition runs. Yes, you are completelly right. But I believe that many people still use PETSc in a sequential way just because PETSc is full featured, well designed, easy to learn, etc. So, despite PETSc being intended for parallel, domain decomposition applications, many people are going to use it for sequential apps and embarassingly parallel operations. To be honest, I've never looked too much at paradigms like Condor. But using them implies to learn yet another framework. Another anecdote: a guy sent me a mail with questions about mpi4py for solving a embarassingly parallel problems. I asked why he was trying to use such a "heavy weight" approach. And then he answered he was tired of the complications and performance of using a Grid-based approach, and that 'mpiexec' a Python script with some coordinating MPI calls was far easier to setup, extend, and maintain and had better overall running times than submitting jobs to "The Grid". > On Mon, Aug 25, 2008 at 10:54 AM, Lisandro Dalcin <dalcinl(a)gmail.com> wrote: >> After working hard on mpi4py, this week I'll spend my time cleaning-up >> and adding features to the new Cython-based petsc4py. Then, I'll be >> asking questions to this list requesting for advise. >> >> In all calls that create new PETSc objects, I've decided to make the >> 'comm' argument optional. If the communicator is not passed, >> PETSC_COMM_WORLD is currently used. This is the approach PETSc uses in >> some C++ calls implemented through PetscPolymorphicFunction(). >> >> But now I believe that is wrong, and that PETSC_COMM_SELF should be >> the default. Or perhaps even better, I should let users set the >> default communicator used by petsc4py to create new (parallel) >> objects. >> >> An anecdote: some time ago, a petsc4py user wrote a sequential code >> and created objects without passing communicator arguments, next he >> wanted to solve many of those problems in different worker processes >> in a "ambarrasingly parallel" fashion and collect results at the >> master process. Of course, he run into trouble. Then I asked him to >> initialize PETSc in such a way that PETSC_COMM_WORLD was actually >> PETSC_COMM_SELF (by setting the world comm before PetscInitalize()). >> This mostly works, but has a problem: we have lost the actual >> PETSC_COMM_WORLD, so we are not able to create a parallel object after >> PetscInitialize(). >> >> Any thoughts? >> >> >> -- >> 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 >> >> > > > > -- > 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 > > -- 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: [petsc4py] What should be the 'default' communicator?
by Matthew Knepley 25 Aug '08

25 Aug '08

I would still maintain that PETSC_COMM_WORLD is the correct default. There are better paradigms for embarassingly parallel operation, like Condor. PETSc is intended for parallel, domain decomposition runs. Matt On Mon, Aug 25, 2008 at 10:54 AM, Lisandro Dalcin <dalcinl(a)gmail.com> wrote: > After working hard on mpi4py, this week I'll spend my time cleaning-up > and adding features to the new Cython-based petsc4py. Then, I'll be > asking questions to this list requesting for advise. > > In all calls that create new PETSc objects, I've decided to make the > 'comm' argument optional. If the communicator is not passed, > PETSC_COMM_WORLD is currently used. This is the approach PETSc uses in > some C++ calls implemented through PetscPolymorphicFunction(). > > But now I believe that is wrong, and that PETSC_COMM_SELF should be > the default. Or perhaps even better, I should let users set the > default communicator used by petsc4py to create new (parallel) > objects. > > An anecdote: some time ago, a petsc4py user wrote a sequential code > and created objects without passing communicator arguments, next he > wanted to solve many of those problems in different worker processes > in a "ambarrasingly parallel" fashion and collect results at the > master process. Of course, he run into trouble. Then I asked him to > initialize PETSc in such a way that PETSC_COMM_WORLD was actually > PETSC_COMM_SELF (by setting the world comm before PetscInitalize()). > This mostly works, but has a problem: we have lost the actual > PETSC_COMM_WORLD, so we are not able to create a parallel object after > PetscInitialize(). > > Any thoughts? > > > -- > 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 > > -- 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: [petsc4py] What should be the 'default' communicator?
by Barry Smith 25 Aug '08

25 Aug '08

What about a PETSC_COMM_DEFAULT that lives in the same place as PETSC_COMM_WORLD and PETSC_COMM_SELF.? You could have it default to PETSC_COMM_whatever and the user could change it immediately after PetscInitialize() to whatever they chose. For simplicity I wouldn't have a setter like PetscSetDefaultComm() rather I would just have the user set it directly, e.g. PETSC_COMM_DEFAULT = mynewcomm Barry If we ever did a proper C++ wrapper, we could use the same idea On Aug 25, 2008, at 10:54 AM, Lisandro Dalcin wrote: > After working hard on mpi4py, this week I'll spend my time cleaning-up > and adding features to the new Cython-based petsc4py. Then, I'll be > asking questions to this list requesting for advise. > > In all calls that create new PETSc objects, I've decided to make the > 'comm' argument optional. If the communicator is not passed, > PETSC_COMM_WORLD is currently used. This is the approach PETSc uses in > some C++ calls implemented through PetscPolymorphicFunction(). > > But now I believe that is wrong, and that PETSC_COMM_SELF should be > the default. Or perhaps even better, I should let users set the > default communicator used by petsc4py to create new (parallel) > objects. > > An anecdote: some time ago, a petsc4py user wrote a sequential code > and created objects without passing communicator arguments, next he > wanted to solve many of those problems in different worker processes > in a "ambarrasingly parallel" fashion and collect results at the > master process. Of course, he run into trouble. Then I asked him to > initialize PETSc in such a way that PETSC_COMM_WORLD was actually > PETSC_COMM_SELF (by setting the world comm before PetscInitalize()). > This mostly works, but has a problem: we have lost the actual > PETSC_COMM_WORLD, so we are not able to create a parallel object after > PetscInitialize(). > > Any thoughts? > > > -- > 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

[petsc4py] What should be the 'default' communicator?
by Lisandro Dalcin 25 Aug '08

25 Aug '08

After working hard on mpi4py, this week I'll spend my time cleaning-up and adding features to the new Cython-based petsc4py. Then, I'll be asking questions to this list requesting for advise. In all calls that create new PETSc objects, I've decided to make the 'comm' argument optional. If the communicator is not passed, PETSC_COMM_WORLD is currently used. This is the approach PETSc uses in some C++ calls implemented through PetscPolymorphicFunction(). But now I believe that is wrong, and that PETSC_COMM_SELF should be the default. Or perhaps even better, I should let users set the default communicator used by petsc4py to create new (parallel) objects. An anecdote: some time ago, a petsc4py user wrote a sequential code and created objects without passing communicator arguments, next he wanted to solve many of those problems in different worker processes in a "ambarrasingly parallel" fashion and collect results at the master process. Of course, he run into trouble. Then I asked him to initialize PETSc in such a way that PETSC_COMM_WORLD was actually PETSC_COMM_SELF (by setting the world comm before PetscInitalize()). This mostly works, but has a problem: we have lost the actual PETSC_COMM_WORLD, so we are not able to create a parallel object after PetscInitialize(). Any thoughts? -- 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: Adding PetscMap getters
by Barry Smith 22 Aug '08

22 Aug '08

Dave, There is a MatGetOwnershipRangeColumn() that might suit your needs. I have also just pushed the missing MatGetOwnershipRangesColumn() Barry On Aug 22, 2008, at 9:04 PM, Dave May wrote: > Hi Barry, > Thanks for your detailed response. It's refreshing to hear that > code in PETSc is introduced for the right reasons. I think the reasons > why I like PETSc so much is the result of the two strict rules you > adhere to. In my experience, many developers don't have such a strong > set of guide lines (or any guide lines at all) in mind when they write > code, and hence often introduce things for all the wrong reasons. > > The actual task I wished to perform was this; > Given a rectangular matrix and an (i,j) location, determine whether > (i,j) was in the matrix's diagonal block or off diagonal block. > > For square matrices this is easy as I can call MatGetOwnershipRanges() > and use the array for the row and column layout. But for rectangular > matrices I needed the column map. I could call MatGetVecs(), generate > a consistent left and right vector but then I still need to figure out > the ownership ranges for the vectors. (Plus I don't think generating > vectors just to determine the parallel layout of the matrix is very > smart.) > > One could say that I should just specify the parallel layout when I > built the matrix and vectors in the first place, then this problem > goes away. But this defeats the purpose of letting PETSc decide it > for me. I should still be able to query the parallel layout of the > matrix independent of whether I asked PETSc to determine the layout or > I prescribed it myself. > > I thought access to the maps would be the simplest approach, ie. > generate the least amount of work for you, as only a very small > getter would need to be written. Given the guide lines you follow > (which I wasn't aware of) I see why this isn't the right thing to do. > > So what about this. > Forget about accessing the PetscMap object directly (meeting the first > rule of PETSc programming) and instead add routines to fully determine > the parallel layout of the vectors and matrices. > > For vectors this would mean introducing > VecGetOwnershipRanges() > > For matrices this would mean introducing > MatGetRowOwnershipRange() ( which is identical to > MatGetOwnershipRange() ) > MatGetRowOwnershipRanges() ( which is identical to > MatGetOwnershipRanges() ) > > MatGetColOwnershipRange() > MatGetColOwnershipRanges() > > Unless I'm mistaken (please correct me if I'm wrong), there is > currently no simple and direct way to determine the parallel layout of > a matrix in terms of both of row and column ranges. Hence adding these > interfaces for the matrices would not violate of the second rule of > PETSc programming. > > Cheers, > Dave > > > > On Sat, Aug 23, 2008 at 1:14 AM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: >> >> Dave, >> >> This has been a controversial issue over the years. We didn't >> have a map >> concept, it was introduced, it was then >> "hidden". The fact that there exist public PetscMap functions at >> all is my >> compromise for a persistent user. >> >> Two key "rules" I use in making decisions in PETSc code are: >> >> 1) try really really hard to minimize the number of abstract >> objects in the >> PETSc interface (this is because I believe >> the difficulty of using a "class" library grows exponentially >> with the >> number of abstract objects in the interface) >> (obviously not everyone believes this, it seems most class >> libraries add >> new classes at the drop of a hat) >> 2) there should be ONE way to do each operation in the package >> (this is >> because I believe having multiple ways to >> do the same thing is confusing to most users). For example, to set >> the >> size of a vector you can only use VecSetSizes(), >> to solve a nonlinear system you call SNESSolve()*. >> >> If we bring PetscMap up as "first class user object" this violates >> both 1 >> and 2. It is a new basic object and manipulating/accessing >> Vec/Mat size information would have TWO ways of doing everything: >> VecGetSizes(), VecGetPetscMap() -> PetscMapGetSizes() etc. >> I actually dream about the opposite; shrinking PetscMap down so it >> is so >> small it can be strangled in a bathtub. I tried once >> to remove the VecGetSizes() etc and have access only through >> PetscMap (to >> satisfy rule 2) but ended up removing the public >> PetscMap because I didn't see having PetscMap as the basic "layout >> definer" >> as a good PETSc interface. >> >> So my answer is, I really do not want to add these getters. you >> really >> have to convince me that the PetscMap object >> is so fundamentally powerful it is worth violating 1 and 2. >> >> Barry >> >> >> Note that the PetscMap() methods are prototyped in vecimpl.h, >> indicating >> they are not public methods for PETSc. >> >> >> >> >> >> *sadly there is also DMMGSolve(), this I view as an interface bug >> and dream >> about someday merging this functionality >> back to the SNES. >> >> >> On Aug 21, 2008, at 5:52 PM, Dave May wrote: >> >>> Hello, >>> Would it be possible to have the following small getters added to >>> petsc-dev? >>> VecGetOwnershipRanges() >>> Since we are just manipulating the maps why not include the >>> following as >>> well >>> VecGetMap() >>> MatGetRowMap() >>> MatGetColumnMap() >>> >>> It seems reasonable to add these (to me at least) given there a few >>> interfaces for >>> PetscMap already made public. >>> >>> Cheers, >>> Dave. >>> >> >> >

1 0

Re: Adding PetscMap getters
by Dave May 22 Aug '08

22 Aug '08

Hi Barry, Thanks for your detailed response. It's refreshing to hear that code in PETSc is introduced for the right reasons. I think the reasons why I like PETSc so much is the result of the two strict rules you adhere to. In my experience, many developers don't have such a strong set of guide lines (or any guide lines at all) in mind when they write code, and hence often introduce things for all the wrong reasons. The actual task I wished to perform was this; Given a rectangular matrix and an (i,j) location, determine whether (i,j) was in the matrix's diagonal block or off diagonal block. For square matrices this is easy as I can call MatGetOwnershipRanges() and use the array for the row and column layout. But for rectangular matrices I needed the column map. I could call MatGetVecs(), generate a consistent left and right vector but then I still need to figure out the ownership ranges for the vectors. (Plus I don't think generating vectors just to determine the parallel layout of the matrix is very smart.) One could say that I should just specify the parallel layout when I built the matrix and vectors in the first place, then this problem goes away. But this defeats the purpose of letting PETSc decide it for me. I should still be able to query the parallel layout of the matrix independent of whether I asked PETSc to determine the layout or I prescribed it myself. I thought access to the maps would be the simplest approach, ie. generate the least amount of work for you, as only a very small getter would need to be written. Given the guide lines you follow (which I wasn't aware of) I see why this isn't the right thing to do. So what about this. Forget about accessing the PetscMap object directly (meeting the first rule of PETSc programming) and instead add routines to fully determine the parallel layout of the vectors and matrices. For vectors this would mean introducing VecGetOwnershipRanges() For matrices this would mean introducing MatGetRowOwnershipRange() ( which is identical to MatGetOwnershipRange() ) MatGetRowOwnershipRanges() ( which is identical to MatGetOwnershipRanges() ) MatGetColOwnershipRange() MatGetColOwnershipRanges() Unless I'm mistaken (please correct me if I'm wrong), there is currently no simple and direct way to determine the parallel layout of a matrix in terms of both of row and column ranges. Hence adding these interfaces for the matrices would not violate of the second rule of PETSc programming. Cheers, Dave On Sat, Aug 23, 2008 at 1:14 AM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Dave, > > This has been a controversial issue over the years. We didn't have a map > concept, it was introduced, it was then > "hidden". The fact that there exist public PetscMap functions at all is my > compromise for a persistent user. > > Two key "rules" I use in making decisions in PETSc code are: > > 1) try really really hard to minimize the number of abstract objects in the > PETSc interface (this is because I believe > the difficulty of using a "class" library grows exponentially with the > number of abstract objects in the interface) > (obviously not everyone believes this, it seems most class libraries add > new classes at the drop of a hat) > 2) there should be ONE way to do each operation in the package (this is > because I believe having multiple ways to > do the same thing is confusing to most users). For example, to set the > size of a vector you can only use VecSetSizes(), > to solve a nonlinear system you call SNESSolve()*. > > If we bring PetscMap up as "first class user object" this violates both 1 > and 2. It is a new basic object and manipulating/accessing > Vec/Mat size information would have TWO ways of doing everything: > VecGetSizes(), VecGetPetscMap() -> PetscMapGetSizes() etc. > I actually dream about the opposite; shrinking PetscMap down so it is so > small it can be strangled in a bathtub. I tried once > to remove the VecGetSizes() etc and have access only through PetscMap (to > satisfy rule 2) but ended up removing the public > PetscMap because I didn't see having PetscMap as the basic "layout definer" > as a good PETSc interface. > > So my answer is, I really do not want to add these getters. you really > have to convince me that the PetscMap object > is so fundamentally powerful it is worth violating 1 and 2. > > Barry > > > Note that the PetscMap() methods are prototyped in vecimpl.h, indicating > they are not public methods for PETSc. > > > > > > *sadly there is also DMMGSolve(), this I view as an interface bug and dream > about someday merging this functionality > back to the SNES. > > > On Aug 21, 2008, at 5:52 PM, Dave May wrote: > >> Hello, >> Would it be possible to have the following small getters added to >> petsc-dev? >> VecGetOwnershipRanges() >> Since we are just manipulating the maps why not include the following as >> well >> VecGetMap() >> MatGetRowMap() >> MatGetColumnMap() >> >> It seems reasonable to add these (to me at least) given there a few >> interfaces for >> PetscMap already made public. >> >> Cheers, >> Dave. >> > >

1 0

Re: Adding PetscMap getters
by Barry Smith 22 Aug '08

22 Aug '08

Dave, This has been a controversial issue over the years. We didn't have a map concept, it was introduced, it was then "hidden". The fact that there exist public PetscMap functions at all is my compromise for a persistent user. Two key "rules" I use in making decisions in PETSc code are: 1) try really really hard to minimize the number of abstract objects in the PETSc interface (this is because I believe the difficulty of using a "class" library grows exponentially with the number of abstract objects in the interface) (obviously not everyone believes this, it seems most class libraries add new classes at the drop of a hat) 2) there should be ONE way to do each operation in the package (this is because I believe having multiple ways to do the same thing is confusing to most users). For example, to set the size of a vector you can only use VecSetSizes(), to solve a nonlinear system you call SNESSolve()*. If we bring PetscMap up as "first class user object" this violates both 1 and 2. It is a new basic object and manipulating/accessing Vec/Mat size information would have TWO ways of doing everything: VecGetSizes(), VecGetPetscMap() -> PetscMapGetSizes() etc. I actually dream about the opposite; shrinking PetscMap down so it is so small it can be strangled in a bathtub. I tried once to remove the VecGetSizes() etc and have access only through PetscMap (to satisfy rule 2) but ended up removing the public PetscMap because I didn't see having PetscMap as the basic "layout definer" as a good PETSc interface. So my answer is, I really do not want to add these getters. you really have to convince me that the PetscMap object is so fundamentally powerful it is worth violating 1 and 2. Barry Note that the PetscMap() methods are prototyped in vecimpl.h, indicating they are not public methods for PETSc. *sadly there is also DMMGSolve(), this I view as an interface bug and dream about someday merging this functionality back to the SNES. On Aug 21, 2008, at 5:52 PM, Dave May wrote: > Hello, > Would it be possible to have the following small getters added to > petsc-dev? > VecGetOwnershipRanges() > Since we are just manipulating the maps why not include the > following as well > VecGetMap() > MatGetRowMap() > MatGetColumnMap() > > It seems reasonable to add these (to me at least) given there a few > interfaces for > PetscMap already made public. > > Cheers, > Dave. >

1 0

Adding PetscMap getters
by Dave May 21 Aug '08

21 Aug '08

Hello, Would it be possible to have the following small getters added to petsc-dev? VecGetOwnershipRanges() Since we are just manipulating the maps why not include the following as well VecGetMap() MatGetRowMap() MatGetColumnMap() It seems reasonable to add these (to me at least) given there a few interfaces for PetscMap already made public. Cheers, Dave.

1 0

Re: New approach for external direct solvers in PETSc
by Barry Smith 18 Aug '08

18 Aug '08

Keita, I don't see anything but very minor fixes to the current config/ configure.py and install process. Barry On Aug 18, 2008, at 1:57 PM, Keita Teranishi wrote: > Hi Barry, > > I'd like to grab the dev version with new direct solver interface. > Just a question. Will the packaging (bmake and configure script) > form of the final release be similar to the dev version? > > Thanks, > > ================================ > Keita Teranishi > Math Software Group > Cray, Inc. > keita(a)cray.com > ================================ > > > -----Original Message----- > From: owner-petsc-dev(a)mcs.anl.gov [mailto:owner-petsc- > dev(a)mcs.anl.gov] On Behalf Of Barry Smith > Sent: Sunday, August 17, 2008 9:35 AM > To: petsc-dev(a)mcs.anl.gov > Subject: Re: New approach for external direct solvers in PETSc > > > We will most definitely make a release by December 1 2008. This > means any > changes like this will have to be made pretty soon. Any changes will > actually > have only a small effect on the user interface (but a bigger effect on > PETSc's simplicity). > > Barry > > On Aug 16, 2008, at 11:55 PM, Chen Shen wrote: > >> Hi, >> >> May I know if this change will appear in petsc-devel soon? We are >> trying to interface Petsc to a proprietary direct solver package. If >> we should expect the change soon, we may consider push back the plan >> a little bit. Thank you. >> >> regards, >> shenchen >> >> On Tue, Aug 12, 2008 at 12:40 AM, Barry Smith <bsmith(a)mcs.anl.gov> >> wrote: >> >> Mat in PETSc is actually a short hand notation for "Linear >> transformation" >> http://en.wikipedia.org/wiki/Linear_transformation. Essentially the >> only >> operation defined by a linear transformation is A*x. Matrices are >> simply one >> way of representing finite dimensional linear transformations. So a >> Mat is suppose to ALWAYS provide a MatMult() everything else is >> optional :-). >> >> Now you can certainly say that we should have called linear >> transformations >> in PETSc something like: LT or LinearTrans or many other things. We >> chose >> Mat since almost everyone understands linear transformations as >> matrices so >> we do not need to provide a high-falutin explanation of LT that >> would simply >> alienate and scare most potential users. Everyone thinks they >> understand >> Mat immediately. >> >> >> Now back to the issue at hand. You had a very valid point in the >> previous >> email: for a factored matrix if we define MatMult() as the >> triangular solves >> then the factored matrix defines a linear transformation (whose >> representation >> is not simply a two dimensional array, but is slightly more >> complicated beast). >> So we could define a MATFACTOR as a new Mat object, like MatMFFD is >> a Mat object and then subclass the various packages factors >> undernether >> it. >> >> Barry >> >> >> >> >> >> >> >> On Aug 11, 2008, at 9:11 AM, Lisandro Dalcin wrote: >> >> On Fri, Aug 8, 2008 at 2:01 PM, Barry Smith <bsmith(a)mcs.anl.gov> >> wrote: >> This is just at hack to avoid so many changes in PETSc. >> MatMFFD is truly a matrix (it has matrix vector product). >> >> Well, I'm not so sure about that. Do it has MatSetValues() ? >> >> >> >> Barry >> >> On Aug 8, 2008, at 8:23 AM, Lisandro Dalcin wrote: >> >> Barry, after thinking a bit more, I believe I've found an approach >> that is not so radical. >> >> Why not to follow the approach for MatMFFD? That is, introduce a new >> matrix type, instead of a brand new object type. Let me call this new >> matrix type MATFACTOR. This new matrix type contains a hidden >> PetscObject, and that object is filled with appropriate methods >> depending on the MatSolverPackage. >> >> This new way will not introduce a new object type in the user-level >> API. An I believe that we will still simplify the code. Now a >> MatSolve >> in a MATSEQAIJ matrix will fail just because it will not have the >> operation filled in the 'virtual table'. MATFACTOR will fill the >> MatXXXFactor{Symbolic/Numeric} options in the vtable, but it will not >> fill things like MatMult and MatSetValues. This way, we will not need >> any more to 'check' if a standard matrix type like SEQAIJ, MPIAIJ, >> etc >> are or are not factored (and perhaps we can handle dense matrices >> in a >> special way, just because inplace factorizations do really make >> sense). So the all the 'mat->factor' checks can finally go away! >> >> What do you think? >> >> >> >> >> On Thu, Aug 7, 2008 at 5:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> >> wrote: >> >> On Jul 25, 2008, at 5:47 PM, Lisandro Dalcin wrote: >> >> BTW, Have you ever considered the posibility of not using the 'Mat' >> class for factored matrices?? I believe that normal matrices and >> factored matrices share very few in common as to being instances of >> the same class. Looking at the source code in src/mat, there are many >> checks like 'if (mat-->factor)' or 'if (!mat->factor)' that seems to >> support my concern. What do you think? >> >> After a (tiny) bit more thought I think we really should make this >> change: >> >> pros >> 1) it will simplify the code a good amount hence less bugs >> 2) conceptually it makes sense >> cons >> 1) requires a large reorganization of the Mat code (will introduce >> bugs) >> 2) requires the user to deal with one more PETSc object (MatFactor) >> But since most users deal with SNES/KSP/PC most users will >> never see this new object. >> >> The question is when/how to make this massive change. It is as big >> or bigger than the change I already made with MatGetFactor(). >> >> Barry >> >> >> >> >> >> On Fri, Jul 25, 2008 at 7:02 PM, Barry Smith <bsmith(a)mcs.anl.gov> >> wrote: >> >> petsc-dev users, >> >> I have finally completed my changes to PETSc for a new approach to >> accessing external direct solvers >> in PETSc like Spooles, MUMPS etc. I will be pushing it to petsc-dev >> sometime >> the middle of next week. >> If you are using the direct solvers you might consider cloning from >> >> http://petsc.cs.iit.edu/hg/petsc/petsc-dev-new-solvers or >> >> ssh://[email protected]//hg/petsc/petsc-dev-new-solvers >> >> and trying it out before then. Please report problem to >> petsc-maint(a)mcs.anl.gov as you hit them. >> >> >> Barry >> >> From the changes file >> >> The "parallel direct solver" matrix types like >> MATAIJSPOOLES are ALL gone. Now you use -pc_factor_mat_solver_package >> spooles etc or PCFactorSetMatSolverPackage() or if working directly >> with >> matrices, MatGetFactor(A,MATSPOOLES,...) >> >> >> >> >> >> -- >> 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 >> >> >> >> >> >> >> -- >> 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 >> >> >> >> >> >> >> -- >> 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: New approach for external direct solvers in PETSc
by Keita Teranishi 18 Aug '08

18 Aug '08

Hi Barry, I'd like to grab the dev version with new direct solver interface. Just a question. Will the packaging (bmake and configure script) form of the final release be similar to the dev version? Thanks, ================================ Keita Teranishi Math Software Group Cray, Inc. keita(a)cray.com ================================ -----Original Message----- From: owner-petsc-dev(a)mcs.anl.gov [mailto:[email protected]] On Behalf Of Barry Smith Sent: Sunday, August 17, 2008 9:35 AM To: petsc-dev(a)mcs.anl.gov Subject: Re: New approach for external direct solvers in PETSc We will most definitely make a release by December 1 2008. This means any changes like this will have to be made pretty soon. Any changes will actually have only a small effect on the user interface (but a bigger effect on PETSc's simplicity). Barry On Aug 16, 2008, at 11:55 PM, Chen Shen wrote: > Hi, > > May I know if this change will appear in petsc-devel soon? We are > trying to interface Petsc to a proprietary direct solver package. If > we should expect the change soon, we may consider push back the plan > a little bit. Thank you. > > regards, > shenchen > > On Tue, Aug 12, 2008 at 12:40 AM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > Mat in PETSc is actually a short hand notation for "Linear > transformation" > http://en.wikipedia.org/wiki/Linear_transformation. Essentially the > only > operation defined by a linear transformation is A*x. Matrices are > simply one > way of representing finite dimensional linear transformations. So a > Mat is suppose to ALWAYS provide a MatMult() everything else is > optional :-). > > Now you can certainly say that we should have called linear > transformations > in PETSc something like: LT or LinearTrans or many other things. We > chose > Mat since almost everyone understands linear transformations as > matrices so > we do not need to provide a high-falutin explanation of LT that > would simply > alienate and scare most potential users. Everyone thinks they > understand > Mat immediately. > > > Now back to the issue at hand. You had a very valid point in the > previous > email: for a factored matrix if we define MatMult() as the > triangular solves > then the factored matrix defines a linear transformation (whose > representation > is not simply a two dimensional array, but is slightly more > complicated beast). > So we could define a MATFACTOR as a new Mat object, like MatMFFD is > a Mat object and then subclass the various packages factors > undernether > it. > > Barry > > > > > > > > On Aug 11, 2008, at 9:11 AM, Lisandro Dalcin wrote: > > On Fri, Aug 8, 2008 at 2:01 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > This is just at hack to avoid so many changes in PETSc. > MatMFFD is truly a matrix (it has matrix vector product). > > Well, I'm not so sure about that. Do it has MatSetValues() ? > > > > Barry > > On Aug 8, 2008, at 8:23 AM, Lisandro Dalcin wrote: > > Barry, after thinking a bit more, I believe I've found an approach > that is not so radical. > > Why not to follow the approach for MatMFFD? That is, introduce a new > matrix type, instead of a brand new object type. Let me call this new > matrix type MATFACTOR. This new matrix type contains a hidden > PetscObject, and that object is filled with appropriate methods > depending on the MatSolverPackage. > > This new way will not introduce a new object type in the user-level > API. An I believe that we will still simplify the code. Now a MatSolve > in a MATSEQAIJ matrix will fail just because it will not have the > operation filled in the 'virtual table'. MATFACTOR will fill the > MatXXXFactor{Symbolic/Numeric} options in the vtable, but it will not > fill things like MatMult and MatSetValues. This way, we will not need > any more to 'check' if a standard matrix type like SEQAIJ, MPIAIJ, etc > are or are not factored (and perhaps we can handle dense matrices in a > special way, just because inplace factorizations do really make > sense). So the all the 'mat->factor' checks can finally go away! > > What do you think? > > > > > On Thu, Aug 7, 2008 at 5:05 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > On Jul 25, 2008, at 5:47 PM, Lisandro Dalcin wrote: > > BTW, Have you ever considered the posibility of not using the 'Mat' > class for factored matrices?? I believe that normal matrices and > factored matrices share very few in common as to being instances of > the same class. Looking at the source code in src/mat, there are many > checks like 'if (mat-->factor)' or 'if (!mat->factor)' that seems to > support my concern. What do you think? > > After a (tiny) bit more thought I think we really should make this > change: > > pros > 1) it will simplify the code a good amount hence less bugs > 2) conceptually it makes sense > cons > 1) requires a large reorganization of the Mat code (will introduce > bugs) > 2) requires the user to deal with one more PETSc object (MatFactor) > But since most users deal with SNES/KSP/PC most users will > never see this new object. > > The question is when/how to make this massive change. It is as big > or bigger than the change I already made with MatGetFactor(). > > Barry > > > > > > On Fri, Jul 25, 2008 at 7:02 PM, Barry Smith <bsmith(a)mcs.anl.gov> > wrote: > > petsc-dev users, > > I have finally completed my changes to PETSc for a new approach to > accessing external direct solvers > in PETSc like Spooles, MUMPS etc. I will be pushing it to petsc-dev > sometime > the middle of next week. > If you are using the direct solvers you might consider cloning from > > http://petsc.cs.iit.edu/hg/petsc/petsc-dev-new-solvers or > > ssh://[email protected]//hg/petsc/petsc-dev-new-solvers > > and trying it out before then. Please report problem to > petsc-maint(a)mcs.anl.gov as you hit them. > > > Barry > > From the changes file > > The "parallel direct solver" matrix types like > MATAIJSPOOLES are ALL gone. Now you use -pc_factor_mat_solver_package > spooles etc or PCFactorSetMatSolverPackage() or if working directly > with > matrices, MatGetFactor(A,MATSPOOLES,...) > > > > > > -- > 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 > > > > > > > -- > 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 > > > > > > > -- > 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