Great ! ----- Mail original -----
> De: "Hong" <[email protected]> > À: "Franck Houssen" <[email protected]> > Cc: "Stefano Zampini" <[email protected]>, "petsc-dev" > <[email protected]> > Envoyé: Vendredi 23 Juin 2017 18:08:08 > Objet: Re: [petsc-dev] [petsc-users] How to compute RARt with A and R as > distributed (MPI) matrices ? > Franck: > See > https://bitbucket.org/petsc/petsc/commits/f3dcd100076caffe0af33a424f0d06cb26a76989 > Hong > > > De: "Hong" < [email protected] > > > > > > > À: "Stefano Zampini" < [email protected] > > > > > > > Cc: "Franck Houssen" < [email protected] >, "petsc-dev" < > > > [email protected] > > > > > > > Envoyé: Vendredi 23 Juin 2017 16:59:42 > > > > > > Objet: Re: [petsc-dev] [petsc-users] How to compute RARt with A and R as > > > distributed (MPI) matrices ? > > > > > > Due to petsc sparse matrix block row distribution among processes, > > > > > > we do not have an efficient way to implement parallel R * A * R^T > > > > > > (explicit transpose of R is very expensive). > > > > > This sentence above is a perfect candidate to add in the doc : 1 line, 1 > > sentence that enables the user to understand WHY he must not use matRARt if > > he has a distributed matrix in hand (= to understand WHY this is not > > implemented and will throw an error). > > > > For P^t * A * P=P^t *(AP), we do local transpose P^t and sum (partical > > > P^t_local *(AP). > > > > > > I agree that a completed documentation would help users; but in reality, > > > with > > > so few developers, so many PETSc functions collected over the years, each > > > of > > > function has various implementations, implementations are changes > > > frequently > > > .... > > > > > > We keep object operation tables in PETSc. When user invokes a > > > non-supported > > > implementation, the execution throughs an error message. User is informed > > > from testing instead of reading user manual to find out which detailed > > > implementation is supported. We also provide instant email support for > > > user's request. > > > > > > We'll try to maintain an update user manual and welcome bug/error report. > > > I'll check the manual on MatRARt/PtAP to make it clear to users. > > > > > > Hong > > > > > > On Fri, Jun 23, 2017 at 7:13 AM, Stefano Zampini < > > > [email protected] > > > > wrote: > > > > > > > > On Jun 23, 2017, at 10:10 AM, Franck Houssen < > > > > > [email protected] > > > > > > > > > > > wrote: > > > > > > > > > > > > > > > > > > > > > > By the way, I noticed this is exactly the same for > > > > > MatRARt/PtAPNumeric/Symbolic: pages of the documentation are > > > > > identical, > > > > > so > > > > > the user can not understand why to use the RARt or the PtAP. > > > > > > > > > > MatRARt computes: R * A * R^T > > > > > > > > > > MatPtAP computes: P^t * A * P > > > > > > > > > > they perform different operations > > > > > > > > > > > I believe adding a short note in MatRARt/PtAP pages and > > > > > MatRARt/PtAPNumeric/Symbolic would be welcome to say "this method is > > > > > meant > > > > > to do/support this but not that”. > > > > > > > > > > An operation like RARt (or PtAP) may have a specialized implementation, > > > > or > > > > be > > > > available though composition of operatotions, i.e. PtAP = P^t * (A * P) > > > > or > > > > compute transpose of P (explicitly) and then call the triple matrix > > > > product. > > > > > > > > > > It’s not possible to keep track of all the possible cases. > > > > > > > > > > > > > > > > > > > > > > Franck > > > > > > > > > > > > > > > > > > > > > > ----- Mail original ----- > > > > > > > > > > >> De: "Jed Brown" < [email protected] > > > > > > > > > > > >> À: "Hong" < [email protected] > > > > > > > > > > > >> Cc: "petsc-dev" < [email protected] >, "PETSc users list" < > > > > >> [email protected] > > > > > > > > > > > >> Envoyé: Jeudi 22 Juin 2017 18:17:33 > > > > > > > > > > >> Objet: Re: [petsc-dev] [petsc-users] How to compute RARt with A and > > > > >> R > > > > >> as > > > > >> distributed (MPI) matrices ? > > > > > > > > > > >> > > > > > > > > > > >> Hong < [email protected] > writes: > > > > > > > > > > >> > > > > > > > > > > >>> Jed: > > > > > > > > > > >>>> > > > > > > > > > > >>>>>> Is it better this way or as a fallback when !A->ops->rart? > > > > >>>>>> MatPtAP > > > > > > > > > > >>>>>> handles other combinations like MAIJ. > > > > > > > > > > >>>>>> > > > > > > > > > > >>>>> > > > > > > > > > > >>>>> Do you mean > > > > > > > > > > >>>>> if ( !A->ops->rart) { > > > > > > > > > > >>>>> Mat Rt; > > > > > > > > > > >>>>> ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&Rt);CHKERRQ(ierr); > > > > > > > > > > >>>>> ierr = MatMatMatMult(R,A,Rt,scall,fill,C);CHKERRQ(ierr); > > > > > > > > > > >>>>> ierr = MatDestroy(&Rt);CHKERRQ(ierr); > > > > > > > > > > >>>>> } > > > > > > > > > > >>>>> This does NOT work for most matrix formats because we do not have > > > > > > > > > > >>>> fallbacks > > > > > > > > > > >>>>> for MatTranspose() and MatMatMult(). > > > > > > > > > > >>>> > > > > > > > > > > >>>> That's fine; they'll trigger an error and we'll be able to see > > > > >>>> from > > > > >>>> the > > > > > > > > > > >>>> stack that it can be made to work by either implementing the > > > > >>>> appropriate > > > > > > > > > > >>>> MatRARt or MatTranspose and MatMatMatMult. > > > > > > > > > > >>>> > > > > > > > > > > >>> > > > > > > > > > > >>> You prefer adding this default, even though it gives error in > > > > >>> either > > > > > > > > > > >>> MatTranspose() or MatMatMatMult() depends on input matrix format? > > > > > > > > > > >> > > > > > > > > > > >> Yeah, in the sense that it gives more opportunities to succeed. > > > > > > > > > > >> > > > > > > > > > > >>> If so, we need add this type of 'default' to all mat operations -- > > > > > > > > > > >>> currently, all routines do > > > > > > > > > > >>> if (!mat->ops-> ) > > > > > > > > > > >>> SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Mat type > > > > > > > > > > >>> %s",((PetscObject)mat)->type_name); > > > > > > > > > > >> > > > > > > > > > > >> Probably. > > > > > > > > > > >> > > > > > >
