[bddcml-users] Complex number support

Jakub Sistek sistek at math.cas.cz
Mon Jan 25 09:28:16 CET 2021 

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
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/a6ed2127/attachment.html>

More information about the bddcml-users mailing list