<div dir="ltr"><div dir="ltr">On Wed, Sep 16, 2020 at 10:57 AM Abhyankar, Shrirang G via petsc-dev <<a href="mailto:petsc-dev@mcs.anl.gov">petsc-dev@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div lang="EN-US">
<div class="gmail-m_-6499821183657313774WordSection1">
<div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(181,196,223);padding:3pt 0in 0in">
<p class="MsoNormal" style="margin-left:0.5in"><b><span style="font-size:12pt;color:black">From:
</span></b><span style="font-size:12pt;color:black">Pierre Jolivet <<a href="mailto:pierre.jolivet@enseeiht.fr" target="_blank">pierre.jolivet@enseeiht.fr</a>><br>
<b>Date: </b>Wednesday, September 16, 2020 at 12:59 AM<br>
<b>To: </b>Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>><br>
<b>Cc: </b>"Abhyankar, Shrirang G" <<a href="mailto:shrirang.abhyankar@pnnl.gov" target="_blank">shrirang.abhyankar@pnnl.gov</a>>, PETSc Development <<a href="mailto:petsc-dev@mcs.anl.gov" target="_blank">petsc-dev@mcs.anl.gov</a>><br>
<b>Subject: </b>Re: [petsc-dev] PDIPDM questions<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
</div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
<div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<p class="MsoNormal" style="margin-left:0.5in">On 15 Sep 2020, at 11:18 PM, Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<u></u><u></u></p>
</div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><br>
<span style="font-size:9pt;font-family:Helvetica"> I see now, the TaoSetup_PDIPM() code explicitly builds the large matrix by calling MatGetRow() on the smaller matrices, no matrix abstractions at all and assuming a global AIJ storage of the large matrix. </span><u></u><u></u></p>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> I think this can be made more general by building a MatNest() from the smaller matrices which allows the smaller matrices to have whatever unique structure is
appropriate for them. Of course, if needed, for example, for a direct solver the MatNest can convert itself to a global AIJ matrix and get the current result.<u></u><u></u></span></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in">Before jumping into the PDIPM train, I thought a MatNest was used internally. Like you can solve [[A, b],[b^T, 0]] (b of dimension M x 1) with a MatNest + fieldsplit or just MatConvert() + LU as you mentioned.<u></u><u></u></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif">The current implementation of PDIPM forms a monolithic KKT matrix in AIJ format and then passes that to the linear solver. We did try using fieldsplit by setting ISes for
primal and dual variables. But, the performance of the solver using Fieldsplit was not that great. We’ll probably try it again after the current development of PDIPM.</span></p></div></div></div></div></blockquote><div><br></div><div>What solver are you using? There is a nice paper by Yin Zhang from a while ago about these. There must be something more recent.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt </div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div lang="EN-US"><div class="gmail-m_-6499821183657313774WordSection1"><div><p class="MsoNormal" style="margin-left:0.5in">
<u></u><u></u></p>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> As Pierre rightly points out, MPIAIJ is just not a good format for many constrained optimization problems. This has been something seriously neglected in PETSc/TAO.
I don't think any single format is ideal for such problems, likely hybrid formats are needed and this gets super tricky when one needs to use direct solvers such a Cholesky. If MUMPS has non-row oriented storage formats we could leverage that but I don't know
if it does.<u></u><u></u></span></p>
</div>
</div>
</blockquote>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in">Just FYI, MUMPS is the most general (direct solver) when it comes to input storage. It handles non-sorted coordinate format, possibly with duplicated entries which are then summed among processes (so, in theory,
it could be used to “invert" a MatIS/MatNest, without having to first convert the MatIS/MatNest to MatAIJ).<u></u><u></u></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif">Thanks for the pointer on MUMPS’s non-sorted coordinate format.
<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><u></u> <u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in">Thanks,<u></u><u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in">Pierre<u></u><u></u></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif"><u></u> <u></u></span></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif">Thanks,<u></u><u></u></span></p>
<p class="MsoNormal"><span style="font-size:12pt;font-family:"Times New Roman",serif">Shri<u></u><u></u></span></p>
</div>
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<u></u><u></u></p>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> Barry<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRange(tao->hessian,&rjstart,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i<pdipm->nx; i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= cranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - cranges[proc] + Jranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(tao->hessian,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->ng) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Insert grad g' */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRanges(tao->jacobian_equality,&ranges);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* find row ownership of */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= ranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> nx_all = rranges[proc+1] - rranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - ranges[proc] + Jranges[proc] + nx_all;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(jac_equality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Insert Jce_xfixed^T' */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->nxfixed) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRanges(pdipm->Jce_xfixed,&ranges);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* find row ownership of */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= ranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> nx_all = rranges[proc+1] - rranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + ng_all[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(Jce_xfixed_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->nh) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Insert -grad h' */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRanges(tao->jacobian_inequality,&ranges);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* find row ownership of */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= ranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> nx_all = rranges[proc+1] - rranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(jac_inequality_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Insert Jci_xb^T' */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRanges(pdipm->Jci_xb,&ranges);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* find row ownership of */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= ranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> nx_all = rranges[proc+1] - rranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - ranges[proc] + Jranges[proc] + nx_all + nce_all[proc] + nh_all[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(Jci_xb_trans,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* 2nd Row block of KKT matrix: [grad Ce, 0, 0, 0] */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->Ng) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRange(tao->jacobian_equality,&rjstart,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < pdipm->ng; i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_lambdae + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= cranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - cranges[proc] + Jranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad g */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(tao->jacobian_equality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Jce_xfixed */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->Nxfixed) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRange(pdipm->Jce_xfixed,&Jcrstart,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < (pdipm->nce - pdipm->ng); i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_lambdae + pdipm->ng + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> j = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (cols[j] >= cranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = cols[j] - cranges[proc] + Jranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(pdipm->Jce_xfixed,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* 3rd Row block of KKT matrix: [ gradCi, 0, 0, -I] */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (pdipm->Nh) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRange(tao->jacobian_inequality,&rjstart,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < pdipm->nh; i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_lambdai + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (aj[j] >= cranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = aj[j] - cranges[proc] + Jranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr); /* grad h */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(tao->jacobian_inequality,i+rjstart,&nc,&aj,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* -I */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < pdipm->nh; i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_lambdai + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = rstart + pdipm->off_z + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* Jci_xb */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetOwnershipRange(pdipm->Jci_xb,&Jcrstart,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < (pdipm->nci - pdipm->nh); i++){<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_lambdai + pdipm->nh + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatGetRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> if (nc != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"nc != 1");<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> proc = 0;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (j=0; j < nc; j++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> while (cols[j] >= cranges[proc+1]) proc++;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = cols[j] - cranges[proc] + Jranges[proc];<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatRestoreRow(pdipm->Jci_xb,i+Jcrstart,&nc,&cols,NULL);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* -I */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> col = rstart + pdipm->off_z + pdipm->nh + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,1,&col,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> /* 4-th Row block of KKT matrix: Z and Ci */<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> for (i=0; i < pdipm->nci; i++) {<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> row = rstart + pdipm->off_z + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> cols1[0] = rstart + pdipm->off_lambdai + i;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> cols1[1] = row;<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> ierr = MatPreallocateSet(row,2,cols1,dnz,onz);CHKERRQ(ierr);<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"> }<u></u><u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><br>
<br>
<u></u><u></u></span></p>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica">On Sep 15, 2020, at 4:03 PM, Abhyankar, Shrirang G <<a href="mailto:shrirang.abhyankar@pnnl.gov" target="_blank">shrirang.abhyankar@pnnl.gov</a>> wrote:<u></u><u></u></span></p>
</div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:9pt;font-family:Helvetica"><u></u> <u></u></span></p>
<div>
<div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">There is some code in PDIPM that copies over matrix elements from user jacobians/hessians to the big KKT matrix. I am not sure how/if that will
work with MatMPI_Column. We’ll have to try out and we need an example for that<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><span style="font-size:12pt;font-family:"Apple Color Emoji"">😊</span><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
<div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Thanks,</span><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Shri </span><u></u><u></u></p>
</div>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
<div>
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
<div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(181,196,223);padding:3pt 0in 0in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><b><span style="font-size:12pt">From:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span></b><span style="font-size:12pt">Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>><br>
<b>Date:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>Tuesday, September 15, 2020 at 1:21 PM<br>
<b>To:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>Pierre Jolivet <<a href="mailto:pierre.jolivet@enseeiht.fr" target="_blank">pierre.jolivet@enseeiht.fr</a>><br>
<b>Cc:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>"Abhyankar, Shrirang G" <<a href="mailto:shrirang.abhyankar@pnnl.gov" target="_blank">shrirang.abhyankar@pnnl.gov</a>>, PETSc Development <<a href="mailto:petsc-dev@mcs.anl.gov" target="_blank">petsc-dev@mcs.anl.gov</a>><br>
<b>Subject:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>Re: [petsc-dev] PDIPDM questions</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> Pierre,<u></u><u></u></p>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> Based on my previous mail I am hoping that the <span style="font-size:12pt;font-family:"Times New Roman",serif">PDIPM algorithm itself won't need a major refactorization to be scalable, only a custom matrix
type is needed to store and compute with the Hessian in a scalable way.</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> Barry</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<br>
<u></u><u></u></p>
</div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">On Sep 15, 2020, at 12:50 PM, Pierre Jolivet <<a href="mailto:pierre.jolivet@enseeiht.fr" target="_blank">pierre.jolivet@enseeiht.fr</a>> wrote:<u></u><u></u></p>
</div>
</div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<br>
<u></u><u></u></p>
</div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">On 15 Sep 2020, at 5:40 PM, Abhyankar, Shrirang G <<a href="mailto:shrirang.abhyankar@pnnl.gov" target="_blank">shrirang.abhyankar@pnnl.gov</a>> wrote:<u></u><u></u></p>
</div>
</div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Pierre,</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> You are right. There are a few MatMultTransposeAdd that may need conforming layouts for the equality/inequality constraint vectors and equality/inequality
constraint Jacobian matrices. I need to check if that’s the case. We only have ex1 example currently, we need to add more examples. We are currently working on making PDIPM robust and while doing it will work on adding another example.</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
</div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Very naive question, but given that I have a single constraint, how do I split a 1 x N matrix column-wise? I thought it was not possible.<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><i><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span></i><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">When setting the size of the constraint vector, you need to set the local size on one rank to 1 and all others to zero. For the Jacobian, the local
row size on that rank will be 1 and all others to zero. The column layout for the Jacobian should follow the layout for vector x. So each rank will set the local column size of the Jacobian to local size of x.</span><u></u><u></u></p>
</div>
</div>
</div>
</blockquote>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">That is assuming I don’t want x to follow the distribution of the Hessian, which is not my case.<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Is there some plan to make PDIPM handle different layouts?<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">I hope I’m not the only one thinking that having a centralized Hessian when there is a single constraint is not scalable?<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Thanks,<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Pierre<u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Shri</span><u></u><u></u></p>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
</div>
<div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">On 15 Sep 2020, at 2:21 AM, Abhyankar, Shrirang G <<a href="mailto:shrirang.abhyankar@pnnl.gov" target="_blank">shrirang.abhyankar@pnnl.gov</a>> wrote:<u></u><u></u></p>
</div>
</div>
</div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Hello Pierre,</span><u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> PDIPM works in parallel so you can have distributed Hessian, Jacobians, constraints, variables, gradients in any layout you want. If you are
using a DM then you can have it generate the Hessian.<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><u></u><u></u></p>
</div>
</div>
</div>
</div>
</blockquote>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Could you please show an example where this is the case?<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">pdipm->x, which I’m assuming is a working vector, is both used as input for Hessian and Jacobian functions, e.g., <a href="https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L369" target="_blank">https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L369</a> (Hessian)
+ <a href="https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L473" target="_blank">https://gitlab.com/petsc/petsc/-/blob/master/src/tao/constrained/impls/ipm/pdipm.c#L473</a> (Jacobian)<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">I thus doubt that it is possible to have different layouts?<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">In practice, I end up with the following error when I try this (2 processes, distributed Hessian with centralized Jacobian):<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: Nonconforming object sizes<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: Vector wrong size 14172 for scatter 0 (scatter reverse and vector to != ctx from size)<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: #1 VecScatterBegin() line 96 in /Users/jolivet/Documents/repositories/petsc/src/vec/vscat/interface/vscatfce.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: #2 MatMultTransposeAdd_MPIAIJ() line 1223 in /Users/jolivet/Documents/repositories/petsc/src/mat/impls/aij/mpi/mpiaij.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: #3 MatMultTransposeAdd() line 2648 in /Users/jolivet/Documents/repositories/petsc/src/mat/interface/matrix.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: Nonconforming object sizes<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: Vector wrong size 13790 for scatter 27962 (scatter reverse and vector to != ctx from size)<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: #4 TaoSNESFunction_PDIPM() line 510 in /Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #5 TaoSolve_PDIPM() line 712 in /Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[1]PETSC ERROR: #6 TaoSolve() line 222 in /Users/jolivet/Documents/repositories/petsc/src/tao/interface/taosolver.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #1 VecScatterBegin() line 96 in /Users/jolivet/Documents/repositories/petsc/src/vec/vscat/interface/vscatfce.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #2 MatMultTransposeAdd_MPIAIJ() line 1223 in /Users/jolivet/Documents/repositories/petsc/src/mat/impls/aij/mpi/mpiaij.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #3 MatMultTransposeAdd() line 2648 in /Users/jolivet/Documents/repositories/petsc/src/mat/interface/matrix.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #4 TaoSNESFunction_PDIPM() line 510 in /Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #5 TaoSolve_PDIPM() line 712 in /Users/jolivet/Documents/repositories/petsc/src/tao/constrained/impls/ipm/pdipm.c<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">[0]PETSC ERROR: #6 TaoSolve() line 222 in /Users/jolivet/Documents/repositories/petsc/src/tao/interface/taosolver.c<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">I think this can be reproduced by ex1.c by just distributing the Hessian instead of having it centralized on rank 0.<u></u><u></u></p>
</div>
</div>
</div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<br>
<br>
<u></u><u></u></p>
</div>
</div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Ideally, you want to have the layout below to minimize movement of matrix/vector elements across ranks.</span><u></u><u></u></p>
</div>
</div>
</div>
<p class="gmail-m_-6499821183657313774MsoListParagraph" style="margin-right:0in;margin-bottom:0in;margin-left:2in">
<span style="font-size:10pt;font-family:Symbol">·</span><span style="font-size:7pt;font-family:"Times New Roman",serif"> <span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><span style="font-size:12pt;font-family:"Times New Roman",serif">The layout
of vectors x, bounds on x, and gradient is same.</span><u></u><u></u></p>
<p class="gmail-m_-6499821183657313774MsoListParagraph" style="margin-right:0in;margin-bottom:0in;margin-left:2in">
<span style="font-size:10pt;font-family:Symbol">·</span><span style="font-size:7pt;font-family:"Times New Roman",serif"> <span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><span style="font-size:12pt;font-family:"Times New Roman",serif">The row
layout of the equality/inequality Jacobian is same as the equality/inequality constraint vector layout.</span><u></u><u></u></p>
<p class="gmail-m_-6499821183657313774MsoListParagraph" style="margin-right:0in;margin-bottom:0in;margin-left:2in">
<span style="font-size:10pt;font-family:Symbol">·</span><span style="font-size:7pt;font-family:"Times New Roman",serif"> <span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><span style="font-size:12pt;font-family:"Times New Roman",serif">The column
layout of the equality/inequality Jacobian is same as that for x.</span><u></u><u></u></p>
</div>
</blockquote>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Very naive question, but given that I have a single constraint, how do I split a 1 x N matrix column-wise? I thought it was not possible.<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Thanks,<u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Pierre<u></u><u></u></p>
</div>
</div>
</div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<br>
<br>
<u></u><u></u></p>
</div>
</div>
<blockquote style="margin-top:5pt;margin-bottom:5pt">
<div>
<p class="gmail-m_-6499821183657313774MsoListParagraph" style="margin-right:0in;margin-bottom:0in;margin-left:2in">
<span style="font-size:10pt;font-family:Symbol">·</span><span style="font-size:7pt;font-family:"Times New Roman",serif"> <span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span><span style="font-size:12pt;font-family:"Times New Roman",serif">The row
and column layout for the Hessian is same as x.</span><u></u><u></u></p>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">The tutorial example ex1 is extremely small (only 2 variables) so its implementation is very simplistic. I think, in parallel, it ships off constraints
etc. to rank 0. It’s not an ideal example w.r.t demonstrating a parallel implementation. We aim to add more examples as we develop PDIPM. If you have an example to contribute then we would most welcome it and provide help on adding it.</span><u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif"> </span><u></u><u></u></p>
</div>
</div>
</div>
<div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Thanks,</span><u></u><u></u></p>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><span style="font-size:12pt;font-family:"Times New Roman",serif">Shri</span><u></u><u></u></p>
</div>
</div>
</div>
</div>
<div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(181,196,223);padding:3pt 0in 0in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"><b><span style="font-size:12pt">From:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></span></b><span style="font-size:12pt">petsc-dev <<a href="mailto:petsc-dev-bounces@mcs.anl.gov" target="_blank">petsc-dev-bounces@mcs.anl.gov</a>>
on behalf of Pierre Jolivet <<a href="mailto:pierre.jolivet@enseeiht.fr" target="_blank">pierre.jolivet@enseeiht.fr</a>><br>
<b>Date:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>Monday, September 14, 2020 at 1:52 PM<br>
<b>To:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>PETSc Development <<a href="mailto:petsc-dev@mcs.anl.gov" target="_blank">petsc-dev@mcs.anl.gov</a>><br>
<b>Subject:<span class="gmail-m_-6499821183657313774apple-converted-space"> </span></b>[petsc-dev] PDIPDM questions</span><u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Hello,<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">In my quest to help users migrate from Ipopt to Tao, I’ve a new question.<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">When looking at src/tao/constrained/tutorials/ex1.c, it seems that almost everything is centralized on rank 0 (local sizes are 0 but on rank 0).<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">I’d like to have my Hessian distributed more naturally, as in (almost?) all other SNES/TS examples, but still keep the Jacobian of my equality constraint, which is of dimension 1 x N (N >> 1), centralized on rank
0.<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Is this possible?<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">If not, is it possible to supply the transpose of the Jacobian, of dimension N x 1, which could then be distributed row-wise like the Hessian?<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Or maybe use some trick to distribute a MatAIJ/MatDense of dimension 1 x N column-wise? Use a MatNest with as many blocks as processes?<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">So, just to sum up, how can I have a distributed Hessian with a Jacobian with a single row?<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in"> <u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Thanks in advance for your help,<u></u><u></u></p>
</div>
</div>
</div>
</div>
<div>
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<div style="margin-left:0.5in">
<p class="MsoNormal" style="margin-left:0.5in">Pierre<u></u><u></u></p>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</blockquote>
</div>
<p class="MsoNormal" style="margin-left:0.5in"><br>
<br>
<u></u><u></u></p>
</div>
</div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>