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

Re: using MUMPS with mg
by Barry Smith 23 Jun '08

23 Jun '08

Lisandro, The way I see it the problem with the current code is that any Mat that one creates "knows" how to do factorizations. For example it has MatCholeskyFactorSymbolic(). Thus if one wants to use a Mumps factorization one needs a different matrix class (MATAIJMUMPS for example), this is difficult for users to understand and also silly since WAY back in the code before I start filling up my matrix, the LU I want to use comes into play. The original "fix" I proposed was MatSetSolverType() that would allow one to set the solver type to the matrix LATER. This would have essentially been like calling a MatConvert() on the matrix to a specialized form like MATAIJMUMPS. A major drawback to this approach is that it keeps all the horrible, horrible, MatConvert() code to the specialized type that is ugly, bug ridden and impossible to maintain extend and understand. Then I realized a clean solution that eliminates ALL the ugly convert code. It comes with the realization that the factor routines belong as methods on the FACTORED matrix (since it is the factored matrix algorithm and data structures that are affected by the MUMPS, SuperLU etc), not the original matrix. Thus ALL the convert crap can be dumped :-). My current test model is MatGetFactor(Mat A,MatSolverType,MatFactorType,Mat *fact); followed by MatLUFactorSymbolic(fact,A,.....) Once I have it all working I am, of course, open to name suggestions and tweaks. Once you see it you'll like it :-) Barry On Jun 23, 2008, at 2:29 PM, Lisandro Dalcin wrote: > On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> >> It will be gone. >> > > Ups! Was this decided in some private discussion? I remember the post > were you proposed this; I also remember there were some objections, > but you started it anyway, so I guessed that this was already > accepted. > > IMHO, the whole idea was +1 for me (though I did not commented on > this). Did you change your mind? Was this rejected for the other > developers in the PETSc team? Was a better mechanism devised? > > > -- > 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: #define XXXType const char* ----> #define XXXType char* ?
by Barry Smith 23 Jun '08

23 Jun '08

> Mmm, just a bit ugly you have to add the 'const' keyword. I insist: > What would be the problem with ... ?? > > typedef const char* VecType; Two problems. The VecType is a const which means one cannot build them on the fly and copy into them etc Second with the typedef one cannot do all the string manipulation on them. For example even the PetscStrcpy() in PetscHeaderCreate_Private() of type_name does not work. It would be nice to use a typedef instead of a define to build it but it looks impossible. (Try it, just change one of them and try to compile PETSc, the screaming is relentless.) Barry On Jun 23, 2008, at 2:38 PM, Lisandro Dalcin wrote: > On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> >> A while ago I made XXType (MatType, VecType etc) const char* >> > > Well, you actually did > > #define VecType const char* > > instead of > > typedef const char* VecType; > > which IMHO is the right way. > >> I now think that was a mistake. > > Yep, the '#define' way can cause problems ;-) > > >> I'd like to change them to just char* > > But why not the 'typedef' way ?? > >> Can anyone suggest a technical reason why this is a bad idea >> before I push >> the code? >> >> Note this means if you have code like >> >> VecType t = VECMPI; >> >> you would need to change it to >> >> const VecType t = VECMPI; >> > > > > -- > 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: #define XXXType const char* ----> #define XXXType char* ?
by Lisandro Dalcin 23 Jun '08

23 Jun '08

On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > A while ago I made XXType (MatType, VecType etc) const char* > Well, you actually did #define VecType const char* instead of typedef const char* VecType; which IMHO is the right way. > I now think that was a mistake. Yep, the '#define' way can cause problems ;-) > I'd like to change them to just char* But why not the 'typedef' way ?? > Can anyone suggest a technical reason why this is a bad idea before I push > the code? > > Note this means if you have code like > > VecType t = VECMPI; > > you would need to change it to > > const VecType t = VECMPI; > Mmm, just a bit ugly you have to add the 'const' keyword. I insist: What would be the problem with ... ?? typedef const char* VecType; -- 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: using MUMPS with mg
by Lisandro Dalcin 23 Jun '08

23 Jun '08

On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > It will be gone. > Ups! Was this decided in some private discussion? I remember the post were you proposed this; I also remember there were some objections, but you started it anyway, so I guessed that this was already accepted. IMHO, the whole idea was +1 for me (though I did not commented on this). Did you change your mind? Was this rejected for the other developers in the PETSc team? Was a better mechanism devised? -- 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

#define XXXType const char* ----> #define XXXType char* ?
by Barry Smith 23 Jun '08

23 Jun '08

A while ago I made XXType (MatType, VecType etc) const char* I now think that was a mistake. I'd like to change them to just char* XXType and change the prototypes XXXSetType(XXXXX x,XXXType t) ------> XXXSetType(XXXXXX x,const XXXType t) XXXGetType(XXXXX x,XXXType *t) -----> XXXGetType(XXXXXX x,const XXTYPE* t) Can anyone suggest a technical reason why this is a bad idea before I push the code? I will warn before I push so that if you are in the middle of some delicate coding and are afraid this will break things you can avoid pulling for a while. Thanks Barry Note this means if you have code like VecType t = VECMPI; you would need to change it to const VecType t = VECMPI;

1 0

Re: using MUMPS with mg
by Barry Smith 23 Jun '08

23 Jun '08

On Jun 23, 2008, at 1:44 PM, Lisandro Dalcin wrote: > On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: >> >> Even better than this setsolvertype nonsense. >> > > Sorry, Barry... Perhaps I did not understand you. What is the role of > MatSetSolverType() in future PETSc releases ? It will be gone. Barry > > > > -- > 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: using MUMPS with mg
by Lisandro Dalcin 23 Jun '08

23 Jun '08

On 6/23/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Even better than this setsolvertype nonsense. > Sorry, Barry... Perhaps I did not understand you. What is the role of MatSetSolverType() in future PETSc releases ? -- 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: using MUMPS with mg
by Barry Smith 23 Jun '08

23 Jun '08

On Jun 23, 2008, at 1:16 PM, Lisandro Dalcin wrote: > On 6/23/08, Thomas Geenen <geenen(a)gmail.com> wrote: > >> i will check if just adding >> >> ierr = MatConvert(A, MATAIJMUMPS, MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); >> >> will do the trick >> > > Indeed. That would do the trick. But perhaps better is to add: > > MatConvert(A, MATSAME, MAT_REUSE_MATRIX,&A) > > Then, you can pass '-matconvert_type aijmumps' (check for the actual > option name, to busy right now to look at the source) to actually use > MATAIJMUMPS. If you do not pass the option, the the MatConvert() call > is just a non-op . > > Perhaps this is a candidate for petsc-dev? Or MatSetSolverType() will > handle this in the near future? Even better than this setsolvertype nonsense. Barry > > > > > -- > 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: using MUMPS with mg
by Lisandro Dalcin 23 Jun '08

23 Jun '08

On 6/23/08, Thomas Geenen <geenen(a)gmail.com> wrote: > i will check if just adding > > ierr = MatConvert(A, MATAIJMUMPS, MAT_REUSE_MATRIX,&A);CHKERRQ(ierr); > > will do the trick > Indeed. That would do the trick. But perhaps better is to add: MatConvert(A, MATSAME, MAT_REUSE_MATRIX,&A) Then, you can pass '-matconvert_type aijmumps' (check for the actual option name, to busy right now to look at the source) to actually use MATAIJMUMPS. If you do not pass the option, the the MatConvert() call is just a non-op . Perhaps this is a candidate for petsc-dev? Or MatSetSolverType() will handle this in the near future? -- 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: some fixes for AIJ and BAIJ matrices
by Lisandro Dalcin 11 Jun '08

11 Jun '08

On 6/10/08, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > Beats the heck out of me, do whatever fixes you feel are needed. Well, I've finally decided not push the fixes directly in petsc-dev. Core AIJ preallocation routines are not good place for making fixes for release, I'm really worried any change there could cause problems to users. I did not do anything about supporting MAT_SKIP_ALLOCATION for MPIAIJ matrices. After looking at the use cases, it seems that this is not intended for end users, but for implementing MPI matrices containing Seq matrices. If this is the case, then nothing remains to be done for this. Note that I've not managed the Fortran interfaces for Mat{Seq|MPI}BAIJSetPreallocationCSR. I'm not completelly sure if the automatic generation of Fortran wrappers is safe in this case (though it seems it is). Feel free to add it if any of you consider that appropriate. I've tested the fixes from petsc4py with many variants of sequential/parallel (2,3,4 procs), AIJ/BAIJ, grid sizes (including zero) and block sizes (1 to 5) (suming up, about 200 test cases in a Python unittest script) to believe that all is working OK. Anyway, I ran 'make testexamples' and got the following two diffs: testexamples_C in: /u/dalcinl/Devel/PETSc/petsc-dev/src/mat/examples/tests (1) The one below should perhaps need a bit more of attention. Looking at the actual code, I would say this does not have nothing to do with my fixes. 42c42 < nonzeros = 48, allocated nonzeros = 120 --- > nonzeros = 48, allocated nonzeros = 84 44c44 < nonzeros = 25, allocated nonzeros = 60 --- > nonzeros = 25, allocated nonzeros = 42 Possible problem with ex9_1, diffs above (2) The one below is surelly outdated. BDIAG matrices are gone. 26,80c26,48 < type=seqbdiag, rows=20, cols=20 ... <more ouput> ... --- > type=seqaij, rows=20, cols=20 ... <more ouput> ... Possible problem with ex20_1, diffs above Regards, > > Some fixes are needed for BAIJ matrices are: > > > > (1) MatMPIBAIJSetPreallocationCSR_MPIBAIJ -> BROKEN > > (2) MatSeqBAIJSetPreallocationCSR_SeqBAIJ -> MISSING > > > > Regarding (1), the fix needs to call PetscMapSetUp() for setting-up > > the row and column maps (rmap&cmap fields) at the very begining, > > before anything else. If this is not done, the stuff is the row/column > > maps are not properly initialized at the point (1) is called. But > > then, PetscMapSetUp() will be called again as (1) calls > > MatMPIBAIJSetPreallocation_MPIBAIJ with nnz information. > However, I > > believe this is not an issue, as the second call will be a no-op (at > > least for optimized builds). Additionally, there is a bug at the point > > MatSetValuesBlocked_MPIBAIJ() is being called, as the block-size is > > not being taking into accout for the 'values' array (when it is > > provided). > > > > Regarding (2), the issues in (1) for rmap/cmap are the same. > > > > It is fine to push the fixes in release-2.3.3, particularly regarding > > that (2) could be seen a 'new functionality'? I would prefer to add it > > to release-2.3.3. > > > > ------------------------- > > > > Finally, for MPIAIJ, something I'm not sure what to do is the following: > > > > - for MatSeq[B]AIJ one can pass nz=MAT_SKIP_PREALLOCATION > > - for MatMPI[B]AIJ, passing > d_nz/o_nz=MAT_SKIP_PREALLOCATION generates > > an error. > > > > Should MatMPI[B]AIJ support MAT_SKIP_PREALLOCATION for d_nz/o_nz > arguments? > > > > Again, it is fine to push this to release-2.3.3 ? > > > > > > > > -- > > 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