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

[Stefano Zampini]

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