[petsc-dev] [petsc-users] How to compute RARt with A and R as distributed (MPI) matrices ?

[Franck Houssen]

Great ! 

----- Mail original -----

> De: "Hong" <hzhang at mcs.anl.gov>
> À: "Franck Houssen" <franck.houssen at inria.fr>
> Cc: "Stefano Zampini" <stefano.zampini at gmail.com>, "petsc-dev"
> <petsc-dev at mcs.anl.gov>
> 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" < hzhang at mcs.anl.gov >
> > 
> 
> > > À: "Stefano Zampini" < stefano.zampini at gmail.com >
> > 
> 
> > > Cc: "Franck Houssen" < franck.houssen at inria.fr >, "petsc-dev" <
> > > petsc-dev at mcs.anl.gov >
> > 
> 
> > > 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 <
> > > stefano.zampini at gmail.com
> > > > wrote:
> > 
> 

> > > > > On Jun 23, 2017, at 10:10 AM, Franck Houssen <
> > > > > franck.houssen at inria.fr
> > > > > >
> > > > > 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" < jed at jedbrown.org >
> > > 
> > 
> 
> > > > >> À: "Hong" < hzhang at mcs.anl.gov >
> > > 
> > 
> 
> > > > >> Cc: "petsc-dev" < petsc-dev at mcs.anl.gov >, "PETSc users list" <
> > > > >> petsc-users at mcs.anl.gov >
> > > 
> > 
> 
> > > > >> 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 < hzhang at mcs.anl.gov > 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.
> > > 
> > 
> 
> > > > >>
> > > 
> > 
> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170625/e371418c/attachment.html>