----- Mail original ----- > 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. > > > >> >
