<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<div class="moz-cite-prefix">Hi Sebastian,<br>
<br>
I have added you as a developer of BDDCML. Just do not push to the
master, please. Another branch, say "complex", is fine!<br>
<br>
I will add Juan Calvo into the discussion on Nedelec elements vs.
BDDC(ML). He has studied these things. I have never worked with
these elements too closely really, but would be happy to learn.<br>
<br>
I think that I now understand your question about the inner
products and weights! It is a slight abuse of the concept of
partition of unity. In DD methods, the weighting operator plays a
crucial role and there is a number of strategies available like
the diagonal stiffness etc. Do you actually plan to use the
diagonal stiffness at this side?<br>
<br>
In the Krylov method, on the other hand, the operator is just
reused to get the right outcome from the duplicated entries at the
interface. There are multiple ways how to do this, but I decided
to go for this weighted inner products to avoid the need to assign
unknowns to processes in a unique way as this is not needed
anywhere else in the code. However, ANY partition of unity would
do this job without affecting the performance of the
preconditioner. So, the simple weights by cardinality would help
here just fine. If you plan to use the diagonal stiffness in the
preconditioner, then we will perhaps need to use separate weights
here, it is not a problem although we would need to think how to
do that. Otherwise, you could just use the<br>
wights_type = 0<br>
straightaway.<br>
Even if you needed the diagonal stiffness weights and thus the new
partition of unity for inner products, that should be a rather
simple thing to be done. I think I would not put the complex
numbers as weights into these inner products.<br>
<br>
Let us keep in touch on this! I will contact the other developers
to get their points of view and put you in the loop.<br>
<br>
Do not hesitate to ask! The code can be quite hard to read at many
places so a quick question can be worth hours of decyphering the
code. I am more than happy to help!<br>
<br>
Best wishes,<br>
<br>
Jakub<br>
<br>
<br>
<br>
<br>
<br>
On 25/01/2021 16:40, Sebastian Grimberg wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CANX_vp=+HH=paQNF+1co+2LVaMZ24m9Y-hTtiZwMy4tXbuKLrQ@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<div dir="auto">Hi Jakub,</div>
<div dir="auto"><br>
</div>
<div dir="auto">That’s good to know about the Nedelec element
implementation. There is a paper on edge elements for FETI by
Tosseli that was the precursor for the H(curl) BDDC
implementations it seems, which suggested a third alternative
algorithm in which all the subdomain edge degrees of freedom are
included in the primal basis rather than just corners (again my
understanding is not complete but this was the basic grasp I
got). This removes the need for a change of basis at the expense
of a slightly larger coarse space. I’m hoping that I may be able
to test with this approach in BDDCML and see if I can make any
progress. It seems like this is similar to the default coarse
space that BDDCML uses with averages over the edges included as
constraints, but correct me if I'm wrong.</div>
<div dir="auto"><br>
</div>
<div dir="auto">Otherwise, I agree that this sounds like a good
plan. My GitHub username is sebastiangrimberg. I have a local
branch called “complex” so could push the changes to the repo
after I’ve tested them.</div>
<div dir="auto"><br>
</div>
<div dir="auto">For the other stuff, I use GMRES as a Krylov
solver usually but would expect similar performance from
BiCGSTAB which is already implemented. I agree that changing the
inner products to just be the standard complex inner product
should work fine, that is where my question about the weights
came from because I wasn’t sure where the complex conjugate
would go if they were complex and not Hermitian when computing
inner products.</div>
<div dir="auto"><br>
</div>
<div dir="auto">I agree that this is a good idea where I can get
started with the complex stuff and hopefully you can get some
info from other project developers for the best approach to
handle multiple types. I’ll keep you posted how things go on my
end!</div>
<div dir="auto"><br>
</div>
<div dir="auto">Thanks again!</div>
<div dir="auto"><br>
</div>
<div dir="auto">Sebastian </div>
<div dir="auto"><br>
</div>
<div><br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Mon, Jan 25, 2021 at
12:28 AM Jakub Sistek <<a
href="mailto:sistek@math.cas.cz" moz-do-not-send="true">sistek@math.cas.cz</a>>
wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-style:solid;padding-left:1ex;border-left-color:rgb(204,204,204)">
<div>
<div>Hi Sebastian,<br>
<br>
nice to meet you as well! :-)<br>
<br>
I will react to your points in the text below...<br>
<br>
On 22/01/2021 23:47, Sebastian Grimberg wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">Hi Jakub,
<div><br>
</div>
<div>Thanks for getting back to me, it's nice to
meet you!</div>
<div><br>
</div>
<div>The applications I'm interested in
involve high-frequency time-harmonic
electromagnetics, so typically dealing with
complex symmetric indefinite operators. I know
there is quite a bit of H(curl) specific BDDC
literature out there involving better selection of
coarse constraints using change of basis for
Nedelec elements, but I figure if I end up needing
this I could use the general user-defined
constraints in BDDCML (maybe you have some
thoughts on this?).</div>
<div><br>
</div>
</div>
</div>
</blockquote>
A few years ago, we though of supporting Nedelec elements
for Maxwell's equations together with Juan Calvo, who was
a postdoc at the Institute of Mathematics in 2017.
Unfortunately, this was when I was on leave, so we could
only collaborate via Skype and email, and we did not
finish with that. Juan is an expert on the BDDC literature
in this regard. However, these approaches rely on changing
the basis so that the coarse unknows appear as primal
unknowns. I do not remember the details, but there were
some limits to what one can do through the user-defined
constraints.<br>
<br>
The BDDCML does not use the change of basis at the moment,
and in my opinion, it would be a major change if it
should... I may be wrong. Consequently, it also does not
support the so called deluxe scaling which relies on it.<br>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">
<div>I'd be happy to work on implementing complex
number support in the library and see where it
goes. I have a few thoughts:</div>
<div><br>
</div>
<div>1. At a high level, I could see two ways to
manage the different scalar types.</div>
<div><br>
</div>
</div>
<blockquote style="margin:0px 0px 0px
40px;border:none;padding:0px">
<div dir="ltr">
<div>- First, something like what MUMPS does, is
to just maintain two versions of every
type-specific file so the user would include the
correct module for double/complex support. This
could be simple for just double/complex, but if
later you also want to add an option for single
precision or 64-bit integers it becomes rather
involved.<br>
</div>
<div>- Alternatively, we could use conditional
compilation in a more "C" approach to define at
compile time the scalar type to be real(kind=8)
or complex(kind=8) and use this throughout the
code, wrapping type specific commands in #if
defined/#else pragmas. I'm not even sure if this
is possible in Fortran but would be the way to
do it in C/C++ I think. Do you know anything
about this?</div>
</div>
</blockquote>
</div>
</blockquote>
Yes, this is a major decision we should discuss properly
before we start. I have not extended many projects in this
way, but I have worked on the PLASMA library, which
supports four precisions (double real, single real, double
complex, and single complex) and uses a third way, a
textual substitution done by a Python script. I think that
the MUMPS source codes are a product of something like
that as well. Only one version of the code is developed
(complex double) and the others are generated and compiled
automatically. The drawback of this approach would be that
you compile a different code than what you actually see,
which makes searching for bugs one step harder. And a new
type of strange bugs as well.<br>
<br>
The ifdefs are possible in Fortran as well. When one uses
capital F or F90 in the suffix, the file is processed by a
C preprocessor before a compilation. Yes, what you suggest
makes sense, with the substitution done either by a
preprocessor or the generating script. However, the script
could also rename all the functions in a semi-automatic
way.<br>
<br>
I think that we may end up with a combination of the two
approaches. In particular because certain operations make
sense only in complex and not in double. So a define line
COMPLEX may still be useful, however just with one branch.
This is what PLASMA is doing. But it is a C code.<br>
<br>
In Fortran, things like Parametrized types could be also
useful. This seems to be a mechanism to bring some of the
C++ templates into Fortran. <br>
<a
href="https://blog.hpc.qmul.ac.uk/fortran-parameterized-derived-types-1.html"
target="_blank" moz-do-not-send="true">https://blog.hpc.qmul.ac.uk/fortran-parameterized-derived-types-1.html</a><br>
I am not sure if we could use this. It may mean rewriting
all the array oriented code, which I would like to avoid,
not sure.<br>
<br>
<blockquote type="cite">
<div dir="ltr">
<blockquote style="margin:0px 0px 0px
40px;border:none;padding:0px">
<div dir="ltr">
<div> </div>
</div>
</blockquote>
<div dir="ltr">
<div><br>
</div>
<div>2. The blopex object files for double versus
complex share the same names (see the files
in blopex_serial_double or blopex_serial_complex
which would be needed). I presume this would cause
some build issues, which further motivates
building only a single library for real OR complex
numbers and then just linking to the correct
blopex version. But again, I don't know how this
is possible in Fortran.</div>
</div>
</div>
</blockquote>
Yes, this is a very good point! Perhaps we should not
combine all the versions into just one big library but
rather make just one version with the data types one
wants.<br>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">
<div><br>
</div>
<div>3. In the case of complex matrices, the
partition of unity weights for interface degrees
of freedom could be complex (when using diagonal
stiffness scaling, for example). Is this allowed?
Do we require them to be Hermitian at least? I see
that for indefinite problems they are constrained
to be positive. This is a bit beyond my expertise
but I am not sure of the implications for
requiring them to always be real and positive or
allowing them to be complex and how to handle that
in the dot products required for Krylov methods.
I'm not sure what PETSc's PCBDDC does about this.</div>
</div>
</div>
</blockquote>
Well, this would probably require studying the DD
literature devoted to these problems. PETSc will probably
use the deluxe scaling. I would say that algorithmically,
the weights with a diagonal are very general and could be
most likely generalized to complex numbers. I do not see
an issue here at the moment, but may be wrong. One can
usually be quite creative at this point.<br>
<br>
To summarize, my idea would be to start simply with a new
branch and devising just the version you need, using
complex numbers. Once we would make sure it is doable and
the code can work, we would convert it to the generic
version. It is difficult to develop all the versions at
once.<br>
<br>
In the meantime, we can discuss with the developers of the
other related codes about what they do for devising the
different precision versions.<br>
<br>
To start with, what is your GitHub name? I will create a
branch for this purpose and add you as a developer.<br>
<br>
It would be good to have also a simple benchmark problem
for testing the methods. Can you think of something?<br>
<br>
I would probably start with changing the number type
somewhere high (like module_bddcml) and followed the
compiler errors with respect to types. This should lead to
the necessary changes. There is some dead code in
module_sm, which will be otherwise full of changes, which
we do not need to even look at in this approach.<br>
It will be quite tedious work before we change that on all
the required places, but it will not require us to look at
all the code.<br>
<br>
I also do not have much experience with Krylov methods
with respect to complex numbers, but I know methods such
as CG work for Hermitian problems. What methods are used
for your class of problems? In scalar products, it should
be enough to consider complex conjugation in transpose,
but I have never tried to implement it actually... so
there will be a lot of new stuff for both of us!<br>
<br>
One last point to get something working quickly. If one
uses only one process, the code falls back to just calling
MUMPS. The same if one does "just_direct_solve" which just
takes the loaded data and passes it to the MUMPS solver.
This may be a shortcut for having something working before
we process the whole code in this way...<br>
<br>
I am excited about this project and I am looking forward
to our collaboration on this! Let us hope we will be able
to make it!<br>
<br>
Best wishes,<br>
<br>
Jakub</div>
<div><br>
<br>
<br>
<blockquote type="cite">
<div dir="ltr">
<div dir="ltr">
<div><br>
</div>
<div>Sorry for the long email, I hope this makes
sense to you and I'm curious to know your input
before getting started. Thanks again!</div>
<div><br>
</div>
<div>Sebastian</div>
</div>
</div>
</blockquote>
<br>
<blockquote type="cite">
<div dir="ltr"><br>
<div class="gmail_quote">
<div dir="ltr" class="gmail_attr">On Fri, Jan 22,
2021 at 7:28 AM Jakub Sistek <<a
href="mailto:sistek@math.cas.cz" target="_blank"
moz-do-not-send="true">sistek@math.cas.cz</a>>
wrote:<br>
</div>
<blockquote class="gmail_quote" style="margin:0px
0px 0px
0.8ex;border-left-width:1px;border-left-style:solid;padding-left:1ex;border-left-color:rgb(204,204,204)">
<div>
<div>Hi Sebastian,<br>
<br>
first of all, thank you for your interest in
BDDCML!<br>
<br>
So far, there has been no need for supporting
complex numbers in BDDCML, so you are the
first one asking for this :-) What
applications do you have in mind?<br>
<br>
Implementation-wise, this may not be too bad,
as you suggest, although the "double" type is
hard-coded in a lot of subroutines where it
would require changing. Even blopex has a
version for complex numbers, although one
should start without the adaptivity. And yes,
it would be applicable only to Hermitian
matrices. MUMPS supports all the types one
needs, so this should be straightforward.<br>
<br>
More flexible variable types would be helpful
in other scenarios I am interested in, such as<br>
* using long integers for certain indices for
very large problems<br>
* using single precision instead of double<br>
<br>
Also algorithmically, I am not aware of any
principal issues. I know other DD codes
support complex numbers, such as HPDDM or the
PCBDDC preconditioner in PETSc. I know the
developers, so I would be able to ask for
advice.<br>
<br>
This would actually be a nice extension of the
functionality of BDDCML! If you were
interested in helping with this, I would be
happy to give you access to the project on
Github and we could try to devise such version
in a new branch together.<br>
<br>
Let me know what you think.<br>
<br>
Best wishes,<br>
<br>
Jakub<br>
<br>
<br>
<br>
On 22/01/2021 05:29, Sebastian Grimberg wrote:<br>
</div>
<blockquote type="cite">
<div><span
style="word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)">Hello
everyone,</span>
<div dir="auto"
style="word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)"><br>
</div>
<div dir="auto"
style="font-size:1rem;word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)">I’m
a new user and am wondering if there is
any interest or thoughts on supporting
complex numbers in BDDCML? I’m unsure
about the blopex dependency for
adaptivity, and further it may not even be
applicable unless the input materials is
Hermitian, but it seems like it could be
doable by calling the correct MUMPS
library and adjusting the interface a bit.
Would this be worth spending time to
develop, and are there any other
foreseeable issues?</div>
<div dir="auto"
style="word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)"><br>
</div>
<div dir="auto"
style="font-size:1rem;word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)">Cheers,</div>
<div dir="auto"
style="word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)"><br>
</div>
<div dir="auto"
style="font-size:1rem;word-spacing:1px;border-color:rgb(49,49,49);color:rgb(49,49,49)">Sebastian
Grimberg </div>
</div>
<br>
<fieldset></fieldset>
<pre style="font-family:monospace">_______________________________________________
bddcml-users mailing list
<a href="mailto:bddcml-users@math.cas.cz" target="_blank" style="font-family:monospace" moz-do-not-send="true">bddcml-users@math.cas.cz</a>
<a href="https://list.math.cas.cz/listinfo/bddcml-users" target="_blank" style="font-family:monospace" moz-do-not-send="true">https://list.math.cas.cz/listinfo/bddcml-users</a>
</pre>
</blockquote>
<br>
<br>
<pre cols="72" style="font-family:monospace">--
Jakub Sistek, Ph.D.
Researcher
Institute of Mathematics
Czech Academy of Sciences
<a href="mailto:sistek@math.cas.cz" target="_blank" style="font-family:monospace" moz-do-not-send="true">sistek@math.cas.cz</a>
<a href="http://www.math.cas.cz/~sistek" target="_blank" style="font-family:monospace" moz-do-not-send="true">http://www.math.cas.cz/~sistek</a>
+420 222 090 710</pre>
</div>
</blockquote>
</div>
</div>
</blockquote>
<br>
<br>
<pre cols="72" style="font-family:monospace">--
Jakub Sistek, Ph.D.
Researcher
Institute of Mathematics
Czech Academy of Sciences
<a href="mailto:sistek@math.cas.cz" target="_blank" style="font-family:monospace" moz-do-not-send="true">sistek@math.cas.cz</a>
<a href="http://www.math.cas.cz/~sistek" target="_blank" style="font-family:monospace" moz-do-not-send="true">http://www.math.cas.cz/~sistek</a>
+420 222 090 710</pre>
</div>
</blockquote>
</div>
</div>
</blockquote>
<br>
<br>
<pre class="moz-signature" cols="72">--
Jakub Sistek, Ph.D.
Researcher
Institute of Mathematics
Czech Academy of Sciences
<a class="moz-txt-link-abbreviated" href="mailto:sistek@math.cas.cz">sistek@math.cas.cz</a>
<a class="moz-txt-link-freetext" href="http://www.math.cas.cz/~sistek">http://www.math.cas.cz/~sistek</a>
+420 222 090 710</pre>
</body>
</html>