红黑高斯塞德尔和OpenMP
Red-Black Gauss Seidel and OpenMP
我试图证明OpenMP与MPICH相比的一个观点,我制作了下面的示例来演示在OpenMP中实现高性能是多么容易。
高斯-赛德尔迭代被分成两个独立的运行,这样在每次扫描中,每个操作都可以以任何顺序执行,并且每个任务之间应该没有依赖关系。因此,从理论上讲,每个处理器都不应该等待另一个进程执行任何类型的同步。
我遇到的问题是,独立于问题大小,我发现只有2个处理器的微弱加速,超过2个处理器的加速甚至可能更慢。许多其他的线性并行程序,我可以获得非常好的缩放,但这一个是棘手的。
我担心的是,我无法向编译器"解释"我在数组上执行的操作是线程安全的,因此它无法真正有效。
请看下面的例子。
谁有任何线索,如何使这更有效与OpenMP?
void redBlackSmooth(std::vector<double> const & b,
std::vector<double> & x,
double h)
{
// Setup relevant constants.
double const invh2 = 1.0/(h*h);
double const h2 = (h*h);
int const N = static_cast<int>(x.size());
double sigma = 0;
// Setup some boundary conditions.
x[0] = 0.0;
x[N-1] = 0.0;
// Red sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 1; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2/2.0)*(b[i] - sigma);
}
// Black sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 2; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2/2.0)*(b[i] - sigma);
}
}
添加:我现在也尝试了一个原始指针实现,它具有与使用STL容器相同的行为,因此可以排除它是来自STL的一些伪临界行为。
首先,确保x
向量与缓存边界对齐。我做了一些测试,如果我强制内存对齐,我在我的机器(core duo)上得到了100%的改进:
double * x;
const size_t CACHE_LINE_SIZE = 256;
posix_memalign( reinterpret_cast<void**>(&x), CACHE_LINE_SIZE, sizeof(double) * N);
第二,您可以尝试为每个线程分配更多的计算(这样您可以保持缓存行分开),但我怀疑openmp已经在底层做了类似的事情,所以对于大n来说它可能毫无价值。
在我的例子中,当x
不是缓存对齐时,这个实现要快得多。
const int workGroupSize = CACHE_LINE_SIZE / sizeof(double);
assert(N % workGroupSize == 0); //Need to tweak the code a bit to let it work with any N
const int workgroups = N / workGroupSize;
int j, base , k, i;
#pragma omp parallel for shared(b, x) private(sigma, j, base, k, i)
for ( j = 0; j < workgroups; j++ ) {
base = j * workGroupSize;
for (int k = 0; k < workGroupSize; k+=2)
{
i = base + k + (redSweep ? 1 : 0);
if ( i == 0 || i+1 == N) continue;
sigma = -invh2* ( x[i-1] + x[i+1] );
x[i] = ( h2/2.0 ) * ( b[i] - sigma );
}
}
总之,你肯定有缓存争夺的问题,但考虑到openmp的工作方式(遗憾的是我不熟悉它),它应该足以使用正确分配的缓冲区。
我认为主要问题是关于你使用的数组结构类型。让我们尝试将结果与向量和数组进行比较。(数组=使用new操作符的c-数组).
向量和数组的大小为N = 10000000。我强制平滑函数重复,以保持运行时间> 0.1秒。
Vector Time: 0.121007 Repeat: 1 MLUPS: 82.6399
Array Time: 0.164009 Repeat: 2 MLUPS: 121.945
MLUPS = ((N-2)*repeat/runtime)/1000000 (Million Lattice Points Update per second)
当涉及到网格计算时,MFLOPS会产生误导。在基本方程中做一些改变就可以考虑在相同的运行时实现高性能。
修改后的代码:
double my_redBlackSmooth(double *b, double* x, double h, int N)
{
// Setup relevant constants.
double const invh2 = 1.0/(h*h);
double const h2 = (h*h);
double sigma = 0;
// Setup some boundary conditions.
x[0] = 0.0;
x[N-1] = 0.0;
double runtime(0.0), wcs, wce;
int repeat = 1;
timing(&wcs);
for(; runtime < 0.1; repeat*=2)
{
for(int r = 0; r < repeat; ++r)
{
// Red sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 1; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2*0.5)*(b[i] - sigma);
}
// Black sweep.
#pragma omp parallel for shared(b, x) private(sigma)
for (int i = 2; i < N-1; i+=2)
{
sigma = -invh2*(x[i-1] + x[i+1]);
x[i] = (h2*0.5)*(b[i] - sigma);
}
// cout << "In Array: " << r << endl;
}
if(x[0] != 0) dummy(x[0]);
timing(&wce);
runtime = (wce-wcs);
}
// cout << "Before division: " << repeat << endl;
repeat /= 2;
cout << "Array Time:t" << runtime << "t" << "Repeat:t" << repeat
<< "tMLUPS:t" << ((N-2)*repeat/runtime)/1000000.0 << endl;
return runtime;
}
我没有改变代码中的任何东西,除了数组类型。为了获得更好的缓存访问和阻塞,您应该查看数据对齐(_mm_malloc)。
相关文章:
- 在C++上实现高斯赛德尔迭代方法
- C++:矩阵高斯消除不起作用:使用单维数组来存储元素
- (C++)(Visual Studio) 将高斯模糊滤镜应用于 RGB 中的灰度图像
- 用于创建高斯随机数的 c++ 函数
- 将高斯模糊应用于灰度图像
- 高斯雅各比法
- 用高斯Seidel红色黑色求解1D泊松方程
- 如何改变高斯分布(提升)中的种子
- 并行化高斯模糊链
- 为什么在 C++ 中实现高斯勒让德算法没有产生结果
- 生成高斯噪声
- 部分透视/高斯消除 - 交换列而不是产生错误输出的行
- 使用高斯模糊进行图像卷积,可以加速
- 如何使用 1D 高斯内核在 Filter2D 上创建自定义 2D 内核
- ***检测到堆栈粉碎***:在高斯消除中
- 使用C 加入高斯噪声
- 高斯滤波器核值
- 高斯-塞德尔法计算3个系统的线性方程
- 红黑高斯塞德尔和OpenMP
- 求解非线性方程组的高斯-塞德尔方法