SIMD Vectors - Complex numbers?

Question

I am wanting to use the System.Numerics.Vector namespace in .NET Core, but I have run into the issue of there being no support for vectors of complex numbers. As it stands, the Vector type will work with any primitive type, from byte to double. I am not a programmer by trade, so I can be ignorant of some low-level/conceptual stuff, but is there a reason why there isn't support for complex numbers? As far as I can see, perhaps the only issue would be that the Complex type is a managed struct. Could I not just extend the Register type to include [FieldOffset(0)] internal Complex complex_0 and build new functions around that?

I am willing to work on the type extension myself, but wanted to ask if there is some reason why this wasn't included in the first place, since it seems like a lot of signal processing work would benefit greatly from complex vectors and SIMD.

Answer 1

The System.Numerics.Vector is useless for your goal. You could rely on the Vector256 class from the System.Runtime.Intrinsics namespace. The necessary instructions are exposed in the System.Runtime.Intrinsics.X86.Avx class. (Vector128 would be of little help here, as it represents a single Complex number. I failed to locate any magic that would help me re-interpreting the Complex as Vector128 without copying, otherwise we would be able to boost up the built-in operators provided for the Complex class)

So, the idea is to

1. Have your Complex data expressed in terms of Span
2. Have the arithmetic operators implemented on those Spans
3. Internally, cast those Spans to Vector256
4. And use the Avx operations

See the sample code below:

public static class ComplexAvx
    {
        public static Span<Complex> Add(ReadOnlySpan<Complex> left, ReadOnlySpan<Complex> right)
        {
            var result = new Complex[Math.Min(left.Length, right.Length)].AsSpan();
            var vectorRes = MemoryMarshal.Cast<Complex, Vector256<double>>(result);
            var vectorLeft = MemoryMarshal.Cast<Complex, Vector256<double>>(left);
            var vectorRight = MemoryMarshal.Cast<Complex, Vector256<double>>(right);
            for (int i = 0; i < vectorRes.Length; i++)
                vectorRes[i] = Avx.Add(vectorLeft[i], vectorRight[i]);

            for (int i = 2 * vectorRes.Length; i < result.Length; i++)
                result[i] = left[i] + right[i];
            return result;
        }

        public static Span<Complex> Subtract(ReadOnlySpan<Complex> left, ReadOnlySpan<Complex> right)
        {
            var result = new Complex[Math.Min(left.Length, right.Length)].AsSpan();
            var vectorRes = MemoryMarshal.Cast<Complex, Vector256<double>>(result);
            var vectorLeft = MemoryMarshal.Cast<Complex, Vector256<double>>(left);
            var vectorRight = MemoryMarshal.Cast<Complex, Vector256<double>>(right);
            for (int i = 0; i < vectorRes.Length; i++)
                vectorRes[i] = Avx.Subtract(vectorLeft[i], vectorRight[i]);

            for (int i = 2 * vectorRes.Length; i < result.Length; i++)
                result[i] = left[i] - right[i];
            return result;
        }

        public static Span<Complex> Multiply(ReadOnlySpan<Complex> left, ReadOnlySpan<Complex> right)
        {
            var result = new Complex[Math.Min(left.Length, right.Length)].AsSpan();
            var vectorRes = MemoryMarshal.Cast<Complex, Vector256<double>>(result);
            var vectorLeft = MemoryMarshal.Cast<Complex, Vector256<double>>(left);
            var vectorRight = MemoryMarshal.Cast<Complex, Vector256<double>>(right);
            for (int i = 0; i < vectorRes.Length; i++)
            {
                var l = vectorLeft[i];  
                var r = vectorRight[i]; 
                vectorRes[i] = Avx.HorizontalAdd(
                    Avx.Multiply(
                        Avx.Multiply(l, r), 
                        Vector256.Create(1.0, -1.0, 1.0, -1.0)), 
                    Avx.Multiply(
                        l, 
                        Avx.Permute(r, 0b0101)
                        )); 
            }
            for (int i = 2 * vectorRes.Length; i < result.Length; i++)
                result[i] = left[i] * right[i];
            return result;
        }

    }

This code performs 5 AVX operations for multiplying two complex numbers, whereas using the System.Numerics.Complex.Multiply would consume 8 multiplications and 4 additions. (Of course, the code should check for Avx.IsEnabled - skipped that for brevity) I haven't bothered with the benchmarks, so not sure how much (if at all) gain would it produce on the modern CPUs. Avx512 would allow using same 5 instructions for handling 4 complex at once; but that's pending the CLR support yet.