Half factorization

[yhmtsai]

this pr adds the factorization with half support.

Hip does not support atomic on the 16bits type currently

TODO:

• add the fix of tri solve with half

[MarcelKoch]

Generally LGTM. I have a question regarding atomics and hip. The latest ROCm shows support for fp16 atomic operations: https://rocm.docs.amd.com/en/latest/reference/precision-support.html#atomic-operations-support, but TBH I can't figure out what operations exactly they mean with that. Did you try anything in that regard?

[MarcelKoch]

Many of the tests here are missing SKIP_HALF if compiling for HIP.

[yhmtsai]

we do not support compute_l_u_factors in hip, but the others still works with half precision in HIP

[yhmtsai]

I got your meaning now

[MarcelKoch]

nit:

Suggested change

	// it does not work in cuda110/100 images
	// Explicitly cast `nullptr` to `const ValueType*` to prevent compiler issues with cuda 10/11

[yhmtsai]

I think it is more on the host compiler side because it goes through our binding first with specfic type

[upsj]

nit:

Suggested change

	// automatically Explicitly cast `nullptr` to `const ValueType*` to
	// automatically explicitly cast `nullptr` to `const ValueType*` to

[upsj]

Could SharedValueType be deduced inside, instead of making it an additional template parameter? You should be able to pull the code from the kernel launch into here and add a type alias. Otherwise it is easier to accidentally call the kernel with inconsistent types.

[yhmtsai]

good idea

[upsj]

I'm not sure we should go this far to accommodate broken compilers. We have workarounds for compilation issues, but not really for this degree of broken-ness.

[upsj]

For HIP 16 bit atomics, as long as you only use load and store, you could implement them as

• a 32 bit load_* plus a memcpy and,
• a 32 bit load_* plus an atomic CAS, similar to how we did atomicAdd in the past. For safety, we need to execute this twice, assuming that every memory location only ever gets written once (which is true for all algorithms that use atomics), because your write can either fail because the upper half changed, or the lower half changed, one of which belongs to you and can't change without your knowledge.

[yhmtsai]

using 32 bit memory operation for 16 bit, it will have illegal memory access in the tail or head if we do not handle it in a upper level.

[upsj]

Theoretically that would be an easy fix: Make sure all allocations are at least 32 bits and rounded up to multiples of 4. But I believe most allocators already silently fulfill that assumption, and GPUs are unlikely to have 16 bit allocation boundaries for alignment purposes.

[yhmtsai]

I do not like slight guarantee unless we have a way to ensure or at least check.
However, I would suggest we do not consider it for this pr and release such that we have enough time ensure that it works correctly on hip.

[upsj]

I can give you a somewhat technical justification for this: cudaMalloc returns correctly aligned memory for thrust::complex<double>, despite not knowing anything about the type. So that means that the allocator is not using any space between those 16 byte-aligned allocations. Whether this is special-cased for allocations divisible by 16 or not I'm not sure (I would assume not, since people also allocate memory pools themselves), but again, we have an easy fix, which I would honestly consider useful in any case: round up the sizes raw_allow uses to at least be divisible by 4.

[yhmtsai]

I know the idea, sometimes it is necessary for optimized half precision by packing them (so, we will have kind of natively 32 bit by enforcing packing structure requirement)
I will still say it is not easy and confident to say it will be correct in this short period.
For example, user allocates some memory with 16 bits type but only pass the odd number to array_view.
should we accept or throw the error? Of course, these memory operation will not change the value out of the actual array, but it is still illegal memory operation.