Mixed precision functionality

[nbeams]

This issue is for discussions related to adding mixed precision capabilities to libCEED, including:

• Which user-facing functions should have mixed precision capability?
• How should the libCEED interface support this (what will the user need to do)?
• What combinations should be available (e.g. mix-and-match input and output precisions, ability to decouple input/output precision from precision of computations, ...)
• Potential for backend specializations, and how this will fit with the interface changes

An idea from @jedbrown was to handle precision conversions through the CeedVector objects, which could be modified to store multiple pointers for different precisions, along with a tag to indicate which precision that vector's data currently uses.

[jedbrown]

I think qfunction variation is the hardest part. Do we need to support different arguments to the qfunction having different precisions or can they be homogeneous (all arguments in matching precision)? Let's defer this part for later.

As a first step, I would add functions like CeedVectorGetArray_f32(CeedVector, CeedMemType, float **) and have the implementation convert as needed. We can use C11 _Generic so C11 users don't need to specify these types and C++ overloading to do the same for C++ users.

We can have a macro CEED_SCALAR_F32 that users can define before including ceed.h so that C99 callers get suitable CeedScalar and appropriate specialization. Perhaps better is to provide #include <ceed/ceed_f32.h> that sets the macro and includes ceed/ceed.h). This would be for callers who want to choose precision at compile time.

We then add a CeedElemRestrictionSetScalarType(CeedElemRestriction, CeedScalarType evec_type) that converts if necessary (when provided with an L-vector). This will require visiting some backends for the cross-type packing, but I don't think will be a lot of code.

I kinda think CeedBasis public interfaces (setting matrices and doing linear algebra) could just stay double precision all the time, and find it more compelling to possibly add __float128 than to support direct specification in single precision (they can be converted as needed behind the scenes).

We can have CeedVectorGetArrayUntyped(CeedVector, CeedMemType, CeedScalarType, void *) to make life easier for dynamic callers who prefer to cast somewhere else.

[nbeams]

I think qfunction variation is the hardest part. Do we need to support different arguments to the qfunction having different precisions or can they be homogeneous (all arguments in matching precision)? Let's defer this part for later.

Just to make sure I'm clear: are you referring to the user-provided per-point QFunction, the libCEED QFunction interface, or both?

Regarding CeedBasis, when users can tolerate lower precision for their application(s), they could potentially get some nice speedup from single precision CeedBasisApply in the non-tensor case. Specifically, for the MAGMA backend, with the GEMM/batch GEMM computations, the ability to use SGEMM instead of DGEMM could be useful. I don't know if this is considered too narrow of a use case (non-tensor basis and no need for double-precision accuracy in the basis computations)...also, I guess CeedBasisApply is different than the basis creation functions since it doesn't take any CeedScalar arguments, only CeedVectors. If CeedVector knows about its type, a lot of the "overloading" could be handled by the backends inside the backend-specific implementations. I guess this goes back to the question of whether we want to be able to specify computation precision, or have it deduced from the types of the CeedVectors? I'd argue there would likely be interest in being able to do computations in double while input/output are single, but also to be able to do everything in single. Edit: was the intention that _f[something] would be added to routines to specify computation precision, but not necessarily input/output?

[jedbrown]

CeedBasis has several interfaces that take CeedScalar *, and I think there's no harm in making those all double (perhaps with complex variants at some point). It can be converted internally. We have several right now

1. Outer solver/Krylov vector storage (PETSc, MFEM, etc.) which gets placed into a CeedVector as an L-vector.
2. E-vector representation (ephemeral)
3. CeedBasisApply
4. Q-vectors processed by user-provided QFunction

We could make CeedElemRestriction have no opinion about types, in which case we would convert as needed. Then L-vectors, CeedBasis, and CeedQFunction would have a designated precision. The precision of E-vectors would be inferred to match the compute precision of CeedBasis (converting as needed to/from the L-vector precision).

A simplification would be to make CeedBasis precision always match the QFunction precision.

To make this incremental, I think it will be simpler to start with the L-vectors and get that merged, then work inward.

[nbeams]

I'm not sure I understand what your list is intended to represent. Cases where we use CeedScalar * directly as a type? CeedBasisApply only takes CeedVectors for input/output. E-Vectors are also represented as CeedVectors. So I must be missing something.

I like the idea of the element restriction being able to infer the output (or input, in transpose mode) type, but I'm not sure how it would work in practice. The CeedElemRestriction doesn't know about the CeedBasis or CeedQFunction, so it can't know which precision they use unless we tell it. I think your initial suggestion of a way to set the precision of the element restriction may be necessary to keep the current separation between/independence of CeedOperator "pieces" (element restriction, basis, qfunction).

[jedbrown]

1. The list included internal representations (E-vectors), and I intended it to identify each place where precision can be varied in G^T B^T D B G. I think it's a significant simplification in the code if we can make the precision of B match the precision of the Q-function D -- it still decouples L-vector precision from compute precision. Note that each (active or passive, input or output) field provided to a CeedOperator could have a different precision. This does restrict algorithmic space somewhat, but I don't know how serious it is. On CPUs, I think it's not a big deal because E-vectors only ever exist in cache. It may be a lot more significant for GPU kernels that need to store complete E-vectors.

2. Backends can always fuse CeedElemRestriction into whatever else they do. The CPU backends, for example, never have an actual E-vector, just a bit of workspace for one element (/serial) or a batch (/blocked). That said, we have two vectors, so if those vectors have preferred precision set, then all the necessary information is there.

int CeedElemRestrictionApply(CeedElemRestriction rstr,
    CeedTransposeMode t_mode, CeedVector u, CeedVector ru, CeedRequest *request);

The CeedOperator backend gets to choose how to represent the E-vector data so even if it just calls CeedElemRestrictionApply, that function knows what precisions to work with.

[nbeams]

Right, I guess this raises another important question: how much do we want to limit the multiprecision functionality to things that happen inside a CeedOperatorApply? I guess if CeedVector provides a function for converting to other precisions, that would cover any cases where a user (or library) wants to do a stand-alone CeedElemRestrictionApply for something and have the output in a different precision -- should such a situation arise. Or we could specify that you should set the precision of the input/output vectors ahead of time, as suggested.

Though even if we are only considering the case where the Operator is setting the preferred element restriction output precision based on all the information it has, the function that lets you set the element restriction precision might still make more sense than relying on the input/output CeedVectors having their precisions set. I'm thinking this would be more consistent across backends, including those that never set up E-Vectors as CeedVector objects, like cuda-gen -- which does still get information from the element restrictions while building the code string, despite not creating E-Vectors (or using CeedElemRestrictionApply).

[jedbrown]

I figured CeedVectorGetArray() and the like would be able to mirror data to any requested precision, but CeedElemRestrictionApply would inspect the available precisions of the input vector and the preferred precision of the output vector and apply without needing to convert the input.

I'd think cuda-gen could still propagate the desired (ephemeral) E-vector precisions from the Basis. One thing I'd like to facilitate is having a residual in double precision and Jacobian action in single, or the Jacobian in double and a preconditioner in single.

[nbeams]

I'd think cuda-gen could still propagate the desired (ephemeral) E-vector precisions from the Basis.

Yes, it could, I just thought it might be preferable to have more consistency across the CUDA/HIP backends, such that they could all get/set this information in the same way. But cuda(hip)-gen is so different from the other backends already that perhaps that's not really a selling point.

[jedbrown]

Do you think it'll be hard to handle consistently if CeedElemRestriction is precision-agnostic (uses its arguments rather than set as state)?

My concern is that when we wire up several different operators, we'd like for the data (integers) of a CeedElemRestriction to be shared between different precisions. That's easier if it's just the same object, rather than needing to reference count the backing storage.