<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">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">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">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">bddcml-users@math.cas.cz</a>
<a href="https://list.math.cas.cz/listinfo/bddcml-users" target="_blank" style="font-family:monospace">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">sistek@math.cas.cz</a>
<a href="http://www.math.cas.cz/~sistek" target="_blank" style="font-family:monospace">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">sistek@math.cas.cz</a>
<a href="http://www.math.cas.cz/~sistek" target="_blank" style="font-family:monospace">http://www.math.cas.cz/~sistek</a>
+420 222 090 710</pre>
</div>
</blockquote></div></div>