小编Z b*_*son的帖子

两个数组的点积

for(int i=0; i<n; i++) {
    sum += x[i]*y[i];
}

不重用数据,因此它应该是一个内存绑定操作.因此,我应该能够从点积测量内存带宽.

使用代码 为什么 - 矢量化 - 循环 - 没有性能改进 我的系统带宽为9.3 GB/s.但是,当我尝试使用点积计算带宽时,我获得单个线程的速率的两倍以及使用多个线程的速率超过三倍(我的系统有四个核心/八个超线程).这对我没有意义,因为内存绑定操作不应该受益于多个线程.以下代码的输出如下:

Xeon E5-1620, GCC 4.9.0, Linux kernel 3.13
dot 1 thread:      1.0 GB, sum 191054.81, time 4.98 s, 21.56 GB/s, 5.39 GFLOPS
dot_avx 1 thread   1.0 GB, sum 191043.33, time 5.16 s, 20.79 GB/s, 5.20 GFLOPS
dot_avx 2 threads: 1.0 GB, sum 191045.34, time 3.44 s, 31.24 GB/s, 7.81 GFLOPS
dot_avx 8 threads: 1.0 GB, sum 191043.34, time 3.26 …

我想解决Ax = b,其中A是一个非常大的正方形正定对称块矩阵和x和b是矢量.当我说大的时候我指的是nxn一个n大到的矩阵300,000.

这是一个我想要解决的更小但有代表性的矩阵的例子.

这里是相同的矩阵放大显示它由密集的matricies块组成.

我以前(看到这里,这里,并且在这里)使用了Eigen的Cholesky解算器,它工作得很好,n<10000但与n=300000Cholesky求解器相比太慢了.但是,我没有利用我有块矩阵的事实.显然,存在用于求解稀疏块矩阵的算法(例如,块cholesky分解).

我想特别知道Eigen是否已经使用因子分解或迭代方法对我可以使用的稀疏密集块矩阵进行了优化算法？

你也可以建议其他算法可能是解决我的矩阵的理想选择吗？我的意思是,据我所知,至少对于因子分解发现排列矩阵是NP难以存在很多不同的启发式方法,并且据我所知,人们建立了不同矩阵结构的直觉(例如带状矩阵)以及哪种算法效果最好跟他们.我还没有这种直觉.

我目前正在研究使用共轭梯度法.我自己用C++实现了这个,但仍然没有利用矩阵是一个块矩阵的事实.

//solve A*rk = xk
//Eigen::SparseMatrix<double> A;
//Eigen::VectorXd rk(n);
Eigen::VectorXd pk = rk;

double rsold = rk.dot(rk);
int maxiter = rk.size();
for (int k = 0; k < maxiter; k++) {
    Eigen::VectorXd Ap = A*pk;  
    double ak = rsold /pk.dot(Ap);
    xk += …

我正在测试台式机和服务器上的内存带宽。

Sklyake desktop 4 cores/8 hardware threads
Skylake server Xeon 8168 dual socket 48 cores (24 per socket) / 96 hardware threads

系统的峰值带宽为

Peak bandwidth desktop = 2-channels*8*2400 = 38.4 GB/s
Peak bandwidth server  = 6-channels*2-sockets*8*2666 = 255.94 GB/s

我正在使用自己的STREAM三合一函数来测量带宽（稍后将提供完整代码）

void triad(double *a, double *b, double *c, double scalar, size_t n) {
  #pragma omp parallel for
  for(int i=0; i<n; i++) a[i] = b[i] + scalar*c[i];
}

这是我得到的结果

         Bandwidth (GB/s)
threads  Desktop  Server         
1             28      16
2(24)         29     146 …

可以在CUDA文件中使用openMP pragma(不在内核代码中)吗？

我将结合gpu和cpu计算.但是nvvc编译器失败,"无法找到未知选项'openmp'",如果我将porgram与openmp选项链接(在linux下)

一种解决方法是仅在c/c ++文件中使用openMP-statments.

我用AVX一次计算八个点产品.在我目前的代码中,我做了类似的事情(在展开之前):

常春藤桥/桑迪桥

__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {        
    __m256 breg0 = _mm256_load_ps(&b[8*i]);
    tmp0 = _mm256_add_ps(_mm256_mul_ps(arge0,breg0), tmp0); 
}

Haswell的

__m256 areg0 = _mm256_set1_ps(a[m]);
for(int i=0; i<n; i++) {      
    __m256 breg0 = _mm256_load_ps(&b[8*i]);
    tmp0 = _mm256_fmadd_ps(arge0, breg0, tmp0);
}

我需要多少次为每个案例展开循环以确保最大吞吐量？

对于使用FMA3的Haswell,我认为答案是每个循环的FLOPS用于沙桥和haswell SSE2/AVX/AVX2.我需要将循环展开10次.

对于Ivy Bridge,我认为它是8.这是我的逻辑.AVX添加的延迟为3,延迟乘以5.Ivy Bridge可以使用不同的端口同时进行一次AVX乘法和一次AVX添加.使用符号m进行乘法,a表示加法,x表示无操作,以及表示部分和的数字(例如m5表示乘以第5部分和)我可以写:

port0:  m1  m2  m3  m4  m5  m6  m7  m8  m1  m2  m3  m4  m5  ... 
port1:   x   x   x   x   x  a1  a2  a3  a4  a5  a6  a7  a8  ...

因此,通过在9个时钟周期后使用8个部分和(4个来自负载,5个来自乘法),我可以在每个时钟周期提交一个AVX负载,一个AVX加法和一个AVX乘法.

我想这意味着在Ivy …

我正在开发大密集矩阵乘法代码.当我对代码进行分析时,它有时会获得我的四个核心系统中大约75%的峰值触发器,而其他时间则大约为36%.在执行代码之间效率不会改变.它要么以75%开始并继续保持效率,要么从36%开始并继续保持这种效率.

我已经将问题追溯到超线程以及我将线程数设置为4而不是默认值8的事实.当我在BIOS中禁用超线程时,我一致地获得大约75%的效率(或者至少我从未看到大幅下降到36%).

在我调用任何并行代码之前omp_set_num_threads(4).我export OMP_NUM_THREADS=4在运行代码之前也尝试过,但它似乎是等效的.

我不想在BIOS中禁用超线程.我想我需要将四个线程绑定到四个核心.我已经测试了一些不同的情况,GOMP_CPU_AFFINITY但到目前为止,我仍然存在效率为36%的问题. 超线程和核心的映射是什么？ 例如,线程0和线程1对应于同一个核心和线程2,线程3对应另一个核心吗？

如何在没有线程迁移的情况下将线程绑定到每个核心,这样我就不必在BIOS中禁用超线程？ 也许我需要研究使用sched_setaffinity？

我当前系统的一些细节:Linux内核3.13,GCC 4.8,Intel Xeon E5-1620(四个物理内核,八个超线程).

编辑:到目前为止,这似乎运作良好

export GOMP_CPU_AFFINITY="0 1 2 3 4 5 6 7"

要么

export GOMP_CPU_AFFINITY="0-7"

编辑:这似乎也运作良好

export OMP_PROC_BIND=true

编辑: 这些选项也很好用(gemm是我的可执行文件的名称)

numactl -C 0,1,2,3 ./gemm

和

taskset -c 0,1,2,3 ./gemm

给定n部分和,可以将log2并行步骤中的所有部分和相加.例如,假设有八个线程与八个部分和:s0, s1, s2, s3, s4, s5, s6, s7.这可以在这样的log2(8) = 3连续步骤中减少;

thread0     thread1    thread2    thread4
s0 += s1    s2 += s3   s4 += s5   s6 +=s7
s0 += s2    s4 += s6
s0 += s4

我想用OpenMP做这个,但我不想使用OpenMP的reduction子句.我想出了一个解决方案,但我认为可以使用OpenMP的task子句找到更好的解决方案.

这比标量加法更通用.让我选择一个更有用的情况:一个数组减少(见这里,这里,并在这里为更多关于阵列减少).

假设我想在阵列上进行数组缩减a.下面是一些代码,它们为每个线程并行填充私有数组.

int bins = 20;
int a[bins];
int **at;  // array of pointers to arrays
for(int i = 0; i<bins; i++) a[i] = 0;
#pragma omp …

fma(a,b,c)相当于a*b+c除了不舍入中间结果.

你能不能给我一些算法的例子,这些算法可以从避免这种舍入中获益？

这并不明显,因为我们避免的乘法后的舍入往往比加法后的舍入更少有问题,而我们没有.

我正在尝试实现2维n体仿真的OpenMP版本.

但是存在一个问题:我假设每个粒子的初始速度和加速度都为零.当颗粒首先聚集在一起时,它们会高速分散,不再聚集.

这似乎与牛顿法不一致,对吧？

有人可以解释为什么会发生这种情况以及如何解决错误？

这是我的代码的一部分:

/* update one frame */
void update() {
    int i, j;

    omp_set_num_threads(8);
    # pragma omp parallel private(j)
    {

    # pragma omp for schedule(static)
        for ( i = 0; i < g_N; ++i ) {
            g_pv[i].a_x = 0.0;
            g_pv[i].a_y = 0.0;
            for ( j = 0; j < g_N; ++j ) {
                if (i == j)
                    continue;
                double r_2 = pow((g_pv[i].pos_x - g_pv[j].pos_x),2) + pow((g_pv[i].pos_y - g_pv[j].pos_y),2);
                g_pv[i].a_x += (-1) * G * g_pv[j].m …

使用GCC 5.3,以下代码符合 -O3 -fma

float mul_add(float a, float b, float c) {
  return a*b + c;
}

生成以下程序集

vfmadd132ss     %xmm1, %xmm2, %xmm0
ret

我注意到GCC -O3已经在GCC 4.8中这样做了.

Clang 3.7带-O3 -mfma产品

vmulss  %xmm1, %xmm0, %xmm0
vaddss  %xmm2, %xmm0, %xmm0
retq

但Clang 3.7与-Ofast -mfmaGCC生成的代码相同-O3 fast.

我很惊讶GCC的确如此,-O3因为从这个答案来看

除非允许使用宽松的浮点模型,否则不允许编译器融合分离的加法和乘法.

这是因为FMA只有一个舍入,而ADD + MUL有两个舍入.因此,编译器将通过融合违反严格的IEEE浮点行为.

但是,从这个链接说

无论FLT_EVAL_METHOD的值如何,任何浮点表达式都可以收缩,即,计算好像所有中间结果都具有无限范围和精度.

所以现在我感到困惑和担忧.

1. GCC是否有理由使用FMA -O3？
2. 融合是否违反了严格的IEEE浮点行为？
3. 如果融合确实违反了IEEE浮点运算,那么GCC的回归__STDC_IEC_559__不是一个矛盾吗？

由于FMA 可以在软件中进行仿真,因此似乎应该有两个用于FMA的编译器开关:一个用于告诉编译器在计算中使用FMA,一个用于告诉编译器硬件具有FMA.

显然,这可以通过选项进行控制-ffp-contract.对于GCC,默认是-ffp-contract=fast和Clang不一样.其他选项例如 …