<div dir="ltr">I added following to the help menu of MatRARt():<div><div>   "This routine is currently only implemented for pairs of SeqAIJ matrices and classes</div><div>   which inherit from SeqAIJ. For parallel matrices, it is recommended to use MatPtAP()."</div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jun 23, 2017 at 9:59 AM, Hong <span dir="ltr"><<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Due to petsc sparse matrix block row distribution among processes,<div>we do not have an efficient way to implement parallel R * A * R^T </div><div>(explicit transpose of R is very expensive).</div><div><br></div><div>For P^t * A * P=P^t *(AP), we do local transpose P^t and sum (partical P^t_local *(AP). </div><div><br></div><div>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 ....</div><div><br></div><div>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. </div><div><br></div><div>We'll try to maintain an update user manual and welcome bug/error report. I'll check the manual on <span style="font-size:12.8px">MatRARt/PtAP to make it clear to users.</span></div><span class="HOEnZb"><font color="#888888"><div><span style="font-size:12.8px"><br></span></div></font></span><div><span class="HOEnZb"><font color="#888888"><span style="font-size:12.8px">Hong<br></span></font></span><div><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jun 23, 2017 at 7:13 AM, Stefano Zampini <span dir="ltr"><<a href="mailto:stefano.zampini@gmail.com" target="_blank">stefano.zampini@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><span class="m_4031888867791915119gmail-"><br>
> On Jun 23, 2017, at 10:10 AM, Franck Houssen <<a href="mailto:franck.houssen@inria.fr" target="_blank">franck.houssen@inria.fr</a>> wrote:<br>
><br>
> 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.<br>
<br>
</span>MatRARt computes: R * A * R^T<br>
MatPtAP computes: P^t * A * P<br>
<br>
they perform different operations<br>
<span class="m_4031888867791915119gmail-"><br>
> 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”.<br>
<br>
</span>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.<br>
It’s not possible to keep track of all the possible cases.<br>
<div class="m_4031888867791915119gmail-HOEnZb"><div class="m_4031888867791915119gmail-h5"><br>
><br>
<br>
> Franck<br>
><br>
> ----- Mail original -----<br>
>> De: "Jed Brown" <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>><br>
>> À: "Hong" <<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>><br>
>> Cc: "petsc-dev" <<a href="mailto:petsc-dev@mcs.anl.gov" target="_blank">petsc-dev@mcs.anl.gov</a>>, "PETSc users list" <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>><br>
>> Envoyé: Jeudi 22 Juin 2017 18:17:33<br>
>> Objet: Re: [petsc-dev] [petsc-users] How to compute RARt with A and R as     distributed (MPI) matrices ?<br>
>><br>
>> Hong <<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>> writes:<br>
>><br>
>>> Jed:<br>
>>>><br>
>>>>>> Is it better this way or as a fallback when !A->ops->rart?  MatPtAP<br>
>>>>>> handles other combinations like MAIJ.<br>
>>>>>><br>
>>>>><br>
>>>>> Do you mean<br>
>>>>> if ( !A->ops->rart) {<br>
>>>>>    Mat Rt;<br>
>>>>>    ierr = MatTranspose(R,MAT_INITIAL_MAT<wbr>RIX,&Rt);CHKERRQ(ierr);<br>
>>>>>    ierr = MatMatMatMult(R,A,Rt,scall,fil<wbr>l,C);CHKERRQ(ierr);<br>
>>>>>    ierr = MatDestroy(&Rt);CHKERRQ(ierr);<br>
>>>>> }<br>
>>>>> This does NOT work for most matrix formats because we do not have<br>
>>>> fallbacks<br>
>>>>> for MatTranspose() and MatMatMult().<br>
>>>><br>
>>>> That's fine; they'll trigger an error and we'll be able to see from the<br>
>>>> stack that it can be made to work by either implementing the appropriate<br>
>>>> MatRARt or MatTranspose and MatMatMatMult.<br>
>>>><br>
>>><br>
>>> You prefer adding this default, even though it gives error in either<br>
>>> MatTranspose() or MatMatMatMult() depends on input matrix format?<br>
>><br>
>> Yeah, in the sense that it gives more opportunities to succeed.<br>
>><br>
>>> If so, we need add this type of 'default' to all mat operations --<br>
>>> currently, all routines do<br>
>>> if (!mat->ops-> )<br>
>>> SETERRQ1(PetscObjectComm((Pets<wbr>cObject)mat),PETSC_ERR_SUP,"<wbr>Mat type<br>
>>> %s",((PetscObject)mat)->type_n<wbr>ame);<br>
>><br>
>> Probably.<br>
>><br>
<br>
</div></div></blockquote></div><br></div></div></div></div></div>
</blockquote></div><br></div>