General floating point type support

[v-dobrev]

Support for floating point types other than double is going to be introduced in MFEM v4.7 with PR #3922.

Summary

This is a summary of the main changes in PR #3922:

• The type double is replaced throughout most of the library with a new type, real_t, which by default is just a typedef to double.
• MFEM can be built with double or single (float type) precision only. Other types may be added in the future.
• When MFEM is built with single precision, real_t is a typedef to float.
• The desired floating-point precision can be specified at configuration time with the option MFEM_PRECISION with one of the strings: single/Single/SINGLE/double/Double/DOUBLE.
  • In the configuration header config/_config.hpp, MFEM_PRECISION is represented by defining exactly one of the macros: MFEM_USE_DOUBLE, or MFEM_USE_SINGLE.
  • In the exported config files config.mk and MFEMConfig.cmake, the option is represented by the variables MFEM_USE_DOUBLE and MFEM_USE_SINGLE defined as YES/NO (make) or ON/OFF (cmake).
• As a convenience, MFEM defines the user-defined literal (https://en.cppreference.com/w/cpp/language/user_literal) operators real_t operator""_r() for integer and floating-point arguments.
  • This allows us to write real_t constants as 1_r, 0.25_r, etc.
  • A more general, standard way to ensure a constant or expression is converted to real_t type is to use c-style casting, e.g. real_t(1), real_t(0.25), real_t(M_PI*M_PI), reat_t(alpha(beta)), where alpha has the signature double alpha(double), etc.
• All examples, except ex33/ex33p, were updated and support single precision, on both CPUs and GPUs.
• Most miniapps were updated and support single precision, on both CPUs and GPUs.
• Most unit tests are NOT supported with single precision yet.
• Only some external libraries are supported with single precision:
  • HYPRE, LAPACK, MUMPS (Generalize floating point type #3922)
    • Note that HYPRE needs to be built with --enable-single when using MFEM_PRECISION=single.
  • PETSc (Changes to support PETSc 3.21 (to be released March 29th) and single precision builds #4203).
    • Note that PETSc needs to be built with --with-precision=single when using MFEM_PRECISION=single.
  • Pardiso, CPardiso (Fix Pardiso for single precision #4225)

Discussion

Replacing double with real_t in user code is NOT necessary, as long as MFEM is built with double precision -- in this case real_t is just an alias of double, i.e. they are the same type, just different names.

The only reason to convert double to real_t in user code is to add support for single precision to the user code -- and that will be just the first step. Depending on the application, adding support for single precision may require changes like: converting double constants to real_t constants to prevent undesirable computations in higher precision; or increasing tolerances (when MFEM_USE_SINGLE is defined) in iterative methods to ensure convergence is possible in the lower precision, see e.g. this code in the updated ex10:

mfem/examples/ex10.cpp

Lines 462 to 472 in 6f37514

	#if defined(MFEM_USE_DOUBLE)
	const real_t rel_tol = 1e-8;
	const real_t newton_abs_tol = 0.0;
	#elif defined(MFEM_USE_SINGLE)
	const real_t rel_tol = 1e-3;
	const real_t newton_abs_tol = 1e-4;
	#else
	#error "Only single and double precision are supported!"
	const real_t rel_tol = real_t(1);
	const real_t newton_abs_tol = real_t(0);
	#endif

[mlstowell]

I just ran into a compile time issue building with both MFEM_PRECISION = single and MFEM_USE_GSLIB = YES. There were several errors in fem/gslib.cpp of the form:
candidate function not viable: no known conversion from 'double *' to 'const mfem::real_t *' (aka 'const float *') for 1st argument

[kmittal2]

How about we add something like this to the makefile for now:

ifeq ($(MFEM_USE_SINGLE)$(MFEM_USE_GSLIB),YESYES)
   $(error GSLIB supported with double precision only)
endif

[v-dobrev]

If we add a check, it is better to do it for all libraries that are not supported with single precision. For example, we can add something like this in the makefile:

PKGS_NO_SINGLE = GSLIB ...
ifeq ($(MFEM_USE_SINGLE),YES)
   $(foreach pkg,$(PKGS_NO_SINGLE),$(if $(MFEM_USE_$(pkg):NO=),\
     $(error Package $(pkg) is NOT supported with single precision)))
endif

[tzanio]

We can maybe add this in the mfem-4.7-dev branch?

[tzanio]

@v-dobrev, does 18cf9d7 look OK?

[v-dobrev]

Looks good to me. Thanks!

[IdoAkkerman]

Any ideas whether the real_t could be substituted for a CodiPack alias, so automatic differentiation would be available?
I assume linear solves would be next to impossible, however derivatives of linear and bilinear forms might be possible.
Although computationally costly, this might be interesting form a academic and/or debug perspective.
Any thoughts?

[v-dobrev]

@bslazarov probably knows best if this is possible.

[IdoAkkerman]

I gave Codipack it a try and after some hacking got it to compile. I can run simple examples/miniapps - no derivatives yet.
There are basically three of issues:

1. There are still a couple of implicit conversions in the code. When using double of single the compiler takes care of this. For more complicated types, such as the codi alias this does not work. I can probably do a low-impact PR of this.

There is probably also a desire for explicit functions from real_t to double and int, and/or vice-versa. This might streamline interfacing with the "outside world", namely: input, output (datacollection for instance) and perhaps even external libraries.

2. The codi alias is not a constant so most/all constexpr keywords gave errors. For now I simply got rid of the keyword.
   A potential solution for this would be to use some kind of MACRO instead of the literal constexpr keywords. This might be considered to obfuscating or invasive.

3. The code alias is not a Trivial type. This gives issues in the array class (and perhaps in 1 or 2 other locations). I side stepped this by simple commenting out the static_assert. I do not know the implications of this... is the assertion necessary?

There was also some small by-catch: in-advert conversions from int to real_t in 2 locations I believe.

[v-dobrev]

Just to clarify, are you building MFEM with MFEM_USE_CODIPACK=YES (and MFEM_PRECISION=single or double) and trying to build and test the miniapps in miniapps/autodiff? I think this is the only place in MFEM where CoDiPack can be used at the moment.

3. The code alias is not a Trivial type. This gives issues in the array class (and perhaps in 1 or 2 other locations). I side stepped this by simple commenting out the static_assert. I do not know the implications of this... is the assertion necessary?

In class Array, data is copied (e.g. when resizing the array) using a simple memcpy which is only legal with trivial types. The same is true for class Memory where we also copy data between host and device memory.

[IdoAkkerman]

Indeed I use MFEM_USE_CODIPACK=YES, but in config/config.hpp I added the option to define using real_t = code::RealForward or using real_t = code::RealReverse depending on the additional codipack flag.

I got all kinds of compile errors letting me to believe I am actually using a codipack wrapper instead of a naked double or real.
This also lead me to find the error corrected in PR #4265.

[jandrej]

Hi @IdoAkkerman using real_t = code::RealForward is quite excessive. This will try to differentiate through MPI calls and other third party libraries which is probably not the intended usage.