多线程Red-Black Gauss Seidel求解不可压缩性步骤运行缓慢问题

问题描述 #

在C++实现的PIC/FLIP流体模拟中，尝试通过线程池将基于Red-Black Gauss Seidel算法的不可压缩性求解步骤多线程化。采用水平拆分网格的方式分配任务给多线程，分红、黑两轮迭代计算，但多线程版本性能远低于单线程：10线程版本耗时46613毫秒，单线程仅226毫秒，且扩展性能极差。更换为原生线程实现后，问题依然存在。

核心代码 #

多线程求解核心函数 #

void solveIncompressibilityRedBlackForColumns(const float dt, const float overRelaxation, const int startColumn, const int endColumn, const bool red) {
    int n = this->numY;
    for (int i = startColumn; i < endColumn; ++i) {
        for (int j = 1; j < numY - 1; ++j) {
            if (red) { 
                if ((i + j) % 2 != 0) {
                    continue;
                }
            }
            if (!red) { 
                if ((i + j) % 2 == 0) {
                    continue;
                }
            } 
            if (this->cellType[i * n + j] == FLUID_CELL) {
                float leftType = cellType[(i - 1) * n + j] <= AIR_CELL ? 1: 0;
                float rightType = cellType[(i + 1) * n + j] <= AIR_CELL ? 1: 0;
                float topType = cellType[i * n + j - 1] <= AIR_CELL ? 1 : 0;
                float bottomType = cellType[i * n + j + 1] <= AIR_CELL ? 1: 0;
                float divideBy = leftType + rightType + topType + bottomType;
                if (divideBy != 0.f) {
                    float divergence;
                    divergence = this->u[(i + 1) * n + j] - this->u[i * n + j]+ this->v[i * n + j + 1] - this->v[i * n + j];
                    if (this->particleRestDensity > 0.f) {
                        float k = 10.f; // 3, 10
                        float compression = this->particleDensity[i * n + j] -this->particleRestDensity;
                        if (compression > 0.f) {
                            divergence = divergence - k * compression;
                        }
                    }
                    float p = divergence / divideBy;
                    p *= overRelaxation;
                    this->u[i * n + j] += leftType * p;
                    this->u[(i + 1) * n + j] -= rightType * p;
                    this->v[i * n + j] += topType * p;
                    this->v[i * n + j + 1] -= bottomType * p;
                }
            }
        }
    }
}

多线程调度代码 #

for (int _ = 0; _ < 10; ++_) { 

    for (int i = 0; i < numThreads; ++i) {
        if (i != numThreads - 1) {
            thread_pool.addTask([&, this, i]() {
                this->solveIncompressibilityRedBlackForRows(sdtoverRelaxation, i * numRowsEachThread + 1, i*numRowsEachThread + numRowsEachThread + 1, true);
            });
        }
        else {
            thread_pool.addTask([&, this, i]() {
                this->solveIncompressibilityRedBlackForRows(sdtoverRelaxation, i * numRowsEachThread + 1, i*numRowsEachThread + numRowsEachThread + numMissedRows + 1, true);
            });
        }
    }
    thread_pool.waitForCompletion();
    
    for (int i = 0; i < numThreads; ++i) {
        if (i != numThreads - 1) {
            thread_pool.addTask([&, this, i]() {
                this->solveIncompressibilityRedBlackForRows(sdtoverRelaxation, i * numRowsEachThread + 1, i*numRowsEachThread + numRowsEachThread + 1, false);
            });
        }
        else {
            thread_pool.addTask([&, this, i]() {
                this->solveIncompressibilityRedBlackForRows(sdt, overRelaxation, i * numRowsEachThread + 1, i * numRowsEachThread + numRowsEachThread + numMissedRows + 1, false);
            });
        }
    }
    thread_pool.waitForCompletion();
}

复现代码及运行结果 #

复现代码 #

#include <thread>
#include <vector>
#include <chrono>
#include <iostream>

void solveIncompressibilityRedBlackForRows(const float  overRelaxation, const int startColumn, const int endColumn, const bool red, const int gridHeight, std::vector<float>& u, std::vector<float>& v) {
    int n = gridHeight;
    for (int i = startColumn; i < endColumn; ++i) {
        for (int j = 1; j < n - 1; ++j) {
        
            if (red) { 
                if ((i + j) % 2 != 0) {
                    continue;
                }
            }
            if (!red) { 
                if ((i + j) % 2 == 0) {
                    continue;
                }
            } 
            float leftType = 1;
            float rightType = 1;
            float topType = 1;
            float bottomType = 1;
            float divideBy = 4;
            if (divideBy != 0.f) {
                float divergence;
                divergence = u[(i + 1) * n + j] - u[i * n + j] + v[i * n + j + 1] - v[i * n + j];
                float p = divergence / divideBy;
                p *= overRelaxation;
                u[i * n + j] += leftType * p;
                u[(i + 1) * n + j] -= rightType * p;
                v[i * n + j] += topType * p;
                v[i * n + j + 1] -= bottomType * p;
            }
        }
    }
}

int main() {
    int gridHeight = 100;
    int gridWidth = 100;
    std::vector<float> u(gridHeight * gridWidth, 10);
    std::vector<float> v(gridHeight * gridWidth, 10);

    int numThreads = 10;

    std::vector<std::thread> threads;

    float overRelaxation = 1.9;
    int numRowsEachThread = (gridWidth - 2) / numThreads;
    int numMissedRows = (gridWidth - 2) - numRowsEachThread * numThreads;

    auto start = std::chrono::high_resolution_clock::now();
    // main loop to run the incompressibility stuff
    for (int i = 0; i < 1000; ++i) {
        for (int i = 0; i < numThreads; ++i) {
            if (i != numThreads - 1) {
                threads.push_back(std::thread(solveIncompressibilityRedBlackForRows, overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + 1, true, gridHeight,  std::ref(u), std::ref(v)));
            }
            else {
                threads.push_back(std::thread(
                    solveIncompressibilityRedBlackForRows, overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + numMissedRows+ 1, true, gridHeight,  std::ref(u), std::ref(v)));
            }
        }
        
        for (auto& thread : threads) {
            thread.join();
        }

        threads.clear();
        
        for (int i = 0; i < numThreads; ++i) {
            if (i != numThreads - 1) {
                threads.push_back(std::thread(
                    solveIncompressibilityRedBlackForRows, overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + 1, false, gridHeight,  std::ref(u), std::ref(v)));
            }
            else {
                threads.push_back(std::thread(
                    solveIncompressibilityRedBlackForRows, overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + numMissedRows+ 1, false, gridHeight,  std::ref(u), std::ref(v)));
            }
        }
            
        for (auto& thread : threads) {
            thread.join();
        }

        threads.clear();
    }

    auto end = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    std::cout << "Time taken for 10 threads: " << duration.count() << " milliseconds" << std::endl;

    // single threaded program to compare
    start = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < 1000; ++i) {
        for (int i = 0; i < numThreads; ++i) {
            if (i != numThreads - 1) {
                solveIncompressibilityRedBlackForRows(overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + 1, true, gridHeight,  std::ref(u), std::ref(v));
            }
            else {
                solveIncompressibilityRedBlackForRows(overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + numMissedRows+ 1, true, gridHeight,  std::ref(u), std::ref(v));
            }
        }

        for (int i = 0; i < numThreads; ++i) {
            if (i != numThreads - 1) {
                solveIncompressibilityRedBlackForRows(overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + 1, false, gridHeight,  std::ref(u), std::ref(v));
            }
            else {
                solveIncompressibilityRedBlackForRows(overRelaxation, i * numRowsEachThread + 1, i *numRowsEachThread + numRowsEachThread + numMissedRows+ 1, false, gridHeight,  std::ref(u), std::ref(v));
            }
        }
    }
    end = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);

    std::cout << "Time taken for single thread: " << duration.count() << " milliseconds" << std::endl;
}

运行结果 #

Time taken for 10 threads: 46613 milliseconds
Time taken for single thread: 226 milliseconds

性能低下原因分析 #

1. 线程创建销毁开销爆炸：复现代码中每轮迭代都创建10个新线程，完成后销毁，1000次迭代共创建销毁20000个线程。线程的创建、调度、销毁涉及内核态操作，开销极大，完全覆盖了并行计算带来的收益。
2. 伪共享（False Sharing）：多个线程同时操作u、v数组中相邻的元素，这些元素大概率处于同一个CPU缓存行。当一个线程修改缓存行内的元素时，其他线程的缓存行会被标记为失效，需要重新从内存读取，频繁的缓存同步导致性能暴跌。
3. 任务粒度不合理：100x100的网格拆分给10个线程，每个线程仅处理约10行数据，任务粒度太小，线程调度开销占比远高于计算本身，并行效率极低。
4. 循环分支效率低下：循环内的红黑判断使用模运算和嵌套分支，影响CPU分支预测效率，增加单线程计算的耗时。

优化方案 #

1. 复用线程池，消除线程创建销毁开销 #

提前创建固定数量的线程（比如等于CPU核心数），每次迭代仅向线程池提交任务，避免重复创建销毁线程。示例自定义线程池：

class ThreadPool {
public:
    ThreadPool(size_t numThreads) {
        for (size_t i = 0; i < numThreads; ++i) {
            workers.emplace_back([this] {
                while (true) {
                    std::function<void()> task;
                    {
                        std::unique_lock<std::mutex> lock(this->queueMutex);
                        this->condition.wait(lock, [this] { return this->stop || !this->tasks.empty(); });
                        if (this->stop && this->tasks.empty()) return;
                        task = std::move(this->tasks.front());
                        this->tasks.pop();
                    }
                    task();
                }
            });
        }
    }

    template<class F>
    void enqueue(F&& f) {
        {
            std::unique_lock<std::mutex> lock(queueMutex);
            tasks.emplace(std::forward<F>(f));
        }
        condition.notify_one();
    }

    ~ThreadPool() {
        {
            std::unique_lock<std::mutex> lock(queueMutex);
            stop = true;
        }
        condition.notify_all();
        for (std::thread& worker : workers) worker.join();
    }

private:
    std::vector<std::thread> workers;
    std::queue<std::function<void()>> tasks;
    std::mutex queueMutex;
    std::condition_variable condition;
    bool stop = false;
};

2. 优化数据布局，解决伪共享 #

• 缓存行对齐：将u、v数组按缓存行大小（通常64字节）对齐，避免相邻元素被不同线程修改导致缓存失效：

alignas(64) std::vector<float> u(gridHeight * gridWidth, 10);
alignas(64) std::vector<float> v(gridHeight * gridWidth, 10);

• 块划分间隔：划分线程任务时，让每个线程处理的块之间间隔一个缓存行对应的元素数量，减少跨线程的缓存冲突。

3. 调整任务粒度 #

根据网格大小和CPU核心数调整线程数量，比如100x100的网格使用2-4线程即可，避免线程过多导致调度开销过高；当网格增大到500x500以上时，再增加线程数至CPU核心数。

4. 优化循环逻辑 #

• 替换模运算：用位运算(i + j) & 1替代(i + j) % 2，提升计算效率；
• 合并分支判断：将红黑判断合并为一个条件，减少嵌套分支：

if ((red && ((i+j)&1) != 0) || (!red && ((i+j)&1) == 0)) {
    continue;
}

• 提取不变量：将循环内的不变量（如n=gridHeight）提取到循环外，避免重复读取。

5. SIMD指令加速计算 #

使用SSE/AVX等SIMD指令批量处理浮点运算，提升单线程计算效率。示例AVX计算片段：

__m256 u_right = _mm256_loadu_ps(&u[(i+1)*n + j]);
__m256 u_left = _mm256_loadu_ps(&u[i*n + j]);
__m256 v_bottom = _mm256_loadu_ps(&v[i*n + j+1]);
__m256 v_top = _mm256_loadu_ps(&v[i*n + j]);
__m256 divergence = _mm256_add_ps(_mm256_sub_ps(u_right, u_left), _mm256_sub_ps(v_bottom, v_top));

内容的提问来源于stack exchange，提问作者Daniel K