[bddcml-users] Complex number support
Jakub Sistek
sistek at math.cas.cz
Mon Jan 25 22:55:04 CET 2021
Hi Sebastian,
I have added you as a developer of BDDCML. Just do not push to the
master, please. Another branch, say "complex", is fine!
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.
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?
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
wights_type = 0
straightaway.
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.
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.
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!
Best wishes,
Jakub
On 25/01/2021 16:40, Sebastian Grimberg wrote:
> Hi Jakub,
>
> 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.
>
> 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.
>
> 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.
>
> 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!
>
> Thanks again!
>
> Sebastian
>
>
> On Mon, Jan 25, 2021 at 12:28 AM Jakub Sistek <sistek at math.cas.cz
> <mailto:sistek at math.cas.cz>> wrote:
>
> Hi Sebastian,
>
> nice to meet you as well! :-)
>
> I will react to your points in the text below...
>
> On 22/01/2021 23:47, Sebastian Grimberg wrote:
>> Hi Jakub,
>>
>> Thanks for getting back to me, it's nice to meet you!
>>
>> 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?).
>>
> 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.
>
> 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.
>> I'd be happy to work on implementing complex number support in
>> the library and see where it goes. I have a few thoughts:
>>
>> 1. At a high level, I could see two ways to manage the different
>> scalar types.
>>
>> - 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.
>> - 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?
>>
> 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.
>
> 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.
>
> 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.
>
> 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.
> https://blog.hpc.qmul.ac.uk/fortran-parameterized-derived-types-1.html
> 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.
>
>>
>> 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.
> 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.
>>
>> 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.
> 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.
>
> 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.
>
> In the meantime, we can discuss with the developers of the other
> related codes about what they do for devising the different
> precision versions.
>
> To start with, what is your GitHub name? I will create a branch
> for this purpose and add you as a developer.
>
> It would be good to have also a simple benchmark problem for
> testing the methods. Can you think of something?
>
> 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.
> 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.
>
> 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!
>
> 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...
>
> 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!
>
> Best wishes,
>
> Jakub
>
>
>
>>
>> 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!
>>
>> Sebastian
>
>>
>> On Fri, Jan 22, 2021 at 7:28 AM Jakub Sistek <sistek at math.cas.cz
>> <mailto:sistek at math.cas.cz>> wrote:
>>
>> Hi Sebastian,
>>
>> first of all, thank you for your interest in BDDCML!
>>
>> 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?
>>
>> 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.
>>
>> More flexible variable types would be helpful in other
>> scenarios I am interested in, such as
>> * using long integers for certain indices for very large problems
>> * using single precision instead of double
>>
>> 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.
>>
>> 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.
>>
>> Let me know what you think.
>>
>> Best wishes,
>>
>> Jakub
>>
>>
>>
>> On 22/01/2021 05:29, Sebastian Grimberg wrote:
>>> Hello everyone,
>>>
>>> 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?
>>>
>>> Cheers,
>>>
>>> Sebastian Grimberg
>>>
>>> _______________________________________________
>>> bddcml-users mailing list
>>> bddcml-users at math.cas.cz <mailto:bddcml-users at math.cas.cz>
>>> https://list.math.cas.cz/listinfo/bddcml-users
>>
>>
>> --
>> Jakub Sistek, Ph.D.
>>
>> Researcher
>> Institute of Mathematics
>> Czech Academy of Sciences
>>
>> sistek at math.cas.cz <mailto:sistek at math.cas.cz>
>> http://www.math.cas.cz/~sistek
>> +420 222 090 710
>>
>
>
> --
> Jakub Sistek, Ph.D.
>
> Researcher
> Institute of Mathematics
> Czech Academy of Sciences
>
> sistek at math.cas.cz <mailto:sistek at math.cas.cz>
> http://www.math.cas.cz/~sistek
> +420 222 090 710
>
--
Jakub Sistek, Ph.D.
Researcher
Institute of Mathematics
Czech Academy of Sciences
sistek at math.cas.cz
http://www.math.cas.cz/~sistek
+420 222 090 710
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://list.math.cas.cz/pipermail/bddcml-users/attachments/20210125/744ccf7a/attachment-0001.html>
More information about the bddcml-users
mailing list