<html><head><meta http-equiv="content-type" content="text/html; charset=utf-8"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br><div><br><blockquote type="cite"><div>On 6 Sep 2023, at 11:27 AM, Carl-Johan Thore <carljohanthore@gmail.com> wrote:</div><br class="Apple-interchange-newline"><div><div dir="ltr"><div>Ok, here's data from a small instance of my problem. It's in Matlab-format which is what I usually use, but let me know if you want something else.</div><div>I can send smaller versions of the problem if you prefer that (really small versions will take a bit of time to code so that one doesn't just get the solution </div><div>0 everywhere). </div><div><br></div><div>One IS, ISO, is included because in my code PETSc computes the other IS as the complement. Also included is the output from PETSc for aij and </div><div>sbaij, and a Matlab-script to check the data.</div><div><br></div><div>Kind regards,</div><div>Carl-Johan<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Sep 6, 2023 at 10:10 AM Carl-Johan Thore <<a href="mailto:carljohanthore@gmail.com" target="_blank">carljohanthore@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>
<div dir="ltr">"Naïve question, but is your problem really symmetric?</div><div dir="ltr">For example, if you do FEM, how do you impose Dirichlet boundary conditions?"</div></div><div><br></div><div>The structure of the matrix is K = [A B; B^T C] with A and C symmetric, so the problem should be symmetric. Numerically, the maximum relative difference between K and K^T is</div><div>on the order of, let's say, 1e-18 which I guess is a good as can be expected? I'm using FEM with (homogeneous) Dirichlet boundary conditions imposed by zeroing rows and columns </div><div>for the fixed DOFs and adding ones to the corresponding diagonal elements, i.e. same as what MatZeroRowsColumns does.
MatZeroRowsColumns
doesn't work for sbaij but we already have an implementation in which rows and columns are cancelled in the element matrices before assembly which is well-tested for the aij-case. We've also verified numerically that aij and sbaij yields the same upper triangular part, and that the RHS is the same in both cases.<br></div><div><br></div><div>"If you can reproduce this on a smaller example and are at liberty to share the Mat, IS, and RHS Vec, feel free to send them at <a href="mailto:petsc-maint@mcs.anl.gov" target="_blank">petsc-maint@mcs.anl.gov</a>, I can have a look."</div><div><br></div><div>Yes, this is reproducible on smaller problems. In case I send you Mat, IS, and RHS, which format is preferable?<br></div><div><br></div><div>Kind regards,</div><div>Carl-Johan<br></div><div><br></div><div><br></div><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Sep 6, 2023 at 8:21 AM Pierre Jolivet <<a href="mailto:pierre.jolivet@lip6.fr" target="_blank">pierre.jolivet@lip6.fr</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="auto"><div dir="ltr"></div><div dir="ltr">Naïve question, but is your problem really symmetric?</div><div dir="ltr">For example, if you do FEM, how do you impose Dirichlet boundary conditions?</div><div dir="ltr">If you can reproduce this on a smaller example and are at liberty to share the Mat, IS, and RHS Vec, feel free to send them at <a href="mailto:petsc-maint@mcs.anl.gov" target="_blank">petsc-maint@mcs.anl.gov</a>, I can have a look.</div><div dir="ltr"><br></div><div dir="ltr">Thanks,</div><div dir="ltr">Pierre</div><div dir="ltr"><br><blockquote type="cite">On 6 Sep 2023, at 8:14 AM, Carl-Johan Thore <<a href="mailto:carljohanthore@gmail.com" target="_blank">carljohanthore@gmail.com</a>> wrote:<br><br></blockquote></div><blockquote type="cite"><div dir="ltr"><div dir="ltr"><div>Thanks for this!</div><div>I now get convergence with sbaij and fieldsplit. However, as can be seen in the attached file, it requires around three times as many outer iterations as with aij to reach the same point (to within rtol). Switching from -pc_fieldsplit_schur_fact_type LOWER to FULL (which is the "correct" choice?) helps quite a bit, but still sbaij takes many more iterations. So it seems that your sbaij-implementation works but that I'm doing something wrong when setting up the fieldsplit pc, or that some of the subsolvers doesn't work properly with sbaij. I've tried bjacobi as you had in your logs, but not hypre yet. But anyway, one doesn't expect the matrix format to have this much impact on convergence? Any other suggestions?<br></div><div>/Carl-Johan<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Sep 4, 2023 at 9:31 PM Pierre Jolivet <<a href="mailto:Pierre.Jolivet@lip6.fr" target="_blank">Pierre.Jolivet@lip6.fr</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div dir="auto">The branch should now be good to go (<a href="https://gitlab.com/petsc/petsc/-/merge_requests/6841" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6841</a>).<div>Sorry, I made a mistake before, hence the error on PetscObjectQuery().</div><div>I’m not sure the code will be covered by the pipeline, but I have tested this on a Raviart—Thomas discretization with PCFIELDSPLIT.</div><div>You’ll see in the attached logs that:</div><div>1) the numerics match</div><div>2) in the SBAIJ case, PCFIELDSPLIT extract the (non-symmetric) A_{01} block from the global (symmetric) A and we get the A_{10} block cheaply by just using MatCreateHermitianTranspose(), instead of calling another time MatCreateSubMatrix()</div><div>Please let me know if you have some time to test the branch and whether it fails or succeeds on your test cases.</div><div><br></div><div>Also, I do not agree with what Hong said.</div><div>Sometimes, the assembly of a coefficient can be more expensive than the communication of the said coefficient.</div><div>So they are instances where SBAIJ would be more efficient than AIJ even if it would require more communication, it is not a black and white picture.</div><div><br></div><div>Thanks,</div><div>Pierre<br><div><br></div><div></div></div></div></div><div><div dir="auto"><div><div></div></div></div></div><div><div dir="auto"><div><div></div><div><br><blockquote type="cite"><div>On 28 Aug 2023, at 12:12 PM, Pierre Jolivet <<a href="mailto:Pierre.Jolivet@lip6.fr" target="_blank">Pierre.Jolivet@lip6.fr</a>> wrote:</div><br><div><div><br><div><blockquote type="cite"><div>On 28 Aug 2023, at 6:50 PM, Carl-Johan Thore <<a href="mailto:carljohanthore@gmail.com" target="_blank">carljohanthore@gmail.com</a>> wrote:</div><br><div><div dir="ltr"><div>I've tried the new files, and with them, PCFIELDSPLIT now gets set up without crashes (but the setup is significantly slower than for MATAIJ)<br></div></div></div></blockquote><div><br></div><div>I’ll be back from Japan at the end of this week, my schedule is too packed to get anything done in the meantime.</div><div>But I’ll let you know when things are working properly (last I checked, I think it was working, but I may have forgotten about a corner case or two).</div><div>But yes, though one would except things to be faster and less memory intensive with SBAIJ, it’s unfortunately not always the case.</div><div><br></div><div>Thanks,</div><div>Pierre</div><br><blockquote type="cite"><div><div dir="ltr"><div>Unfortunately I still get errors later in the process:<br></div><div><br></div><div>[0]PETSC ERROR: Null argument, when expecting valid pointer<br>[0]PETSC ERROR: Null Pointer: Parameter # 1<br>[0]PETSC ERROR: Petsc Development GIT revision: v3.19.4-1023-ga6d78fcba1d GIT Date: 2023-08-22 20:32:33 -0400<br>[0]PETSC ERROR: Configure options -f --with-fortran-bindings=0 --with-cuda --with-cusp --download-scalapack --download-hdf5 --download-zlib --download-mumps --download-parmetis --download-metis --download-ptscotch --download-hypre --download-spai<br>[0]PETSC ERROR: #1 PetscObjectQuery() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/sys/objects/inherit.c:742<br>[0]PETSC ERROR: #2 MatCreateSubMatrix_MPISBAIJ() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/mat/impls/sbaij/mpi/mpisbaij.c:1414<br>[0]PETSC ERROR: #3 MatCreateSubMatrix() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/mat/interface/matrix.c:8476<br>[0]PETSC ERROR: #4 PCSetUp_FieldSplit() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/pc/impls/fieldsplit/fieldsplit.c:826<br>[0]PETSC ERROR: #5 PCSetUp() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/pc/interface/precon.c:1069<br>[0]PETSC ERROR: #6 KSPSetUp() at /mnt/c/mathware/petsc/petsc-v3-19-4/src/ksp/ksp/interface/itfunc.c:415</div><div><br></div><div>The code I'm running here works without any problems for MATAIJ. To run it with MATSBAIJ I've simply used the command-line option</div><div>-dm_mat_type sbaij</div><br><div><br></div><div>Kind regards,</div><div>Carl-Johan<br></div><div><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, Aug 26, 2023 at 5:21 PM Pierre Jolivet via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><br><div><br><blockquote type="cite"><div>On 27 Aug 2023, at 12:14 AM, Carl-Johan Thore <<a href="mailto:carl-johan.thore@liu.se" target="_blank">carl-johan.thore@liu.se</a>> wrote:</div><br><div><div style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">“Well, your A00 and A11 will possibly be SBAIJ also, so you’ll end up with the same issue.”<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">I’m not sure I follow. Does PCFIELDSPLIT extract further submatrices from these blocks, or is there<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">somewhere else in the code that things will go wrong?</div></div></div></blockquote><div><br></div><div>Ah, no, you are right, in that case it should work.</div><br><blockquote type="cite"><div style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">For the MATNEST I was thinking to get some savings from the block-symmetry at least<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">even if symmetry in A00 and A11 cannot be exploited; using SBAIJ for them would just be a<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">(pretty big) bonus.<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">“I’ll rebase on top of main and try to get it integrated if it could be useful to you (but I’m traveling<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">right now so everything gets done more slowly, sorry).”<u></u><u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">Sound great, Thanks again!</div></div></blockquote><div><br></div><div>The MR is there <a href="https://gitlab.com/petsc/petsc/-/merge_requests/6841" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/6841</a>.</div><div>I need to add a new code path in MatCreateRedundantMatrix() to make sure the resulting Mat is indeed SBAIJ, but that is orthogonal to the PCFIELDSPLIT issue.</div><div>The branch should be usable in its current state.</div><div><br></div><div>Thanks,</div><div>Pierre</div><div><br></div><blockquote type="cite"><div style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div><div><div style="border-width:1pt medium medium;border-style:solid none none;border-color:rgb(225,225,225) currentcolor currentcolor;padding:3pt 0cm 0cm"><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><b>From:</b><span> </span>Pierre Jolivet <<a href="mailto:pierre.jolivet@lip6.fr" target="_blank">pierre.jolivet@lip6.fr</a>><span> </span><br><b>Sent:</b><span> </span>Saturday, August 26, 2023 4:36 PM<br><b>To:</b><span> </span>Carl-Johan Thore <<a href="mailto:carl-johan.thore@liu.se" target="_blank">carl-johan.thore@liu.se</a>><br><b>Cc:</b><span> </span>Carl-Johan Thore <<a href="mailto:carljohanthore@gmail.com" target="_blank">carljohanthore@gmail.com</a>>; <a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a><br><b>Subject:</b><span> </span>Re: [petsc-users] PCFIELDSPLIT with MATSBAIJ<u></u><u></u></div></div></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><br><br><u></u><u></u></div><blockquote style="margin-top:5pt;margin-bottom:5pt"><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">On 26 Aug 2023, at 11:16 PM, Carl-Johan Thore <<a href="mailto:carl-johan.thore@liu.se" style="color:blue;text-decoration:underline" target="_blank">carl-johan.thore@liu.se</a>> wrote:<u></u><u></u></div></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div><div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">"(Sadly) MATSBAIJ is extremely broken, in particular, it cannot be used to retrieve rectangular blocks in MatCreateSubMatrices, thus you cannot get the A01 and A10 blocks in PCFIELDSPLIT.<br>I have a branch that fixes this, but I haven’t rebased in a while (and I’m AFK right now), would you want me to rebase and give it a go, or must you stick to a release tarball?"<br><br>Ok, would be great if you could look at this! I don't need to stick to any particular branch.<br><br>Do you think MATNEST could be an alternative here?<u></u><u></u></div></div></div></blockquote><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">Well, your A00 and A11 will possibly be SBAIJ also, so you’ll end up with the same issue.<u></u><u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">I’m using both approaches (monolithic SBAIJ or Nest + SBAIJ), it was crashing but I think it was thoroughly fixed in <a href="https://gitlab.com/petsc/petsc/-/commits/jolivet/feature-matcreatesubmatrices-rectangular-sbaij/" style="color:blue;text-decoration:underline" target="_blank">https://gitlab.com/petsc/petsc/-/commits/jolivet/feature-matcreatesubmatrices-rectangular-sbaij/</a><u></u><u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">It is ugly code on top of ugly code, so I didn’t try to get it integrated and just used the branch locally, and then moved to some other stuff.<u></u><u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">I’ll rebase on top of main and try to get it integrated if it could be useful to you (but I’m traveling right now so everything gets done more slowly, sorry).<u></u><u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><u></u> <u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">Thanks,<u></u><u></u></div></div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">Pierre<u></u><u></u></div></div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif"><br><br><u></u><u></u></div><blockquote style="margin-top:5pt;margin-bottom:5pt"><div><div><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">My matrix is<br>[A00 A01;<br>A01^t A11]<br>so perhaps with MATNEST I can make use of the block-symmetry at least, and then use MATSBAIJ for<span> </span><br>A00 and A11 if it's possible to combine matrix types which the manual seems to imply.<span> </span><br><br>Kind regards<br>Carl-Johan<br><br><br><br><u></u><u></u></div><blockquote style="margin-top:5pt;margin-bottom:5pt"><div style="margin:0cm;font-size:11pt;font-family:Calibri,sans-serif">On 26 Aug 2023, at 10:09 PM, Carl-Johan Thore <<a href="mailto:carljohanthore@gmail.com" style="color:blue;text-decoration:underline" target="_blank">carljohanthore@gmail.com</a>> wrote:<br><br>Hi,<br><br>I'm trying to use PCFIELDSPLIT with MATSBAIJ in PETSc 3.19.4.<span> </span><br>According to the manual "[t]he fieldsplit preconditioner cannot<span> </span><br>currently be used with the MATBAIJ or MATSBAIJ data formats if the<span> </span><br>blocksize is larger than 1". Since my blocksize is exactly 1 it would seem that I can use PCFIELDSPLIT. But this fails with "PETSC ERROR: For symmetric format, iscol must equal isrow"<br>from MatCreateSubMatrix_MPISBAIJ. Tracing backwards one ends up in<span> </span><br>fieldsplit.c at<br><br>/* extract the A01 and A10 matrices */ ilink = jac->head;<span> </span><br>PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis)); if<span> </span><br>(jac->offdiag_use_amat) { PetscCall(MatCreateSubMatrix(pc->mat,<span> </span><br>ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B)); } else {<br> PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis,<span> </span><br>MAT_INITIAL_MATRIX, &jac->B)); }<br><br>This, since my A01 and A10 are not square, seems to explain why iscol is not equal to isrow.<br>From this I gather that it is in general NOT possible to use<span> </span><br>PCFIELDSPLIT with MATSBAIJ even with block size 1?<br><br>Kind regards,<br>Carl-Johan</div></blockquote></div></div></blockquote></div></div></blockquote></div><br></div></blockquote></div></div>
</div></blockquote></div><br></div></div></blockquote></div><br></div></div></div></blockquote></div>
<div><experiment_with_sbaij.txt></div></div></blockquote></div></blockquote></div></div></blockquote></div></div></blockquote></div><br></div></body></html>