sd_fft: AVX512 implementation for some functions

[user202729]

This PR defines new types vec8dz and vec16dz, corresponding to
AVX512 register types, and use that to speed up some functions in sd_fft.c.

Only a limited number of functions in sd_fft.c are modified.
This is mostly to discuss what's the correct interface that we want in machine_vectors.h.
Note that vec8d is already taken. Changing vec8d from two ymm to one zmm register is not unconditionally profitable, since there would not be enough instruction level parallelism.

To limit the amount of code changes, I keep the original transformed element ordering, which in particular mean sd_fft_store_vec8dz_basecase_order becomes very terrible. (it takes only 24 cycles though.)

Original:

  bits coeffs   min
     5     32   203
     6     64   374
     7    128   697
     8    256  1477
     9    512  3155
    10   1024  6849

After first commit:

  bits coeffs   min
     5     32   203
     6     64   289
     7    128   695
     8    256  1218
     9    512  3145
    10   1024  5980

bits=6 shows 22% improvement and bits=8 shows 17% improvement.
bits=7 is unmodified (will need more work).

After second commit:

  bits coeffs   min
     5     32   203
     6     64   289
     7    128   702
     8    256  1110
     9    512  3294
    10   1024  5288

bits=8 shows 8% additional improvement.

---

The temp benchmark file (AI).

Details

#include <asm/unistd.h>
#include <errno.h>
#include <inttypes.h>
#include <linux/perf_event.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <sys/ioctl.h>
#include <sys/mman.h>
#include <sys/syscall.h>
#include <unistd.h>

#include "ulong_extras.h"
#include "fft_small.h"

typedef struct
{
    int fd;
    struct perf_event_mmap_page * pc;
    uint32_t ecx;
    uint16_t pmc_width;
} perf_cycle_counter_t;

static void
perf_cycle_counter_init(perf_cycle_counter_t * counter)
{
    struct perf_event_attr attr;

    memset(&attr, 0, sizeof(attr));
    attr.type = PERF_TYPE_HARDWARE;
    attr.size = sizeof(attr);
    attr.config = PERF_COUNT_HW_CPU_CYCLES;
    attr.exclude_kernel = 1;
    attr.exclude_hv = 1;

    counter->fd = syscall(__NR_perf_event_open, &attr, 0, -1, -1, 0);
    if (counter->fd < 0)
    {
        flint_printf("perf_event_open failed: %s\n", strerror(errno));
        flint_abort();
    }

    counter->pc = (struct perf_event_mmap_page *)
            mmap(NULL, 4096, PROT_READ, MAP_SHARED, counter->fd, 0);
    if (counter->pc == MAP_FAILED)
    {
        flint_printf("mmap failed: %s\n", strerror(errno));
        close(counter->fd);
        flint_abort();
    }

    if (!counter->pc->cap_user_rdpmc || counter->pc->index == 0)
    {
        flint_printf("userspace rdpmc is unavailable\n");
        munmap(counter->pc, 4096);
        close(counter->fd);
        flint_abort();
    }

    counter->ecx = counter->pc->index - 1;
    counter->pmc_width = counter->pc->pmc_width;

    ioctl(counter->fd, PERF_EVENT_IOC_RESET, 0);
    ioctl(counter->fd, PERF_EVENT_IOC_ENABLE, 0);
}

static void
perf_cycle_counter_clear(perf_cycle_counter_t * counter)
{
    ioctl(counter->fd, PERF_EVENT_IOC_DISABLE, 0);
    munmap(counter->pc, 4096);
    close(counter->fd);
}

static uint64_t
perf_cycle_counter_get_raw(const perf_cycle_counter_t * counter)
{
    uint32_t lo, hi;

    __asm__ volatile("lfence\n\trdpmc\n\tlfence"
                     : "=a"(lo), "=d"(hi)
                     : "c"(counter->ecx)
                     : "memory");

    return ((uint64_t) hi << 32) | lo;
}

static uint64_t
perf_cycle_counter_get_delta(const perf_cycle_counter_t * counter,
        uint64_t before, uint64_t after)
{
    return (after - before) & ((UINT64_C(1) << counter->pmc_width) - 1);
}

static uint64_t
time_sd_fft(perf_cycle_counter_t * counter, sd_fft_ctx_t Q, double * data,
        ulong L, ulong n)
{
    uint64_t before, after;

    before = perf_cycle_counter_get_raw(counter);
    sd_fft_trunc(Q, data, L, n, n);
    after = perf_cycle_counter_get_raw(counter);

    if (data[0] == -1.0)
        flint_printf("\r");

    return perf_cycle_counter_get_delta(counter, before, after);
}

int
main(void)
{
    const ulong reps = 1000;
    sd_fft_ctx_t Q;
    flint_rand_t state;
    perf_cycle_counter_t counter;

    flint_rand_init(state);
    sd_fft_ctx_init_prime(Q, UWORD(0x0003f00000000001));
    perf_cycle_counter_init(&counter);

    flint_printf("bits     coeffs       cycles\n");

    for (ulong i = 5; i <= 10; i++)
    {
        ulong n = UWORD(1) << i;
        uint64_t best = UINT64_MAX;
        double * data;

        sd_fft_ctx_fit_depth(Q, i);
        data = (double *) flint_aligned_alloc(64,
                FLINT_MAX(64, sd_fft_ctx_data_size(i) * sizeof(double)));

        for (ulong j = 0; j < n; j++)
            data[j] = n_randint(state, Q->mod.n);
        sd_fft_trunc(Q, data, i, n, n);

        for (ulong rep = 0; rep < reps; rep++)
        {
            uint64_t cycles;

            for (ulong j = 0; j < n; j++)
                data[j] = n_randint(state, Q->mod.n);

            cycles = time_sd_fft(&counter, Q, data, i, n);
            if (cycles < best)
                best = cycles;
        }

        flint_printf("%4wu %10wu %12" PRIu64 "\n", i, n, best);

        flint_aligned_free(data);
    }

    perf_cycle_counter_clear(&counter);
    sd_fft_ctx_clear(Q);
    flint_rand_clear(state);
    flint_cleanup_master();

    return 0;
}

Use perf events to count number of CPU cycles.
More reliable than rdtscp since clock frequency varies.
Might need sudo sysctl kernel.perf_event_paranoid=-1 to run this.

[user202729]

https://www.texmacs.org/joris/ntt/ntt.html claims 51.4ns (xeon) and 16.5ns (zen4) for r=64... I guess more room for improvement is possible.

[albinahlback]

https://www.texmacs.org/joris/ntt/ntt.html claims 51.4ns (xeon) and 16.5ns (zen4) for r=64... I guess more room for improvement is possible.

I believe Joris optimized through SLPs, i.e. no branching at all, in the framework of Mathemagix (i.e. their own compiler).

[user202729]

The sd_fft_basecase_6_0 and sd_fft_basecase_6_1 are straight-line programs too. There must be some other magic being done here...

[vneiger]

I tested performance through the profile files in fft_small/profile, on two different machines supporting avx512 (zen4 and Cascade Lake).

I observe a performance penalty on the p-sd_fft profile, on both machines, concerning fft (columns 1 vs 3, and 5 vs 7), and for ifft sometimes performance is better, sometimes worse (columns 2 vs 4 and 6 vs 8).

             zen4 BEFORE   zen4 AFTER    casc. BEFORE  casc. AFTER
depth 20 setup: 2 ms       3ms           2ms           2ms
      10.3:  39.10  39.10  42.65  37.04  16.96  17.82  25.59  16.76
      10.6:  42.45  42.45  48.11  42.45  18.51  18.75  27.23  17.60
      10.8:  42.11  42.11  49.13  40.94  17.76  18.20  25.41  16.94
depth 11.0:  45.46  45.46  51.73  44.12  18.52  19.74  26.32  18.75
      11.2:  42.31  39.89  48.15  39.89  17.67  17.03  24.07  16.05
      11.3:  42.89  41.63  48.80  40.44  18.14  17.92  25.27  16.85
      11.5:  43.41  40.93  49.39  40.93  18.13  17.91  25.13  16.46
      11.6:  45.26  42.60  49.94  42.60  18.57  18.81  26.33  17.66
      11.7:  44.33  35.68  48.76  43.02  18.06  18.75  24.38  17.84
      11.8:  39.89  42.17  46.13  42.17  18.68  19.42  25.45  18.45
      11.9:  42.53  45.11  48.02  43.78  18.38  19.33  25.23  18.38
depth 19.0:  46.79  48.25  51.47  46.79  17.55  19.06  20.05  18.83
      19.1:  39.89  38.81  46.33  41.03  15.78  14.96  18.41  14.96
      19.2:  43.43  37.87  47.64  42.19  15.88  15.38  18.46  15.54
      19.3:  44.10  42.80  48.51  41.58  16.17  15.99  18.90  15.99
      19.4:  42.53  42.53  38.17  34.62  16.18  15.67  18.84  15.67
      19.4:  44.71  43.43  42.23  37.08  16.52  16.35  19.24  16.35
      19.5:  43.39  43.39  44.70  43.39  16.58  16.96  19.16  16.96
      19.6:  44.00  42.74  48.26  33.24  16.62  16.62  19.18  16.62
      19.7:  44.50  43.23  44.50  42.03  16.63  16.81  19.15  17.00
      19.8:  44.86  38.13  47.67  43.58  16.95  17.14  19.31  17.53
      19.9:  35.64  35.64  47.89  45.07  16.84  17.41  19.65  17.61
      20.0:  41.44  43.81  49.46  45.10  17.04  18.04  19.41  18.25
depth 20.0:  44.94  44.94  47.66  47.66  16.91  18.50  19.66  18.72

However when launching other benchs, I do not see much difference and mostly no negative impact, sometimes not changed and sometimes a small improvement. For example, p-mul:

 --- fft_small 1 thread  --- 
       CASCADE LAKE BEFORE  CASCADE LAKE AFTER
20.00: 4.219  4.195  4.314  4.080  4.066  4.098
20.10: 4.123  4.113  4.144  3.965  3.959  3.969
20.20: 4.220  4.214  4.236  4.114  4.098  4.146
20.30: 4.068  4.046  4.156  3.943  3.918  4.046
20.40: 5.145  5.136  5.160  5.251  5.243  5.273
20.50: 5.084  5.070  5.108  5.142  5.126  5.162
20.60: 4.270  4.261  4.280  4.170  4.152  4.183
20.70: 4.235  4.219  4.255  4.243  4.220  4.267
20.80: 4.159  4.133  4.180  4.162  4.138  4.211
20.90: 4.282  4.269  4.290  4.238  4.217  4.256
       avg 0.177, max 4.12  avg 0.207, max 5.35

So, I suppose this PR needs further investigation/comments about the performance gains?

[user202729]

Seem pretty problematic, since changes such as this one are obviously architecture-dependent, and I don't have these machines to test on... (any idea?)

The benchmark in the original PR was only for these functions and only count clock cycles, so it's possible that a local improvement somehow (?) translate to a global regression. Seem pretty unlikely though.

Some hypotheses I can think of:

1. the code takes less clock, but clock speed becomes slower. (I think on the machine I test on I ruled this out by changing the setting to performance mode, the difference in clock speed is negligible, though in powersave mode clock cycle can vary significantly)

2. there's some interference, for an example, see "Mixing VEX and SSE code" in https://www.agner.org/optimize/optimizing_assembly.pdf . Seem unlikely though.

3. something else, such as the warmup time. https://www.agner.org/optimize/optimizing_assembly.pdf contains

   Warm up time
   Many processors are able to turn off the upper bits of the 256 bit or 512 bit execution units in order to save power when these units are not used.

   It takes typically 10-20 µs to power up this upper part after an idle period. The throughput of 256-bit vector instructions is much lower during this warm-up period because the processor uses the lower 128-bit units twice to execute a 256-bit operation.

4. This part is much slower on other machines.

   To limit the amount of code changes, I keep the original transformed element ordering, which in particular mean sd_fft_store_vec8dz_basecase_order becomes very terrible. (it takes only 24 cycles though.)