Авто-векторизация в C с GCC 14

// Скалярный код
for (int i = 0; i < 8; i++) {
    array[i] += 5;
}
 
// Эквивалентные векторные операции (упрощённо)
__m256 vec = _mm256_load_ps(array);       // Загрузить 8 float'ов
__m256 constant = _mm256_set1_ps(5.0f);   // Подготовить вектор из восьми пятёрок
__m256 result = _mm256_add_ps(vec, constant); // Сложить векторы
_mm256_store_ps(array, result);           // Сохранить результат

// Тяжело векторизуемый код из-за зависимостей
for (int i = 2; i < n; i++) {
    fib[i] = fib[i-1] + fib[i-2];
}

// Выделение выровненной памяти
float *aligned_array = (float *)aligned_alloc(32, size * sizeof(float));
 
// Указание компилятору, что данные выровнены
__attribute__((aligned(32))) float another_array[1024];

void vector_add(float * restrict a, float * restrict b, 
               float * restrict result, int size) {
    for (int i = 0; i < size; i++) {
        result[i] = a[i] + b[i];
    }
}

// Умножение матриц
for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        float sum = 0;
        for (int k = 0; k < n; k++) {
            sum += a[i][k] * b[k][j];
        }
        c[i][j] = sum;
    }
}

void update(float *a, int n, int stride) {
    for (int i = 0; i < n-stride; i++) {
        a[i] = a[i+stride] * 2.0f;
    }
}

if (stride >= 8) {
    // Векторизованная версия
    ...
} else {
    // Обычная скалярная версия
    ...
}

// Векторизованная часть для блоков по 8 элементов
for (i = 0; i < n/8*8; i += 8) {
    // Векторные операции
}
 
// Остаток обрабатывается скалярно
for (; i < n; i++) {
    // Скалярные операции
}

// Неэффективный доступ к памяти (strided access)
for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        sum += matrix[j][i]; // Неоптимально: перепрыгиваем через строки
    }
}
 
// После loop interchange:
for (int j = 0; j < n; j++) {
    for (int i = 0; i < n; i++) {
        sum += matrix[j][i]; // Оптимально: последовательный доступ
    }
}

// Скалярный код с условием
for (int i = 0; i < n; i++) {
    if (a[i] > threshold) {
        b[i] = a[i] * 2.0f;
    } else {
        b[i] = a[i] * 0.5f;
    }
}
 
// Векторизованный эквивалент (псевдокод)
for (int i = 0; i < n; i += 8) {
    vec_a = load(a + i);
    mask = compare_gt(vec_a, vec_threshold);
    
    vec_result1 = multiply(vec_a, vec_2_0); // Для true-ветки
    vec_result2 = multiply(vec_a, vec_0_5); // Для false-ветки
    
    vec_b = blend(mask, vec_result1, vec_result2);
    store(b + i, vec_b);
}

// Скалярная редукция (сумма элементов)
float sum = 0.0f;
for (int i = 0; i < n; i++) {
    sum += a[i];
}

// Векторизованная редукция (псевдокод)
vec_sum1 = vec_0;
vec_sum2 = vec_0;
vec_sum3 = vec_0;
vec_sum4 = vec_0;
 
for (i = 0; i < n; i += 16) {
    vec_sum1 = vec_add(vec_sum1, load(a + i));
    vec_sum2 = vec_add(vec_sum2, load(a + i + 4));
    vec_sum3 = vec_add(vec_sum3, load(a + i + 8));
    vec_sum4 = vec_add(vec_sum4, load(a + i + 12));
}
 
// Горизонтальное сложение векторов
vec_sum = vec_add(vec_add(vec_sum1, vec_sum2), 
                 vec_add(vec_sum3, vec_sum4));
                 
// Финальное горизонтальное сложение элементов внутри вектора
sum = horizontal_add(vec_sum);

void scale_array(int *arr, int n) {
    for (int i = 0; i < n; i++) {
        arr[i] = arr[i] * 10; // Может переполниться для некоторых значений
    }
}

void scale_array_safe(int *arr, int n) {
    for (int i = 0; i < n; i++) {
        arr[i] = (int)((long long)arr[i] * 10); // Безопасное вычисление
    }
}

// Принудительная векторизация цикла
#pragma GCC ivdep
for (int i = 0; i < n; i++) {
    // Компилятор будет игнорировать потенциальные зависимости
    result[i] = a[i] + b[i];
}
 
// Запрет векторизации
#pragma GCC novector
for (int i = 0; i < n; i++) {
    // Этот цикл не будет векторизован
    complex_function(a[i]);
}

// Использование векторных типов GCC
typedef float float4 __attribute__((vector_size(16))); // 4 floats = 16 bytes
 
void vector_add(float4 *a, float4 *b, float4 *c, int n) {
    for (int i = 0; i < n; i++) {
        c[i] = a[i] + b[i]; // Компилируется в SIMD-инструкции
    }
}

// Обычная реализация на C
void vector_add_c(float *a, float *b, float *result, int size) {
    for (int i = 0; i < size; i++) {
        result[i] = a[i] + b[i];
    }
}
 
// Реализация с ручными SIMD-инструкциями (AVX)
void vector_add_avx(float *a, float *b, float *result, int size) {
    int i = 0;
    
    // Обработка по 8 элементов за раз
    for (; i <= size - 8; i += 8) {
        __m256 va = _mm256_loadu_ps(a + i);
        __m256 vb = _mm256_loadu_ps(b + i);
        __m256 vresult = _mm256_add_ps(va, vb);
        _mm256_storeu_ps(result + i, vresult);
    }
    
    // Оставшиеся элементы
    for (; i < size; i++) {
        result[i] = a[i] + b[i];
    }
}

void rgb_to_grayscale(unsigned char *rgb, unsigned char *gray, int pixel_count) {
    for (int i = 0; i < pixel_count; i++) {
        int r = rgb[i*3];
        int g = rgb[i*3 + 1];
        int b = rgb[i*3 + 2];
        
        // Стандартная формула яркости
        gray[i] = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b);
    }
}

void convolve(float *input, float *output, float *kernel, 
             int width, int height, int kernel_size) {
    int radius = kernel_size / 2;
    
    for (int y = radius; y < height - radius; y++) {
        for (int x = radius; x < width - radius; x++) {
            float sum = 0.0f;
            
            for (int ky = 0; ky < kernel_size; ky++) {
                for (int kx = 0; kx < kernel_size; kx++) {
                    int input_idx = (y + ky - radius) * width + (x + kx - radius);
                    sum += input[input_idx] * kernel[ky * kernel_size + kx];
                }
            }
            
            output[y * width + x] = sum;
        }
    }
}

float array_sum(float *arr, int size) {
    float sum = 0.0f;
    for (int i = 0; i < size; i++) {
        sum += arr[i];
    }
    return sum;
}

// Базовая реализация фильтра скользящего среднего
void moving_average_base(float *input, float *output, int window_size, int data_size) {
for (int i = 0; i < data_size; i++) {
    float sum = 0.0f;
    int count = 0;
    
    for (int j = max(0, i - window_size + 1); j <= i; j++) {
        sum += input[j];
        count++;
    }
    
    output[i] = sum / count;
}
}

// Поиск и обработка выбросов
void process_outliers(float *data, float *result, int size, float threshold) {
for (int i = 0; i < size; i++) {
    if (fabs(data[i]) > threshold) {
        result[i] = 0.0f; // Заменяем выбросы нулями
    } else {
        result[i] = data[i]; // Оставляем нормальные значения без изменений
    }
}
}

// Каждый нуклеотид кодируется 2 битами (A=00, C=01, G=10, T=11)
void find_pattern(uint8_t *sequence, uint8_t *pattern, 
              bool *matches, int seq_length, int pat_length) {
for (int i = 0; i <= seq_length - pat_length; i++) {
    bool match = true;
    for (int j = 0; j < pat_length; j++) {
        // Извлекаем 2 бита для текущего нуклеотида
        uint8_t nucl = (sequence[i+j/4] >> (6 - 2*(j%4))) & 0x03;
        if (nucl != pattern[j]) {
            match = false;
            break;
        }
    }
    matches[i] = match;
}
}

// Этот код может быть векторизован через SLP без циклов
float a = v1 * v2;
float b = v3 * v4;
float c = v5 * v6;
float d = v7 * v8;
float result = a + b + c + d;

void complex_transform(float *data, int size) {
for (int i = 1; i < size - 1; i++) {
    data[i] = (data[i-1] + data[i] + data[i+1]) / 3.0f; // Усреднение
}
}

void complex_transform_vectorizable(float *data, int size) {
float *temp = malloc(size * sizeof(float));
memcpy(temp, data, size * sizeof(float));
 
for (int i = 1; i < size - 1; i++) {
    data[i] = (temp[i-1] + temp[i] + temp[i+1]) / 3.0f;
}
 
free(temp);
}

// Принудительное разворачивание цикла до указанного фактора
#pragma GCC unroll 8
for (...) { ... }
 
// Выравнивание данных для эффективной векторизации
#pragma GCC vector_aligned
for (...) { ... }
 
// Указание ожидаемой длины цикла (помогает планировщику)
#pragma GCC loop_optimize(size, "256")
for (...) { ... }

// Подсказка компилятору, что функция хороший кандидат для встраивания
__attribute__((always_inline)) void hot_function() { ... }
 
// Подсказка, что функция обычно вызывается с определёнными значениями
int calculate(int x) __attribute__((target_clones("avx2","avx512f","default")));

// Неоптимальное чередование векторных и скалярных операций
for (int i = 0; i < size; i++) {
    // Векторизуемая часть
    float temp = a[i] * b[i];
    
    // Скалярная часть (сложные вычисления)
    if (temp > threshold) {
        result[i] = complex_function(temp);
    } else {
        result[i] = temp;
    }
}
 
// Оптимизированная версия с разделением операций
// Сначала выполняем все векторные операции
for (int i = 0; i < size; i++) {
    temp[i] = a[i] * b[i];
}
 
// Затем все скалярные
for (int i = 0; i < size; i++) {
    if (temp[i] > threshold) {
        result[i] = complex_function(temp[i]);
    } else {
        result[i] = temp[i];
    }
}

// Указание компилятору развернуть цикл
#pragma GCC unroll 8
for (int i = 0; i < size; i++) {
    result[i] = a[i] * b[i];
}

// Оригинальный код, сложный для анализа
float calculate(float x) {
    if (x < 0) return 0;
    if (x > 1) return 1;
    return x * x * (3 - 2 * x); // Сглаживание Эрмита
}
 
void process_array(float *input, float *output, int size) {
    for (int i = 0; i < size; i++) {
        output[i] = calculate(input[i]);
    }
}
 
// Версия с "векторной подсказкой"
void process_array_vectorized(float *input, float *output, int size) {
    // Создаём временные массивы для промежуточных вычислений
    float *temp1 = malloc(size * sizeof(float));
    float *temp2 = malloc(size * sizeof(float));
    
    // Разбиваем сложную функцию на простые векторизуемые шаги
    for (int i = 0; i < size; i++) {
        temp1[i] = input[i] < 0 ? 0 : (input[i] > 1 ? 1 : input[i]);
    }
    
    for (int i = 0; i < size; i++) {
        temp2[i] = temp1[i] * temp1[i];
    }
    
    for (int i = 0; i < size; i++) {
        output[i] = temp2[i] * (3 - 2 * temp1[i]);
    }
    
    free(temp1);
    free(temp2);
}

// Array of Structures (AoS) - менее удобно для векторизации
struct Particle {
    float x, y, z;
    float vx, vy, vz;
};
Particle particles[1000];
 
// Structure of Arrays (SoA) - гораздо удобнее для векторизации
struct ParticleSystem {
    float x[1000], y[1000], z[1000];
    float vx[1000], vy[1000], vz[1000];
};

// Стандартное суммирование столбцов матрицы
void sum_columns(float **matrix, float *result, int rows, int cols) {
    for (int j = 0; j < cols; j++) {
        float sum = 0;
        for (int i = 0; i < rows; i++) {
            sum += matrix[i][j]; // Плохая локальность, прыжки по памяти
        }
        result[j] = sum;
    }
}
 
// Транспонированный алгоритм
void sum_columns_transposed(float **matrix, float *result, int rows, int cols) {
    // Инициализируем результаты нулями
    for (int j = 0; j < cols; j++) {
        result[j] = 0;
    }
    
    // Аккумулируем вклад каждой строки в соответствующие столбцы
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            result[j] += matrix[i][j]; // Последовательный доступ по строке
        }
    }
}

// Стандартная реализация поиска подстроки (не векторизуется)
int find_substring(char *text, int text_len, char *pattern, int pattern_len) {
    for (int i = 0; i <= text_len - pattern_len; i++) {
        int j;
        for (j = 0; j < pattern_len; j++) {
            if (text[i + j] != pattern[j]) {
                break;
            }
        }
        if (j == pattern_len) return i; // Нашли совпадение
    }
    return -1; // Не нашли
}
 
// Версия с опережающим вычислением (хорошо векторизуется)
int find_substring_vectorized(char *text, int text_len, char *pattern, int pattern_len) {
    // Предварительная фильтрация: проверяем только первый символ
    int candidates[1024]; // Буфер для потенциальных совпадений
    int num_candidates = 0;
    
    // Этот цикл хорошо векторизуется
    for (int i = 0; i <= text_len - pattern_len; i++) {
        if (text[i] == pattern[0]) {
            candidates[num_candidates++] = i;
            if (num_candidates >= 1024) break;
        }
    }
    
    // Проверка кандидатов (обычно их намного меньше, чем общая длина текста)
    for (int c = 0; c < num_candidates; c++) {
        int pos = candidates[c];
        int j;
        for (j = 1; j < pattern_len; j++) { // Начинаем с j=1, т.к. первый символ уже проверен
            if (text[pos + j] != pattern[j]) {
                break;
            }
        }
        if (j == pattern_len) return pos;
    }
    
    return -1;
}

# В CMakeLists.txt
add_custom_target(vectorization_report
    COMMAND ${CMAKE_COMMAND} -E make_directory ${CMAKE_BINARY_DIR}/vec_reports
    COMMAND ${CMAKE_CXX_COMPILER} -O3 -march=native -fopt-info-vec-all ${CMAKE_CURRENT_SOURCE_DIR}/hot_path.c -o ${CMAKE_BINARY_DIR}/vec_reports/hot_path.s -S
    COMMENT "Generating vectorization report for hot_path.c"
)

// Стандартная версия с разрежеными доступами
void sparse_operation(float *values, int *indices, float *result, int n) {
    for (int i = 0; i < n; i++) {
        result[indices[i]] += values[i]; // Разреженный доступ
    }
}
 
// Версия с векторной предагрегацией
void sparse_operation_vectorized(float *values, int *indices, float *result, int size, int result_size) {
    // Разбиваем на блоки по 1024 элемента
    const int block_size = 1024;
    
    for (int start = 0; start < size; start += block_size) {
        int end = start + block_size < size ? start + block_size : size;
        
        // Локальный буфер для агрегации результатов текущего блока
        float local_result[MAX_RESULT_SIZE] = {0};
        
        // Этот цикл теперь работает с последовательным доступом к памяти
        for (int i = start; i < end; i++) {
            local_result[indices[i]] += values[i];
        }
        
        // Объединяем локальные результаты с глобальными
        // Этот цикл хорошо векторизуется
        for (int j = 0; j < result_size; j++) {
            result[j] += local_result[j];
        }
    }
}

// Гарантируем выравнивание на 32-байтную границу (для AVX)
float *aligned_array = (float *)aligned_alloc(32, size * sizeof(float));
 
// Для больших массивов можно выровнять и на границу страницы памяти
void *page_aligned = aligned_alloc(4096, large_size);

// Сложная функция, плохо поддающаяся векторизации
float complex_function(float x) {
    return sinf(x) * expf(-x*x) / (1.0f + x*x);
}
 
// Векторизуемая декомпозиция
void complex_function_vec(float *x, float *result, int n) {
    float temp1[MAX_SIZE], temp2[MAX_SIZE], temp3[MAX_SIZE];
    
    // Каждый из этих циклов хорошо векторизуется
    for (int i = 0; i < n; i++) temp1[i] = sinf(x[i]);
    for (int i = 0; i < n; i++) temp2[i] = expf(-x[i]*x[i]);
    for (int i = 0; i < n; i++) temp3[i] = 1.0f + x[i]*x[i];
    
    for (int i = 0; i < n; i++) {
        result[i] = temp1[i] * temp2[i] / temp3[i];
    }
}

// Объявляем несколько версий функции для разных наборов инструкций
__attribute__((target("default")))
void process_data_default(float *data, int size) {
    // Базовая реализация без SIMD
}
 
__attribute__((target("sse4.1")))
void process_data_sse41(float *data, int size) {
    // Версия с SSE4.1
}
 
__attribute__((target("avx2")))
void process_data_avx2(float *data, int size) {
    // Версия с AVX2
}
 
__attribute__((target("avx512f")))
void process_data_avx512(float *data, int size) {
    // Версия с AVX-512
}
 
// Функция-диспетчер, которая выбирает подходящую реализацию
void process_data(float *data, int size) {
    // GCC автоматически генерирует код выбора оптимальной версии
    process_data_default(data, size);
}

// Традиционная реализация с неявными зависимостями
float complex_calculation(float *input, int size) {
float sum1 = 0, sum2 = 0;
for (int i = 0; i < size; i++) {
    float temp = sqrt(input[i]);
    sum1 += temp;
    sum2 += temp * temp;
}
return sum1 * sum2;
}
 
// DAG-ориентированная реализация для лучшей векторизации
float complex_calculation_dag(float *input, int size) {
// Промежуточные результаты (вершины графа)
float *sqrt_values = malloc(size * sizeof(float));
float *squared_values = malloc(size * sizeof(float));
float sum1 = 0, sum2 = 0;
 
// Первый этап: вычисление промежуточных значений
for (int i = 0; i < size; i++) {
    sqrt_values[i] = sqrt(input[i]);
}
 
// Второй этап: возведение в квадрат (независимо от первого этапа)
for (int i = 0; i < size; i++) {
    squared_values[i] = sqrt_values[i] * sqrt_values[i];
}
 
// Третий этап: суммирование (хорошо векторизуется)
for (int i = 0; i < size; i++) {
    sum1 += sqrt_values[i];
}
 
// Четвёртый этап: ещё одно суммирование (тоже хорошо векторизуется)
for (int i = 0; i < size; i++) {
    sum2 += squared_values[i];
}
 
free(sqrt_values);
free(squared_values);
return sum1 * sum2;
}

// matrix_benchmark.c - комплексный пример автовекторизации
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
 
// Используем атрибут выравнивания для лучшей векторизации
typedef float matrix_t __attribute__((aligned(32)));
 
// Функция умножения матриц, хорошо поддающаяся векторизации
void matrix_multiply(matrix_t *a, matrix_t *b, matrix_t *c, int n) {
    memset(c, 0, n * n * sizeof(matrix_t));
    
    // Используем оптимальный порядок циклов для локальности данных
    for (int i = 0; i < n; i++) {
        for (int k = 0; k < n; k++) {
            // Извлекаем элемент A[i][k] в регистр для повторного использования
            matrix_t a_ik = a[i*n + k];
            
            // Внутренний цикл отлично векторизуется
            #pragma GCC ivdep
            for (int j = 0; j < n; j++) {
                c[i*n + j] += a_ik * b[k*n + j];
            }
        }
    }
}
 
// Функция фильтрации с условными операциями
void threshold_filter(matrix_t *matrix, int size, float threshold) {
    #pragma GCC unroll 8
    for (int i = 0; i < size; i++) {
        // Тернарный оператор обычно лучше векторизуется, чем if-else
        matrix[i] = (matrix[i] < threshold) ? 0.0f : matrix[i];
    }
}
 
// Конвертация в градации серого с использованием весовых коэффициентов
void rgb_to_grayscale(unsigned char *rgb, unsigned char *gray, int pixel_count) {
    // Используем restrictпрямо указать на отсутствие пересечений
    const unsigned char * restrict rgb_ptr = rgb;
    unsigned char * restrict gray_ptr = gray;
    
    for (int i = 0; i < pixel_count; i++) {
        int r = rgb_ptr[i*3];
        int g = rgb_ptr[i*3 + 1];
        int b = rgb_ptr[i*3 + 2];
        
        // Стандартная формула яркости
        gray_ptr[i] = (unsigned char)(0.299f * r + 0.587f * g + 0.114f * b);
    }
}
 
int main() {
    const int N = 1024;
    const int IMG_SIZE = 1920 * 1080;
    
    // Выделяем выровненную память для лучшей производительности векторных операций
    matrix_t *a = (matrix_t*)aligned_alloc(32, N * N * sizeof(matrix_t));
    matrix_t *b = (matrix_t*)aligned_alloc(32, N * N * sizeof(matrix_t));
    matrix_t *c = (matrix_t*)aligned_alloc(32, N * N * sizeof(matrix_t));
    
    unsigned char *rgb = (unsigned char*)malloc(IMG_SIZE * 3);
    unsigned char *gray = (unsigned char*)malloc(IMG_SIZE);
    
    // Инициализация данных...
    
    // Замер времени выполнения
    clock_t start = clock();
    
    // Основные вычисления
    matrix_multiply(a, b, c, N);
    threshold_filter(c, N*N, 0.5f);
    rgb_to_grayscale(rgb, gray, IMG_SIZE);
    
    clock_t end = clock();
    double elapsed = (double)(end - start) / CLOCKS_PER_SEC;
    
    printf("Время выполнения: %.3f сек\n", elapsed);
    
    // Освобождение памяти
    free(a);
    free(b);
    free(c);
    free(rgb);
    free(gray);
    
    return 0;
}