[bddcml-users] Complex number support

Mon Jan 25 16:40:44 CET 2021

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> 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> 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 listbddcml-users at math.cas.czhttps://list.math.cas.cz/listinfo/bddcml-users
>>
>>
>>
>> --
>> Jakub Sistek, Ph.D.
>>
>> Researcher
>> Institute of Mathematics
>> Czech Academy of Sciences
>> sistek at math.cas.czhttp://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.czhttp://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/cf118c55/attachment-0001.html>