Weaning off std::complex arithmetic

[TysonRayJones]

Experiment with fully custom cpu_qcomp and gpu_qcomp types to avoid...

• compiler-specific flags to overcome qcomp = std::complex performance pitfalls
• HIP-specific functions to overcome cu_qcomp correctness pitfalls

[TysonRayJones]

@JPRichings @otbrown Here's a first draft of the proposed qcomp changes (ignore that it provokes a stack-overflow bug in multithreaded MSVC, grr). One only really needs to look at cpu_types.hpp and gpu_types.cuh, and at an example of how cpu_qcomp is used within cpu_subroutines.cpp, like here.

You'll see most of the cpu_types.hpp and gpu_types.cuh boilerplate is identical, merely expressing the arithmetic operators of (c|g)pu_qcomp via their qreal components. I propose further simplification by defining a new complex.hpp (or similarly named, perhaps basetypes.hpp?) in core/ (rather than in cpu/ or gpu/) which defines all the operator overloads of a templated type, and (c|g)pu_qcomp gets defined in a typedef of that "ancestor type" in cpu/ or gpu/. Then adding a new operator overload to all backends is trivial (just add one function to complex.hpp), as is adding a backend specific function (define it exclusively within cpu_types.hpp or gpu_types.hpp). Noting this comment here in case you object to the unification (so I can see what to roll-back if I find time to implement beforehand - which is alas unlikely!)

[otbrown]

Thanks Tyson!

I have 5(!) hours of meetings tomorrow, so unlikely to get to this then, but will make time later in the week!

[otbrown]

Hi Tyson, one of my meetings got cancelled 🥳

Moving away from the standard library entirely is philosphically upsetting, but clearly has practical benefits!

I propose further simplification by defining a new complex.hpp

I strongly support this proposal. Structurally it would be nice to have the common interface in core to then be specialised/overwritten in cpu/gpu. I am still a bit nervous about needing to maintain a set of maths functions but if we can also contain them in one place in core (relying on the concrete implementations of arithmetic primitives elsewhere if needed), that would be good.

 * TODO:
 * OLD UNPACKERS
 * which I am hestitant to switch to the CPU-style until I better
 * understand why the explicit gpu_qcomp instantiation is necessary
 * (iirc static HIP structs have a different alignment than qcomp?!)

If it's helpful I can volunteer @eessmann as tribute a second pair of eyes and hands to dig into this? We now have a task in QATCH (phase 2 of the Quantum Software Lab) which is broadly "help maintain QuEST (and other scalable emulators)", that Erich is assigned to.

[TysonRayJones]

Moving away from the standard library entirely is philosphically upsetting, but clearly has practical benefits!

Yea it's outrageous but at least it won't keep us up at night about esoteric performance pitfalls! (We'll trade that for alignment nightmares)

I strongly support this proposal. Structurally it would be nice to have the common interface in core to then be specialised/overwritten in cpu/gpu. I am still a bit nervous about needing to maintain a set of maths functions but if we can also contain them in one place in core (relying on the concrete implementations of arithmetic primitives elsewhere if needed), that would be good.

Made common definitions in the previous commit. We can still discretionarily fallback to the std::complex operators within the CPU custom maths, paying the NaN performance penalties, when convenient (as pow in cpu_qcomp.hpp currently does).

If it's helpful I can volunteer @eessmann as tribute a second pair of eyes and hands to dig into this? We now have a task in QATCH (phase 2 of the Quantum Software Lab) which is broadly "help maintain QuEST (and other scalable emulators)", that Erich is assigned to.

That'd be super helpful! I preserved the method (used in relation to that comment) in the GPU backend, so I'm hoping things just work on HIP. I don't have an AMD machine handy to test on, so it would be terrific if Erich can give it a whirl! I can also re-setup the (paid) GPU Github Action runners if billing has been sorted out.

Note the new cpu_qcomp type might introduce a similar issue (misalignment of static arrays); currently, the qcomp[2] and qcomp[4] within a DiagMatr1 and DiagMatr2 are being converted/decayed to cpu_qcomp pointers. An explicit copy might be necessary, like in the GPU case, but I've so far been unable to trigger the alignment problem.

[TysonRayJones]

You can see the custom backend-specific maths/overloads here and here - so far, it's just pow!

[eessmann]

Heya @TysonRayJones,
Happy to help, I have access to HIP on Dirac so I give it a test on MI300s. Out of curiosity, why not go with thrust complex for our gpu complex data type?

[TysonRayJones]

Happy to help, I have access to HIP on Dirac so I give it a test on MI300s.

Wew that'd be brilliant!

Out of curiosity, why not go with thrust complex for our gpu complex data type?

It's a good question; previously, because using thrust complex inside CUDA kernels has some pitfalls and (iirc) didn't help solve the operator overload issues already seen of cuComplex (where HIP/rocThrust defines arithmetic operator overloads but CUDA/Thrust doesn't). It was my original intention to use thrust complex for QuEST 4.0, but was forced to switch to cuComplex (although I've forgotten precisely why).

With this PR however, there is a more compelling reason not to use thrust-complex; we must already define a custom complex type for the CPU backend to avoid compiler-specific performance pitfalls. Since we have to define its arithmetic overloads, we may as well re-use this (base) type for the GPU backend too, for better symmetry. We should now never have any issues related to operator overloads on particular platforms like HIP (but I do live in fear of alignment issues, as alluded to!)

[otbrown]

We can still discretionarily fallback to the std::complex operators within the CPU custom maths, paying the NaN performance penalties, when convenient

Since Erich's identified the relevant fortran-rules flag for all the main compilers, there should be no penalty! Can always just add the flag to the specific relevant source file(s) since it won't be needed elsewhere. Of course if we can't identify the compiler we can just omit.

I can also re-setup the (paid) GPU Github Action runners if billing has been sorted out.

Sadly, we don't have a budget we could use for that, but one of the three of us (Me, Erich, @JPRichings) can be the CI runner 😉

[JPRichings]

+ state:

V3:

QuEST/QuEST/src/GPU/QuEST_gpu.cu

Line 434 in d4f75f7

void statevec_initPlusState(Qureg qureg)

V4.2:

QuEST/quest/src/api/initialisations.cpp

Line 59 in 9d7618d

void initPlusState(Qureg qureg) {

QuEST/quest/src/core/localiser.cpp

Line 729 in 9d7618d

void localiser_statevec_initUniformState(Qureg qureg, qcomp amp) {

QuEST/quest/src/core/accelerator.cpp

Line 1131 in 9d7618d

void accel_statevec_initUniformState_sub(Qureg qureg, qcomp amp) {

QuEST/quest/src/gpu/gpu_subroutines.cpp

Line 1866 in 9d7618d

void gpu_statevec_initUniformState_sub(Qureg qureg, qcomp amp) {

QuEST/quest/src/gpu/gpu_thrust.cuh

Line 1051 in 9d7618d

void thrust_statevec_initUniformState(Qureg qureg, cu_qcomp amp) {

Given that the performance is comparable here and the kernel is very simple and is just setting values in the state vector probably points to how some of the changes to gpu kernels in v4. I assume there is a common pattern for gate application across the gates and something there is adding additional overhead that wasn't present in v3?

I need to check but I guess many of the gates end up here:

QuEST/quest/src/gpu/gpu_kernels.cuh

Line 517 in 9d7618d

__global__ void kernel_statevec_anyCtrlAnyTargDiagMatr_sub(

[JPRichings]

Really dumb guess based on eye balling the code is that:

QuEST/quest/src/core/bitwise.hpp

Line 206 in 9d7618d

INLINE qindex insertBitsWithMaskedValues(qindex number, int* bitInds, int numBits, qindex mask) {

is common in most of the statevector kernels and is called by each thread.

hmmm for loop in GPU kernel...

QuEST/quest/src/core/bitwise.hpp

Line 164 in 9d7618d

INLINE qindex insertBits(qindex number, int* bitIndices, int numIndices, int bitValue) {

[JPRichings]

I also notice that this masking logic is also present in the CPU side changes as noted in the examples here: #638

with insertBitsWithMaskedValues in a parallel region in the v4 CPU side code.

The performance discrepancy in #720 also seems present in the same range of qubits < 20...

Profiling required!

[JPRichings]

hmmm...

QuEST/quest/src/gpu/gpu_subroutines.cpp

Line 700 in 9d7618d

/// @todo

[JPRichings]

Ok so after a fair amount of digging I think we are moving data from host to device just before every kernel call that has ctrls.

Until this line,

QuEST/quest/src/gpu/gpu_subroutines.cpp

Line 708 in 9d7618d

devints deviceCtrls = util_getSorted(ctrls);

I think Ctrls are host side data but on the return of util_getSorted they are placed in a devints type which is a thrust device_vector this causes a host to device copy of data just before the call to the kernel function.

using devints = thrust::device_vector<int>;

int* getPtr(devints& qubits) {

    return thrust::raw_pointer_cast(qubits.data());
}

getPtr is called in the call to the kernel function:

QuEST/quest/src/gpu/gpu_subroutines.cpp

Line 713 in 9d7618d

getPtr(deviceCtrls), ctrls.size(), ctrlStateMask, getPtr(deviceTargs), targs.size(),

Comparing this to V3.7 is the ctrl mask was captured by value in the kernel launch as a long long int so transferred in kernel call to device (I think):

QuEST/QuEST/src/GPU/QuEST_gpu.cu

Line 875 in d4f75f7

void statevec_multiControlledTwoQubitUnitary(Qureg qureg, long long int ctrlMask, int q1, int q2, ComplexMatrix4 u)

and lived in constant memory which will probably have less of a performance impact than relying on thrust to move the data on creation of the vector.

[JPRichings]

report1.nsys-rep.txt
nsight systems profile of a qft circuit quest V4:

cuda malloc: ~ 230 microsecs
cuda memcpyasync: ~ 8 microsec
Kernel call: ~ 9 microsec

Run on grace-hopper node so mem copy from host to device will be really quick and much faster than in above results.

qft output:

Total number of gates: 210
Measured probability amplitude of |0..0> state: 9.53674e-07
Calculated probability amplitude of |0..0>, C0 = 1 / 2^20: 9.53674e-07
Measuring final state: (all probabilities should be 0.5)
Qubit 0 measured in state 1 with probability 0.5
Qubit 1 measured in state 1 with probability 0.5
Qubit 2 measured in state 1 with probability 0.5
Qubit 3 measured in state 1 with probability 0.5
Qubit 4 measured in state 1 with probability 0.5
Qubit 5 measured in state 1 with probability 0.5
Qubit 6 measured in state 1 with probability 0.5
Qubit 7 measured in state 0 with probability 0.5
Qubit 8 measured in state 1 with probability 0.5
Qubit 9 measured in state 1 with probability 0.5
Qubit 10 measured in state 1 with probability 0.5
Qubit 11 measured in state 1 with probability 0.5
Qubit 12 measured in state 1 with probability 0.5
Qubit 13 measured in state 1 with probability 0.5
Qubit 14 measured in state 1 with probability 0.5
Qubit 15 measured in state 1 with probability 0.5
Qubit 16 measured in state 1 with probability 0.5
Qubit 17 measured in state 0 with probability 0.5
Qubit 18 measured in state 1 with probability 0.5
Qubit 19 measured in state 1 with probability 0.5

Final state:
|11111110111111111011>
QFT run time: 0.0691252s
Total run time: 0.386635s
Generating '/tmp/nsys-report-50e2.qdstrm'
[1/1] [========================100%] report1.nsys-rep
Generated:

[TysonRayJones]

Apologies for delay - everything's on fire, I'll likely not have time to return until next weekend.
Your analysis looks right to me!

I think Ctrls are host side data but on the return of util_getSorted they are placed in a devints type which is a thrust device_vector this causes a host to device copy of data just before the call to the kernel function.

Absolutely right - a common Thrust paradigm but a gross pitfall in my opinion. I hate the implicit staging! Eager to change this if this is such an opportunity.

cuda malloc: ~ 230 microsecs
cuda memcpyasync: ~ 8 microsec
Kernel call: ~ 9 microsec

Ahh nicely done! The logic for handling controls differs between v3 and v4; the latter is asymptotically superior.
In v3, we enumerate every local amplitude and skip those failing the control conditions (permits the ctrls to be passed as a simple bitmask).
In v4, we enumerate strictly those amplitudes which satisfy the control condition. This leads to a factor 1/2^num_controls fewer iterations (a big asymptotic saving!) but so far requires the controls are kept in a sorted list (for bitwise insertion). You've pointed out re-allocating the GPU memory to stage that list is extremely non-negligible, eep!

It'd be amazing if we can perform that "enumerate only control-satisfying amps" logic using a control bitmask, possibly using these primitives?? Otherwise, we can easily maintain a persistent negligibly-sized device workspace, like we already do here.

[TysonRayJones]

Noting though that we stil see a performance regression when there are no control qubits - could it be that a superfluous zero-size allocation is happening in that scenario? 😅 Certainly there's no algorithmic overhead because of the NumCtrls templating (I thiiiiiink)

[JPRichings]

I think I have a fix where we move the Ctrls over to device into constant memory using cudaMemcpyToSymbol. Just profiling now. Looks fast could be wrong.

[JPRichings]

I caution this with I could have made a horrible error in how I am doing the data transfer and messed up the correctness but based on the attached profile the performance is much improved:

The call to memcopy to symbol only takes 200ns!

The kernel now takes: ~ 2-3 micro seconds (probably as accessing constant read only memory is much quicker)

There is now no allocation!

Profile: report_symbol.nsys-rep.txt

Great caution until I confirm correctness on monday!

Given there is another large cuda malloc call (~250 microsec) still in the profile I think there is more performance to be had in this.

Hardware target same grace-hopper node as before.

qft output:

Total number of gates: 210
Measured probability amplitude of |0..0> state: 9.53674e-07
Calculated probability amplitude of |0..0>, C0 = 1 / 2^20: 9.53674e-07
Measuring final state: (all probabilities should be 0.5)
Qubit 0 measured in state 0 with probability 0.5
Qubit 1 measured in state 0 with probability 0.5
Qubit 2 measured in state 0 with probability 0.5
Qubit 3 measured in state 1 with probability 0.5
Qubit 4 measured in state 1 with probability 0.5
Qubit 5 measured in state 0 with probability 0.5
Qubit 6 measured in state 0 with probability 0.5
Qubit 7 measured in state 1 with probability 0.5
Qubit 8 measured in state 0 with probability 0.5
Qubit 9 measured in state 1 with probability 0.5
Qubit 10 measured in state 1 with probability 0.5
Qubit 11 measured in state 0 with probability 0.5
Qubit 12 measured in state 1 with probability 0.5
Qubit 13 measured in state 0 with probability 0.5
Qubit 14 measured in state 0 with probability 0.5
Qubit 15 measured in state 1 with probability 0.5
Qubit 16 measured in state 0 with probability 0.5
Qubit 17 measured in state 0 with probability 0.5
Qubit 18 measured in state 0 with probability 0.5
Qubit 19 measured in state 0 with probability 0.5

Final state:
|00011001011010010000>
QFT run time: 0.00877581s
Total run time: 0.318082s
Generating '/tmp/nsys-report-e8a0.qdstrm'
[1/1] [========================100%] report1.nsys-rep
Generated:

[JPRichings]

Same issue in:

QuEST/quest/src/gpu/gpu_subroutines.cpp

Line 289 in 9d7618d

void gpu_statevec_anyCtrlOneTargDenseMatr_subA(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, int targ, CompMatr1 matr) {

22 instances of devints in gpu_subroutines.cpp my guess is that each one will be causing the same issue.

[JPRichings]

Further improvement if 1 additional kernel updated: kernel_statevec_anyCtrlOneTargDenseMatr_subA

Same qft executable same grace-hopper hardware.

Profile:
report_symbol2.nsys-rep.txt

Total number of gates: 210
Measured probability amplitude of |0..0> state: 9.53674e-07
Calculated probability amplitude of |0..0>, C0 = 1 / 2^20: 9.53674e-07
Measuring final state: (all probabilities should be 0.5)
Qubit 0 measured in state 1 with probability 0.5
Qubit 1 measured in state 1 with probability 0.5
Qubit 2 measured in state 1 with probability 0.5
Qubit 3 measured in state 0 with probability 0.5
Qubit 4 measured in state 0 with probability 0.5
Qubit 5 measured in state 1 with probability 0.5
Qubit 6 measured in state 0 with probability 0.5
Qubit 7 measured in state 1 with probability 0.5
Qubit 8 measured in state 0 with probability 0.5
Qubit 9 measured in state 0 with probability 0.5
Qubit 10 measured in state 0 with probability 0.5
Qubit 11 measured in state 0 with probability 0.5
Qubit 12 measured in state 0 with probability 0.5
Qubit 13 measured in state 0 with probability 0.5
Qubit 14 measured in state 1 with probability 0.5
Qubit 15 measured in state 0 with probability 0.5
Qubit 16 measured in state 1 with probability 0.5
Qubit 17 measured in state 0 with probability 0.5
Qubit 18 measured in state 1 with probability 0.5
Qubit 19 measured in state 0 with probability 0.5

Final state:
|11100101000000101010>
QFT run time: 0.00287859s
Total run time: 0.318814s

[JPRichings]

Note to check gpu_thrust.cuh as

QuEST/quest/src/gpu/gpu_thrust.cuh

Line 1015 in 9d7618d

devints devQubits = qubits;

QuEST/quest/src/gpu/gpu_thrust.cuh

Line 787 in 9d7618d

qreal thrust_statevec_calcProbOfMultiQubitOutcome_sub(Qureg qureg, vector<int> qubits, vector<int> outcomes) {

Causing applyQubitMeasurementAndGetProb to be impacted.

[TysonRayJones]

Wew super neat to have this addressed, nice one! 🎉

For this PR, I'll proceed under the assumption there's no CUDA performance impact of gpu_qcomp. Still need to...

• check this is true with HIP also (I'd bet fine!)
• check cpu_qcomp has no performance impact (i.e. do all my tests again for CPU, with and without multithreading).

The pain of these manual cross-platform performance regression checks begs for some automated benchmarker - but it's of course extremely unsexy work. Maybe a good unitaryHack bounty, or a chore for someone interested in experimenting with AI sloppers ¯\(ツ)/¯

[otbrown]

begs for some automated benchmarker

We've been pondering this a bit -- it might be something we can set up through EIDF. I've also been considering reviewing all tests to confirm they're MPI-safe and then running MPI tests through EIDF as well.

All very much post-v4.3 activity though I think!

[JPRichings]

This is the pull request that corrected the performance issues on nvidia not tested on AMD #739

[TysonRayJones]

Note to self: can probably verify cpu_qcomp performance using Meisam's scripts in #720

[TysonRayJones]

Serial MacOS ARM tests confirm this achieves the same performance as the -Ofast flags

[TysonRayJones]

After disabling the relevant v4 middle-end overheads (see issue #720), v4 multithreaded performance (8 threads on a Mac M4 Pro) matches v3. Note it outperforms v3 in the small end because of the v4 autodeployer which auto-disables multithreaded process of <8-qubit statevectors 🎉

All my concerns with this (c|g)pu_qcomp refactor have been addressed, except AMD/HIP correctness/performance. Should really check my tomfoolery with stack-memory casting is okay! I might be able to get my hands on an AMD eGPU quickly though 🙏

[TysonRayJones]

Confirmed it all runs fine on HIP! 🎉 Doing final checks of stack-alloc memory casting (which may have funky alignment pitfalls)