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