<html><head><meta http-equiv="Content-Type" content="text/html; charset=us-ascii"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><br class=""><div> Olivier,</div><div><br class=""></div><div> I believe the code is ready for your use. MatInvertVariableBlockEnvelope() will take an MPIAIJ matrix (this would be your C*C' matrix, determine the block diagonal envelope automatically and construct a new MPIAIJ matrix that is the inverse of the block diagonal envelope (by inverting each of the blocks). For small blocks, as is your case, the inverse should inexpensive. </div><div><br class=""></div><div> I have a test code src/mat/tests/ex178.c that demonstrates its use and tests if for some random blocking across multiple processes.</div><div><br class=""></div><div> Please let me know if you have any difficulties. The branch is the same name as before <b style="color: rgb(200, 20, 201); font-family: Menlo; font-size: 14px;" class="">barry/2020-10-08/invert-block-diagonal-aij </b><span style="color: rgb(200, 20, 201); font-family: Menlo; font-size: 14px;" class=""> </span>but it has been rebased so you will need to delete your local copy of the branch and then recreate it after git fetch.</div><div><br class=""></div><div> Barry</div><div><br class=""></div><div><br class=""></div><div><br class=""><blockquote type="cite" class=""><div class="">On Sep 8, 2021, at 3:09 PM, Olivier Jamond <<a href="mailto:olivier.jamond@cea.fr" class="">olivier.jamond@cea.fr</a>> wrote:</div><br class="Apple-interchange-newline"><div class="">
<meta http-equiv="Content-Type" content="text/html;
charset=windows-1252" class="">
<div class=""><p class="">Ok thanks a lot! I look forward to hearing from you!</p><p class="">Best regards,<br class="">
Olivier<br class="">
</p>
<div class="moz-cite-prefix">On 08/09/2021 20:56, Barry Smith wrote:<br class="">
</div>
<blockquote type="cite" cite="mid:122027EE-C4EC-4DCC-832A-ED0FBD092B30@petsc.dev" class="">
<meta http-equiv="Content-Type" content="text/html;
charset=windows-1252" class="">
<div class=""><br class="">
</div>
Olivier,
<div class=""><br class="">
</div>
<div class=""> I'll refresh my memory on this and see if I can
finish it up.</div>
<div class=""><br class="">
</div>
<div class=""> Barry</div>
<div class=""><br class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Sep 2, 2021, at 12:38 PM, Olivier Jamond
<<a href="mailto:olivier.jamond@cea.fr" class="" moz-do-not-send="true">olivier.jamond@cea.fr</a>>
wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div class=""><p class="">Dear Barry,</p><p class="">First I hope that you and your relatives are
doing well in these troubled times...<br class="">
</p><p class="">I allow myself to come back to you about the
subject of being able to compute something like
'(C*Ct)^(-1)*C', where 'C' is a 'MPC' matrix that is
used to impose some boundary conditions for a
structural finite element problem:</p><p class=""><font class="" face="monospace">[K
C^t][U]=[F]<br class="">
[C 0 ][L] [D]</font></p><p class="">as we discussed some time ago, I would like
to solve such a problem using the Ainsworth method,
which involves this '(C*Ct)^(-1)*C'.</p><p class="">You kindly started some developments to help
me on that, which worked as a 'proof of concept' in
sequential, but not yet in parallel, and also kindly
suggested that you could extend it to the parallel
case (MR: <a class="moz-txt-link-freetext" href="https://gitlab.com/petsc/petsc/-/merge_requests/3544" moz-do-not-send="true">https://gitlab.com/petsc/petsc/-/merge_requests/3544</a>).
Can this be still 'scheduled' on your side?</p><p class="">Sorry again to "harass" you about that...</p><p class="">Best regards,<br class="">
Olivier Jamond<br class="">
</p>
<div class="moz-cite-prefix">On 03/02/2021 08:45,
Olivier Jamond wrote:<br class="">
</div>
<blockquote type="cite" cite="mid:12d9be6f-1aab-5773-f73d-a00f9106be45@cea.fr" class=""><p class="">Dear Barry,</p><p class="">I come back to you about this topic. As I
wrote last time, this is not a "highly urgent"
subject (whereas we will have to deal with it in the
next months), but it is an important one (especially
since the code I am working on just raised
significantly its ambitions). So I just would like
to check with you that this is still 'scheduled' on
your side.</p><p class="">I am sorry, I feel a little embarrassed
about asking you again about your work schedule, but
I need some kind of 'visibility' about this topic
which will become critical for our application.</p><p class="">Many thanks for helping me on that!<br class="">
Olivier<br class="">
</p>
<div class="moz-cite-prefix">On 02/12/2020 21:34,
Barry Smith wrote:<br class="">
</div>
<blockquote type="cite" cite="mid:ACDCA1CC-69A3-4068-A6B6-567F3A15A4DE@petsc.dev" class="">
<div class=""><br class="">
</div>
Sorry I have not gotten back to you quicker, give
me a few days to see how viable it is.
<div class=""><br class="">
</div>
<div class=""> Barry</div>
<div class=""><br class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Nov 25, 2020, at 11:57 AM,
Olivier Jamond <<a href="mailto:olivier.jamond@cea.fr" class="" moz-do-not-send="true">olivier.jamond@cea.fr</a>>
wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div class=""><p class="">Hi Barry,</p><p class="">I come back to you about the
feature to unlock the Ainsworth method
for saddle point problems in parallel.
If I may ask (sorry for that...), is it
still on your schedule (I just checked
the branch, and it seems 'stuck')?</p><p class="">This is not "highly urgent" on
my side, but the ability to handle
efficiently saddle point problems with
iterative solvers will be a critical
point for the software I am working
on...</p><p class="">Many thanks (and sorry again
for asking about your work schedule...)!<br class="">
Olivier<br class="">
</p>
<div class="moz-cite-prefix">On 12/10/2020
16:49, Barry Smith wrote:<br class="">
</div>
<blockquote type="cite" cite="mid:DF5D8D7C-8DD9-4E1A-8617-36A5FF472FA3@petsc.dev" class=""> <br class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Oct 12, 2020, at
6:10 AM, Olivier Jamond <<a href="mailto:olivier.jamond@cea.fr" class="" moz-do-not-send="true">olivier.jamond@cea.fr</a>>
wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div class=""><p class="">Hi Barry,</p><p class="">Thanks for this
work! I tried this branch with
my code and sequential
matrices on a small case: it
does work!</p>
<div class=""><br class="">
</div>
</div>
</div>
</blockquote>
Excellant. I will extend it to the
parallel case and get it into our
master release. </div>
<div class=""><br class="">
</div>
<div class=""> We'd be interested in
hearing about your convergence and
timing experiences when you run
largish jobs (even sequentially) since
this type of problem comes up
relatively frequently and we do need a
variety of solvers that can handle it
while currently we do not have great
tools for it.</div>
<div class=""><br class="">
</div>
<div class=""> Barry</div>
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">
<div class=""><p class="">Thanks a lot,<br class="">
Olivier<br class="">
</p>
<div class="moz-cite-prefix">On
09/10/2020 03:50, Barry Smith
wrote:<br class="">
</div>
<blockquote type="cite" cite="mid:218E7696-2A50-42A3-8CF2-D58FCC17B855@petsc.dev" class="">
<div class=""><br class="">
</div>
Olivier,
<div class=""><br class="">
</div>
<div class=""> The branch <b style="color: rgb(200, 20,
201); font-family: Menlo;
font-size: 14px;" class="">barry/2020-10-08/invert-block-diagonal-aij</b> contains
an example
src/mat/tests/ex178.c that
shows how to compute
inv(CC'). It works for
SeqAIJ matrices.</div>
<div class=""><br class="">
</div>
<div class=""> Please let
us know if it works for you
and then I will implement
the parallel version.</div>
<div class=""><br class="">
</div>
<div class=""> Barry</div>
<div class=""><br class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Oct 8,
2020, at 3:59 PM,
Barry Smith <<a href="mailto:bsmith@petsc.dev" class="" moz-do-not-send="true">bsmith@petsc.dev</a>>
wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div class=""><br class="">
Olivier<br class="">
<br class="">
I am working on
extending the
routines now and
hopefully push a
branch you can try
fairly soon.<br class="">
<br class="">
Barry<br class="">
<br class="">
<br class="">
<blockquote type="cite" class="">On Oct 8,
2020, at 3:07 PM,
Jed Brown <<a href="mailto:jed@jedbrown.org" class="" moz-do-not-send="true">jed@jedbrown.org</a>>
wrote:<br class="">
<br class="">
Olivier Jamond
<<a href="mailto:olivier.jamond@cea.fr" class="" moz-do-not-send="true">olivier.jamond@cea.fr</a>>
writes:<br class="">
<br class="">
<blockquote type="cite" class="">
<blockquote type="cite" class="">
Given the
structure of C
it seems you
should just
explicitly
construct Sp
and use GAMG
(or other
preconditioners,
even a direct
solver)
directly on
Sp. Trying to
avoid
explicitly
forming Sp
will give you
a much slower
performing
solving for
what benefit?
If C was just
some generic
monster than
forming Sp
might be
unrealistic
but in your
case CCt is is
block diagonal
with tiny
blocks which
means
(C*Ct)^(-1) is
block diagonal
with tiny
blocks (the
blocks are the
inverses of
the blocks of
(C*Ct)).<br class="">
<br class="">
Sp = Ct*C
+ Qt * S * Q
= Ct*C + [I
- Ct *
(C*Ct)^(-1)*C]
S [I - Ct *
(C*Ct)^(-1)*C]<br class="">
<br class="">
[Ct *
(C*Ct)^(-1)*C]
will again be
block diagonal
with slightly
larger blocks.<br class="">
<br class="">
You can do D =
(C*Ct) with
MatMatMult()
then write
custom code
that zips
through the
diagonal
blocks of D
inverting all
of them to get
iD then use
MatPtAP
applied to C
and iD to get
Ct *
(C*Ct)^(-1)*C
then
MatShift() to
include the I
then MatPtAP
or MatRAR to
get [I - Ct *
(C*Ct)^(-1)*C]
S [I - Ct *
(C*Ct)^(-1)*C]
then finally
MatAXPY() to
get Sp. The
complexity of
each of the
Mat operations
is very low
because of the
absurdly
simple
structure of C
and its
descendants.
You might
even be able
to just use
MUMPS to give
you the
explicit
inv(C*Ct)
without
writing custom
code to get
iD.<br class="">
</blockquote>
<br class="">
At this time, I
didn't manage to
compute
iD=inv(C*Ct)
without using <br class="">
dense matrices,
what may be a
shame because
all matrices are
sparse . Is <br class="">
it possible?<br class="">
<br class="">
And I get no
idea of how to
write code to
manually zip
through the <br class="">
diagonal blocks
of D to invert
them...<br class="">
</blockquote>
<br class="">
You could use
MatInvertVariableBlockDiagonal(),
which should
perhaps return a
Mat instead of a
raw array.<br class="">
<br class="">
If you have
constant block
sizes,
MatInvertBlockDiagonalMat
will return a Mat.<br class="">
</blockquote>
<br class="">
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</blockquote>
</div>
</div>
</blockquote>
</div>
<br class="">
</blockquote>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</blockquote>
</blockquote>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</blockquote>
</div>
</div></blockquote></div><br class=""></body></html>