基于OpenMP与std::experimental::simd实现曼德博集合的技术咨询

Great question! Using C++'s experimental SIMD extensions is a fantastic middle ground between messy manual intrinsics and underperforming basic #pragma omp simd for complex loops like Mandelbrot. Let's break this down step by step to get your implementation working clearly and efficiently.

First: What is std::experimental::simd<double>? #

Let's clear up the confusion first:

• native_simd<double> is not a single double — it's a vector type that packs multiple double values into a single SIMD register, optimized for your target CPU.
• The number of doubles it holds depends on your CPU's SIMD capabilities:
  • AVX2: 4 doubles per vector
  • AVX-512: 8 doubles per vector
  • SSE: 2 doubles per vector
• The compiler automatically picks the right size at compile time, which makes it perfect for your "compile-time unknown boundary" scenario.

Modified computeMandelbrot Implementation #

Here's how to adapt your existing OpenMP code to use std::experimental::simd, keeping readability a priority:

#include <experimental/simd>
namespace stdx = std::experimental; // Shorthand for readability

void Mandelbrot::computeMandelbrot() {
    using SimdDouble = stdx::native_simd<double>;
    const std::size_t simd_width = SimdDouble::size(); // Number of doubles per vector

    double dx = (d_xrange.second - d_xrange.first) / d_res.first;
    double dy = (d_yrange.second - d_yrange.first) / d_res.second;

    #pragma omp parallel for schedule(dynamic)
    for (std::size_t row = 0; row != d_res.second; ++row) {
        const double c_imag = d_yrange.first + row * dy;
        const std::size_t total_cols = d_res.first;
        std::size_t col = 0;

        // Process full SIMD chunks first
        for (; col + simd_width <= total_cols; col += simd_width) {
            // Create a vector of c_real values for this chunk of columns
            SimdDouble c_real;
            for (std::size_t i = 0; i < simd_width; ++i) {
                c_real[i] = d_xrange.first + (col + i) * dx;
            }

            SimdDouble real = 0.0;
            SimdDouble imag = 0.0;
            SimdDouble realSquared = 0.0;
            SimdDouble imagSquared = 0.0;
            stdx::native_simd<std::size_t> iter = 0;

            // Mandelbrot iteration loop (vectorized)
            while (true) {
                // Check if all elements have exited or hit max iter
                const auto mask = (iter < d_maxIter) && (realSquared + imagSquared < 4.0);
                if (!stdx::any_of(mask)) break;

                realSquared = real * real;
                imagSquared = imag * imag;
                imag = 2 * real * imag + c_imag; // Broadcast c_imag to all vector elements
                real = realSquared - imagSquared + c_real;
                iter += mask; // Only increment iter for elements still in the loop
            }

            // Write vector results back to the array
            for (std::size_t i = 0; i < simd_width; ++i) {
                d_iters[row * total_cols + col + i] = static_cast<double>(iter[i]);
            }
        }

        // Handle remaining columns (non-SIMD chunk)
        for (; col < total_cols; ++col) {
            double real = 0.0;
            double imag = 0.0;
            double realSquared = 0.0;
            double imagSquared = 0.0;
            const double c_real = d_xrange.first + col * dx;
            std::size_t iter = 0;

            while (iter < d_maxIter && realSquared + imagSquared < 4.0) {
                realSquared = real * real;
                imagSquared = imag * imag;
                imag = 2 * real * imag + c_imag;
                real = realSquared - imagSquared + c_real;
                ++iter;
            }
            d_iters[row * total_cols + col] = iter;
        }
    }
}

Key Implementation Notes #

Let's walk through the important parts to make sure you understand:

1. Vector Type Setup: We define SimdDouble as stdx::native_simd<double> and get its width with SimdDouble::size() — this tells us how many columns to process per SIMD chunk.
2. Chunked Column Processing: We split the column loop into two parts:
   • Full SIMD chunks: Process simd_width columns at once using vector operations.
   • Remaining columns: Fall back to the original scalar code for any leftover columns that don't fill a full vector.
3. Vectorized Iteration:
   • We use a mask (mask = (iter < d_maxIter) && (realSquared + imagSquared < 4.0)) to track which elements are still in the Mandelbrot loop.
   • stdx::any_of(mask) checks if any elements are still active — if not, we exit the loop.
   • iter += mask increments only the elements that are still iterating (the mask acts as a vector of booleans converted to 0/1).
4. Broadcasting Scalars: When we do imag = 2 * real * imag + c_imag, the scalar c_imag is automatically broadcast to all elements of the vector — no extra code needed!

Compilation Tips #

To compile this code, you'll need:

• A compiler that supports C++ experimental SIMD extensions (GCC 9+, Clang 12+).
• Flags to enable SIMD:
  • For GCC: -fexperimental-simd -mavx2 (or -mavx512f for AVX-512)
  • For Clang: -fsimd -mavx2
• Make sure your resolution's width is a multiple of simd_width (as you mentioned) to minimize the leftover scalar chunk — this will maximize performance.

Why This Works for HPC Beginners #

This approach balances:

• Readability: No messy intrinsics (like _mm256_add_pd), just familiar C++ syntax with vector types.
• Performance: Automatically leverages your CPU's full SIMD capabilities.
• Future-Proofing: Uses a C++ standard experimental feature, which is more forward-looking than manual intrinsics.

You can easily compare this implementation against your original OpenMP version, or other HPC paradigms (CUDA, MPI) to show tradeoffs in complexity vs. performance for your tutorial.

内容的提问来源于stack exchange，提问作者Nigel Overmars