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

[Hong]

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

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/20170623/29749f63/attachment.html>