tối ưu hóa cột tối đa khôn ngoan với SIMD

Question

Tôi có chức năng này trong đó tôi đã dành một lượng thời gian đáng kể cho mã của mình và tôi muốn tối ưu hóa nó bằng vectơ hóa-SIMD-compiler, nếu có thể.

Về cơ bản, nó tìm giá trị và vị trí của giá trị tối đa trên một ma trận trên các cột và lưu trữ chúng:

• val_ptr: ma trận đầu vào: column-major (kiểu Fortran) n_rows-by-n_cols (trong đó thường là n_rows >> n_cols)
• opt_pos_ptr: int vector có độ dài n_rows nơi lưu trữ vị trí của giá trị lớn nhất. Trên mục nhập có đầy các số không.
• max_ptr: vector float có độ dài n_rows nơi lưu trữ giá trị lớn nhất. Trên mục nhập có các bản sao của cột đầu tiên của val_ptr
• Hàm sẽ được gọi trong một vòng lặp song song
• Vùng bộ nhớ được đảm bảo không chồng chéo
• Tôi không thực sự cần max_ptr để được lấp đầy, hiện tại nó chỉ được sử dụng để ghi sổ và tránh cấp phát bộ nhớ
• Tôi sử dụng MSVC, C ++ 17 trên Windows 10. Có nghĩa là để chạy các CPU Intel hiện đại

Mã, trong đó loại mẫu có nghĩa là float hoặc double:

template <typename eT>
find_max(const int n_cols, 
         const int n_rows, 
         const eT* val_ptr,
         int* opt_pos_ptr,
         eT* max_ptr){
    for (int col = 1; col < n_cols; ++col)
    {
        //Getting the pointer to the beginning of the column
        const auto* value_col = val_ptr + col * n_rows;
        //Looping over the rows
        for (int row = 0; row < n_rows; ++row)
        {
            //If the value is larger than the current maximum, we replace and we store its positions
            if (value_col[row] > max_ptr[row])
            {
                max_ptr[row] = value_col[row];
                opt_pos_ptr[row] = col;
            }
        }
    }
}

Những gì tôi đã thử cho đến nay:

• Tôi đã cố gắng sử dụng song song OpenMP cho vòng lặp bên trong, nhưng chỉ mang lại thứ gì đó trên các hàng rất lớn, lớn hơn một chút so với mức sử dụng hiện tại của tôi.
• Vòng lặp if in inside ngăn #pragma omp simd hoạt động và tôi không thể viết lại nếu không có nó.

Answer 1

Dựa trên mẫu mã bạn đã đăng, có vẻ như bạn muốn tính giá trị tối đa theo chiều dọc, nghĩa là trong trường hợp của bạn, "cột" là theo chiều ngang. Trong C / C ++, dãy phần tử theo chiều ngang (tức là hai phần tử liền kề có khoảng cách bằng một phần tử trong bộ nhớ) thường được gọi là hàng và dọc (trong đó hai phần tử liền kề có khoảng cách bằng kích thước hàng trong bộ nhớ) - cột. Trong câu trả lời của tôi dưới đây, tôi sẽ sử dụng thuật ngữ truyền thống, trong đó các hàng là hàng ngang và các cột là hàng dọc.

Ngoài ra, để ngắn gọn, tôi sẽ tập trung vào một loại phần tử ma trận có thể có - float. Ý tưởng cơ bản là giống nhau double, với sự khác biệt chính là số phần tử trên mỗi vectơ và lựa chọn _ps/ _pdnội dung. Tôi sẽ cung cấp một phiên bản cho doubleở cuối.

---

Ý tưởng là bạn có thể tính toán tối đa theo chiều dọc cho nhiều cột song song bằng cách sử dụng _mm_max_ps/ _mm_max_pd. Để cũng ghi lại vị trí của giá trị lớn nhất được tìm thấy, bạn có thể so sánh giá trị tối đa trước đó với các phần tử hiện tại. Kết quả của phép so sánh là một mặt nạ, trong đó các phần tử là tất cả các phần tử ở đó giá trị tối đa được cập nhật. Mặt nạ đó cũng có thể được sử dụng để chọn vị trí cần được cập nhật.

Tôi phải lưu ý rằng thuật toán dưới đây giả định rằng việc ghi lại vị trí của phần tử max nào là không quan trọng, nếu có nhiều phần tử max bằng nhau trong một cột. Ngoài ra, tôi giả sử ma trận không chứa các giá trị NaN, điều này sẽ ảnh hưởng đến các phép so sánh. Thêm về điều này sau.

void find_max(const int n_cols, 
         const int n_rows, 
         const float* val_ptr,
         int* opt_pos_ptr,
         float* max_ptr){
    const __m128i mm_one = _mm_set1_epi32(1);

    // Pre-compute the number of rows that can be processed in full vector width.
    // In a 128-bit vector there are 4 floats or 2 doubles
    int tail_size = n_rows & 3;
    int n_rows_aligned = n_rows - tail_size;
    int row = 0;
    for (; row < n_rows_aligned; row += 4)
    {
        const auto* col_ptr = val_ptr + row;
        __m128 mm_max = _mm_loadu_ps(col_ptr);
        __m128i mm_max_pos = _mm_setzero_si128();
        __m128i mm_pos = mm_one;
        col_ptr += n_rows;
        for (int col = 1; col < n_cols; ++col)
        {
            __m128 mm_value = _mm_loadu_ps(col_ptr);

            // See if this value is greater than the old maximum
            __m128 mm_mask = _mm_cmplt_ps(mm_max, mm_value);
            // If it is, save its position
            mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, _mm_castps_si128(mm_mask));

            // Compute the maximum
            mm_max = _mm_max_ps(mm_value, mm_max);

            mm_pos = _mm_add_epi32(mm_pos, mm_one);
            col_ptr += n_rows;
        }

        // Store the results
        _mm_storeu_ps(max_ptr + row, mm_max);
        _mm_storeu_si128(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);
    }

    // Process tail serially
    for (; row < n_rows; ++row)
    {
        const auto* col_ptr = val_ptr + row;
        auto max = *col_ptr;
        int max_pos = 0;
        col_ptr += n_rows;
        for (int col = 1; col < n_cols; ++col)
        {
            auto value = *col_ptr;
            if (value > max)
            {
                max = value;
                max_pos = col;
            }

            col_ptr += n_rows;
        }

        max_ptr[row] = max;
        opt_pos_ptr[row] = max_pos;
    }
}

The code above requires SSE4.1 because of the blending intrinsics. You can replace those with a combination of _mm_and_si128/_ps, _mm_andnot_si128/_ps and _mm_or_si128/_ps, in which case the requirements will be lowered to SSE2. See Intel Intrinsics Guide for more details on the particular intrinsics, including which instruction set extensions they require.

---

Lưu ý về giá trị NaN. Nếu ma trận của bạn có thể có NaN, _mm_cmplt_pskiểm tra sẽ luôn trả về false. Về phần _mm_max_ps, nó thường không biết những gì nó sẽ trở lại. Lệnh maxpsmà nội tại dịch để trả về toán hạng thứ hai (nguồn) của nó nếu một trong hai toán hạng là NaN, vì vậy bằng cách sắp xếp các toán hạng của lệnh, bạn có thể đạt được một trong hai hành vi. Tuy nhiên, không có tài liệu nào ghi lại đối số nào của _mm_max_psnội tại đại diện cho toán hạng nào của lệnh và thậm chí có thể trình biên dịch có thể sử dụng kết hợp khác nhau trong các trường hợp khác nhau. Xem câu trả lời này để biết thêm chi tiết.

Để đảm bảo hành vi chính xác wrt. NaNs bạn có thể sử dụng trình hợp dịch nội tuyến để buộc thứ tự chính xác của các maxpstoán hạng. Thật không may, đó không phải là một tùy chọn với MSVC cho mục tiêu x86-64, mà bạn đã nói rằng bạn đang sử dụng, vì vậy thay vào đó bạn có thể sử dụng lại _mm_cmplt_pskết quả cho hỗn hợp thứ hai như sau:

// Compute the maximum
mm_max = _mm_blendv_ps(mm_max, mm_value, mm_mask);

Điều này sẽ ngăn chặn các NaN trong các giá trị tối đa thu được. Nếu bạn muốn giữ lại NaN, bạn có thể sử dụng phép so sánh thứ hai để phát hiện NaN:

// Detect NaNs
__m128 mm_nan_mask = _mm_cmpunord_ps(mm_value, mm_value);

// Compute the maximum
mm_max = _mm_blendv_ps(mm_max, mm_value, _mm_or_ps(mm_mask, mm_nan_mask));

---

Bạn có thể có thể cải thiện hơn nữa hiệu suất của thuật toán ở trên nếu bạn sử dụng vectơ rộng hơn ( __m256hoặc __m512) và bỏ cuộn vòng ngoài bằng một yếu tố nhỏ, để ít nhất một dòng trong bộ nhớ cache của dữ liệu hàng được tải trên mỗi lần lặp lại của vòng lặp bên trong.

---

Đây là một ví dụ về triển khai cho double. Điểm quan trọng cần lưu ý ở đây là bởi vì chỉ có hai doublephần tử trên mỗi vectơ và vẫn có bốn vị trí trên mỗi vectơ, chúng tôi phải mở vòng lặp bên ngoài để xử lý hai vectơ doublecùng một lúc và sau đó nén hai mặt nạ khỏi so sánh với mức tối đa trước đó để trộn các vị trí 32 bit.

void find_max(const int n_cols, 
         const int n_rows, 
         const double* val_ptr,
         int* opt_pos_ptr,
         double* max_ptr){
    const __m128i mm_one = _mm_set1_epi32(1);

    // Pre-compute the number of rows that can be processed in full vector width.
    // In a 128-bit vector there are 2 doubles, but we want to process
    // two vectors at a time.
    int tail_size = n_rows & 3;
    int n_rows_aligned = n_rows - tail_size;
    int row = 0;
    for (; row < n_rows_aligned; row += 4)
    {
        const auto* col_ptr = val_ptr + row;
        __m128d mm_max1 = _mm_loadu_pd(col_ptr);
        __m128d mm_max2 = _mm_loadu_pd(col_ptr + 2);
        __m128i mm_max_pos = _mm_setzero_si128();
        __m128i mm_pos = mm_one;
        col_ptr += n_rows;
        for (int col = 1; col < n_cols; ++col)
        {
            __m128d mm_value1 = _mm_loadu_pd(col_ptr);
            __m128d mm_value2 = _mm_loadu_pd(col_ptr + 2);

            // See if this value is greater than the old maximum
            __m128d mm_mask1 = _mm_cmplt_pd(mm_max1, mm_value1);
            __m128d mm_mask2 = _mm_cmplt_pd(mm_max2, mm_value2);
            // Compress the 2 masks into one
            __m128i mm_mask = _mm_packs_epi32(
                _mm_castpd_si128(mm_mask1), _mm_castpd_si128(mm_mask2));
            // If it is, save its position
            mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, mm_mask);

            // Compute the maximum
            mm_max1 = _mm_max_pd(mm_value1, mm_max1);
            mm_max2 = _mm_max_pd(mm_value2, mm_max2);

            mm_pos = _mm_add_epi32(mm_pos, mm_one);
            col_ptr += n_rows;
        }

        // Store the results
        _mm_storeu_pd(max_ptr + row, mm_max1);
        _mm_storeu_pd(max_ptr + row + 2, mm_max2);
        _mm_storeu_si128(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);
    }

    // Process 2 doubles at once
    if (tail_size >= 2)
    {
        const auto* col_ptr = val_ptr + row;
        __m128d mm_max1 = _mm_loadu_pd(col_ptr);
        __m128i mm_max_pos = _mm_setzero_si128();
        __m128i mm_pos = mm_one;
        col_ptr += n_rows;
        for (int col = 1; col < n_cols; ++col)
        {
            __m128d mm_value1 = _mm_loadu_pd(col_ptr);

            // See if this value is greater than the old maximum
            __m128d mm_mask1 = _mm_cmplt_pd(mm_max1, mm_value1);
            // Compress the mask. The upper half doesn't matter.
            __m128i mm_mask = _mm_packs_epi32(
                _mm_castpd_si128(mm_mask1), _mm_castpd_si128(mm_mask1));
            // If it is, save its position
            mm_max_pos = _mm_blendv_epi8(mm_max_pos, mm_pos, mm_mask);

            // Compute the maximum
            mm_max1 = _mm_max_pd(mm_value1, mm_max1);

            mm_pos = _mm_add_epi32(mm_pos, mm_one);
            col_ptr += n_rows;
        }

        // Store the results
        _mm_storeu_pd(max_ptr + row, mm_max1);
        // Only store the lower two positions
        _mm_storel_epi64(reinterpret_cast< __m128i* >(opt_pos_ptr + row), mm_max_pos);

        row += 2;
    }

    // Process tail serially
    for (; row < n_rows; ++row)
    {
        const auto* col_ptr = val_ptr + row;
        auto max = *col_ptr;
        int max_pos = 0;
        col_ptr += n_rows;
        for (int col = 1; col < n_cols; ++col)
        {
            auto value = *col_ptr;
            if (value > max)
            {
                max = value;
                max_pos = col;
            }

            col_ptr += n_rows;
        }

        max_ptr[row] = max;
        opt_pos_ptr[row] = max_pos;
    }
}