- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
问题
我一直在研究 HPC,特别是使用矩阵乘法作为我的项目(请参阅我在个人资料中的其他帖子)。我在这些方面取得了不错的成绩,但还不够好。我退后一步,看看我可以用点积计算做多好。
点积与矩阵乘法
点积更简单,可以让我测试 HPC 概念,而无需处理打包和其他相关问题。缓存阻塞仍然是一个问题,这构成了我的第二个问题。
算法
乘法 n
两个中的对应元素double
数组 A
和 B
并将它们相加。一个 double
组装中的点积只是一系列movapd
, mulpd
, addpd
.以巧妙的方式展开和排列,可以有 movapd
组/mulpd
/addpd
在不同的 xmm
上运行寄存器,因此是独立的,优化流水线。当然,事实证明这并不重要,因为我的 CPU 是乱序执行的。另请注意,重新安排需要剥离最后一次迭代。
其他假设
我不是在为一般的点积编写代码。该代码适用于特定尺寸,我不处理边缘情况。这只是为了测试 HPC 概念并看看我可以达到什么类型的 CPU 使用率。
结果
编译 gcc -std=c99 -O2 -m32 -mincoming-stack-boundary=2 -msse3 -mfpmath=sse,387 -masm=intel
.我在一台与平时不同的电脑上。这台电脑有 i5 540m
可以获得2.8 GHz * 4 FLOPS/cycle/core = 11.2 GFLOPS/s per core
在两步 Intel Turbo Boost 之后(两个内核现在都处于开启状态,所以它只能进行 2 步……如果我关闭一个内核,则可以进行 4 步增强)。当设置为使用一个线程运行时,32 位 LINPACK 的速度约为 9.5 GFLOPS/s。
N Total Gflops/s Residual
256 5.580521 1.421085e-014
384 5.734344 -2.842171e-014
512 5.791168 0.000000e+000
640 5.821629 0.000000e+000
768 5.814255 2.842171e-014
896 5.807132 0.000000e+000
1024 5.817208 -1.421085e-013
1152 5.805388 0.000000e+000
1280 5.830746 -5.684342e-014
1408 5.881937 -5.684342e-014
1536 5.872159 -1.705303e-013
1664 5.881536 5.684342e-014
1792 5.906261 -2.842171e-013
1920 5.477966 2.273737e-013
2048 5.620931 0.000000e+000
2176 3.998713 1.136868e-013
2304 3.370095 -3.410605e-013
2432 3.371386 -3.410605e-013
n > 2048
,您可以看到性能下降。这是因为我的L1缓存是32KB,当
n = 2048
和
A
和
B
是
double
,它们占用了整个缓存。任何更大的它们都是从内存中流式传输的。
#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <x86intrin.h>
#include <math.h>
#include <omp.h>
#include <stdint.h>
#include <windows.h>
// computes 8 dot products
#define KERNEL(address) \
"movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \
"mulpd xmm7, XMMWORD PTR [edx+48+"#address"] \n\t" \
"addpd xmm2, xmm6 \n\t" \
"movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \
"mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \
"addpd xmm3, xmm7 \n\t" \
"movapd xmm6, XMMWORD PTR [eax+96+"#address"] \n\t" \
"mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \
"addpd xmm0, xmm4 \n\t" \
"movapd xmm7, XMMWORD PTR [eax+112+"#address"] \n\t" \
"mulpd xmm6, XMMWORD PTR [edx+96+"#address"] \n\t" \
"addpd xmm1, xmm5 \n\t"
#define PEELED(address) \
"movapd xmm4, XMMWORD PTR [eax+"#address"] \n\t" \
"mulpd xmm7, [edx+48+"#address"] \n\t" \
"addpd xmm2, xmm6 \n\t" \
"movapd xmm5, XMMWORD PTR [eax+16+"#address"] \n\t" \
"mulpd xmm4, XMMWORD PTR [edx+"#address"] \n\t" \
"addpd xmm3, xmm7 \n\t" \
"mulpd xmm5, XMMWORD PTR [edx+16+"#address"] \n\t" \
"addpd xmm0, xmm4 \n\t" \
"addpd xmm1, xmm5 \n\t"
inline double
__attribute__ ((gnu_inline))
__attribute__ ((aligned(64))) ddot_ref(
int n,
const double* restrict A,
const double* restrict B)
{
double sum0 = 0.0;
double sum1 = 0.0;
double sum2 = 0.0;
double sum3 = 0.0;
double sum;
for(int i = 0; i < n; i+=4) {
sum0 += *(A + i ) * *(B + i );
sum1 += *(A + i+1) * *(B + i+1);
sum2 += *(A + i+2) * *(B + i+2);
sum3 += *(A + i+3) * *(B + i+3);
}
sum = sum0+sum1+sum2+sum3;
return(sum);
}
inline double
__attribute__ ((gnu_inline))
__attribute__ ((aligned(64))) ddot_asm
( int n,
const double* restrict A,
const double* restrict B)
{
double sum;
__asm__ __volatile__
(
"mov eax, %[A] \n\t"
"mov edx, %[B] \n\t"
"mov ecx, %[n] \n\t"
"pxor xmm0, xmm0 \n\t"
"pxor xmm1, xmm1 \n\t"
"pxor xmm2, xmm2 \n\t"
"pxor xmm3, xmm3 \n\t"
"movapd xmm6, XMMWORD PTR [eax+32] \n\t"
"movapd xmm7, XMMWORD PTR [eax+48] \n\t"
"mulpd xmm6, XMMWORD PTR [edx+32] \n\t"
"sar ecx, 7 \n\t"
"sub ecx, 1 \n\t" // peel
"L%=: \n\t"
KERNEL(64 * 0)
KERNEL(64 * 1)
KERNEL(64 * 2)
KERNEL(64 * 3)
KERNEL(64 * 4)
KERNEL(64 * 5)
KERNEL(64 * 6)
KERNEL(64 * 7)
KERNEL(64 * 8)
KERNEL(64 * 9)
KERNEL(64 * 10)
KERNEL(64 * 11)
KERNEL(64 * 12)
KERNEL(64 * 13)
KERNEL(64 * 14)
KERNEL(64 * 15)
"lea eax, [eax+1024] \n\t"
"lea edx, [edx+1024] \n\t"
" \n\t"
"dec ecx \n\t"
"jnz L%= \n\t" // end loop
" \n\t"
KERNEL(64 * 0)
KERNEL(64 * 1)
KERNEL(64 * 2)
KERNEL(64 * 3)
KERNEL(64 * 4)
KERNEL(64 * 5)
KERNEL(64 * 6)
KERNEL(64 * 7)
KERNEL(64 * 8)
KERNEL(64 * 9)
KERNEL(64 * 10)
KERNEL(64 * 11)
KERNEL(64 * 12)
KERNEL(64 * 13)
KERNEL(64 * 14)
PEELED(64 * 15)
" \n\t"
"addpd xmm0, xmm1 \n\t" // summing result
"addpd xmm2, xmm3 \n\t"
"addpd xmm0, xmm2 \n\t" // cascading add
"movapd xmm1, xmm0 \n\t" // copy xmm0
"shufpd xmm1, xmm0, 0x03 \n\t" // shuffle
"addsd xmm0, xmm1 \n\t" // add low qword
"movsd %[sum], xmm0 \n\t" // mov low qw to sum
: // outputs
[sum] "=m" (sum)
: // inputs
[A] "m" (A),
[B] "m" (B),
[n] "m" (n)
: //register clobber
"memory",
"eax","ecx","edx","edi",
"xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7"
);
return(sum);
}
int main()
{
// timers
LARGE_INTEGER frequency, time1, time2;
double time3;
QueryPerformanceFrequency(&frequency);
// clock_t time1, time2;
double gflops;
int nmax = 4096;
int trials = 1e7;
double sum, residual;
FILE *f = fopen("soddot.txt","w+");
printf("%16s %16s %16s\n","N","Total Gflops/s","Residual");
fprintf(f,"%16s %16s %16s\n","N","Total Gflops/s","Residual");
for(int n = 256; n <= nmax; n += 128 ) {
double* A = NULL;
double* B = NULL;
A = _mm_malloc(n*sizeof(*A), 64); if (!A) {printf("A failed\n"); return(1);}
B = _mm_malloc(n*sizeof(*B), 64); if (!B) {printf("B failed\n"); return(1);}
srand(time(NULL));
// create arrays
for(int i = 0; i < n; ++i) {
*(A + i) = (double) rand()/RAND_MAX;
*(B + i) = (double) rand()/RAND_MAX;
}
// warmup
sum = ddot_asm(n,A,B);
QueryPerformanceCounter(&time1);
// time1 = clock();
for (int count = 0; count < trials; count++){
// sum = ddot_ref(n,A,B);
sum = ddot_asm(n,A,B);
}
QueryPerformanceCounter(&time2);
time3 = (double)(time2.QuadPart - time1.QuadPart) / frequency.QuadPart;
// time3 = (double) (clock() - time1)/CLOCKS_PER_SEC;
gflops = (double) (2.0*n*trials)/time3/1.0e9;
residual = ddot_ref(n,A,B) - sum;
printf("%16d %16f %16e\n",n,gflops,residual);
fprintf(f,"%16d %16f %16e\n",n,gflops,residual);
_mm_free(A);
_mm_free(B);
}
fclose(f);
return(0); // successful completion
}
sum += a[i]*b[i]
.
sum
必须初始化为
0
在第一次迭代之前。向量化,你可以一次做 2 个求和,必须在最后求和:
[sum0 sum1] = [a[i] a[i+1]]*[b[i] b[i+1]]
,
sum = sum0 + sum1
.在(英特尔)汇编中,这是 3 个步骤(初始化之后):
pxor xmm0, xmm0 // accumulator [sum0 sum1] = [0 0]
movapd xmm1, XMMWORD PTR [eax] // load [a[i] a[i+1]] into xmm1
mulpd xmm1, XMMWORD PTR [edx] // xmm1 = xmm1 * [b[i] b[i+1]]
addpd xmm0, xmm1 // xmm0 = xmm0 + xmm1
xmm
来获得更好的性能。可用的寄存器(32 位模式下的 8 个寄存器)。因此,如果您将其展开 4 次,则可以使用所有 8 个寄存器
xmm0
通过
xmm7
.您将有 4 个累加器和 4 个寄存器用于存储
movapd
的结果。和
addpd
.同样,编译器可以想出这个。真正的思考部分是试图想出一种管道代码的方法,即让 MOV/MUL/ADD 组中的每条指令都在不同的寄存器上运行,以便所有 3 条指令同时执行(通常情况下大多数 CPU)。这就是你击败编译器的方式。因此,您必须对 4x 展开代码进行模式化才能做到这一点,这可能需要提前加载 vector 并剥离第一次或最后一次迭代。这是什么
KERNEL(address)
是。为方便起见,我制作了一个 4x 展开流水线代码的宏。这样我只需更改
address
就可以轻松地将其展开为 4 的倍数.每个
KERNEL
计算 8 个点积。
最佳答案
要回答您的整体问题,您无法使用点积实现最佳性能。
问题是您的 CPU 可以在每个时钟周期执行一次 128 位加载,而要执行点积,您需要每个时钟周期执行两次 128 位加载。
但它比大 n 更糟糕。你的第二个问题的答案是点积是内存限制而不是计算限制,因此它不能并行化具有快速内核的大 n 。在这里解释得更好 why-vectorizing-the-loop-does-not-have-performance-improvement .这是与快速内核并行化的一个大问题。我花了一段时间才弄明白这一点,但学习非常重要。
实际上,很少有基本算法可以充分受益于快速内核上的并行化。就 BLAS 算法而言,只有 Level-3 算法 (O(n^3))(例如矩阵乘法)才能真正受益于并行化。在慢速内核上情况更好,例如使用 GPU 和 Xeon Phi,因为内存速度和核心速度之间的差异要小得多。
如果您想找到一种算法,它可以接近小 n 的峰值触发器,请尝试例如标量 * vector 或标量 * vector 的总和。第一种情况应该在每个时钟周期执行一次加载、一次乘法和一次存储,而第二种情况应该在每个时钟周期进行一次乘法、一次加法和一次加载。
我在 Knoppix 7.3 32 位的 Core 2 Duo P9600@2.67GHz 上测试了以下代码。 我得到了标量乘积的大约 75% 的峰值和标量乘积总和的峰值的 75%。 标量乘积的触发器/周期为 2,标量乘积的总和为 4。
编译 g++ -msse2 -O3 -fopenmp foo.cpp -ffast-math
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <x86intrin.h>
void scalar_product(double * __restrict a, int n) {
a = (double*)__builtin_assume_aligned (a, 64);
double k = 3.14159;
for(int i=0; i<n; i++) {
a[i] = k*a[i];
}
}
void scalar_product_SSE(double * __restrict a, int n) {
a = (double*)__builtin_assume_aligned (a, 64);
__m128d k = _mm_set1_pd(3.14159);
for(int i=0; i<n; i+=8) {
__m128d t1 = _mm_load_pd(&a[i+0]);
_mm_store_pd(&a[i],_mm_mul_pd(k,t1));
__m128d t2 = _mm_load_pd(&a[i+2]);
_mm_store_pd(&a[i+2],_mm_mul_pd(k,t2));
__m128d t3 = _mm_load_pd(&a[i+4]);
_mm_store_pd(&a[i+4],_mm_mul_pd(k,t3));
__m128d t4 = _mm_load_pd(&a[i+6]);
_mm_store_pd(&a[i+6],_mm_mul_pd(k,t4));
}
}
double scalar_sum(double * __restrict a, int n) {
a = (double*)__builtin_assume_aligned (a, 64);
double sum = 0.0;
double k = 3.14159;
for(int i=0; i<n; i++) {
sum += k*a[i];
}
return sum;
}
double scalar_sum_SSE(double * __restrict a, int n) {
a = (double*)__builtin_assume_aligned (a, 64);
__m128d sum1 = _mm_setzero_pd();
__m128d sum2 = _mm_setzero_pd();
__m128d sum3 = _mm_setzero_pd();
__m128d sum4 = _mm_setzero_pd();
__m128d k = _mm_set1_pd(3.14159);
for(int i=0; i<n; i+=8) {
__m128d t1 = _mm_load_pd(&a[i+0]);
sum1 = _mm_add_pd(_mm_mul_pd(k,t1),sum1);
__m128d t2 = _mm_load_pd(&a[i+2]);
sum2 = _mm_add_pd(_mm_mul_pd(k,t2),sum2);
__m128d t3 = _mm_load_pd(&a[i+4]);
sum3 = _mm_add_pd(_mm_mul_pd(k,t3),sum3);
__m128d t4 = _mm_load_pd(&a[i+6]);
sum4 = _mm_add_pd(_mm_mul_pd(k,t4),sum4);
}
double tmp[8];
_mm_storeu_pd(&tmp[0],sum1);
_mm_storeu_pd(&tmp[2],sum2);
_mm_storeu_pd(&tmp[4],sum3);
_mm_storeu_pd(&tmp[6],sum4);
double sum = 0;
for(int i=0; i<8; i++) sum+=tmp[i];
return sum;
}
int main() {
//_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
//_mm_setcsr(_mm_getcsr() | 0x8040);
double dtime, peak, flops, sum;
int repeat = 1<<18;
const int n = 2048;
double *a = (double*)_mm_malloc(sizeof(double)*n,64);
double *b = (double*)_mm_malloc(sizeof(double)*n,64);
for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;
dtime = omp_get_wtime();
for(int r=0; r<repeat; r++) {
scalar_product_SSE(a,n);
}
dtime = omp_get_wtime() - dtime;
peak = 2*2.67;
flops = 1.0*n/dtime*1E-9*repeat;
printf("time %f, %f, %f\n", dtime,flops, flops/peak);
//for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;
//sum = 0.0;
dtime = omp_get_wtime();
for(int r=0; r<repeat; r++) {
scalar_sum_SSE(a,n);
}
dtime = omp_get_wtime() - dtime;
peak = 2*2*2.67;
flops = 2.0*n/dtime*1E-9*repeat;
printf("time %f, %f, %f\n", dtime,flops, flops/peak);
//printf("sum %f\n", sum);
}
关于c - 如何使用 Dot Product 获得峰值 CPU 性能?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/24102103/
我在网上搜索但没有找到任何合适的文章解释如何使用 javascript 使用 WCF 服务,尤其是 WebScriptEndpoint。 任何人都可以对此给出任何指导吗? 谢谢 最佳答案 这是一篇关于
我正在编写一个将运行 Linux 命令的 C 程序,例如: cat/etc/passwd | grep 列表 |剪切-c 1-5 我没有任何结果 *这里 parent 等待第一个 child (chi
所以我正在尝试处理文件上传,然后将该文件作为二进制文件存储到数据库中。在我存储它之后,我尝试在给定的 URL 上提供文件。我似乎找不到适合这里的方法。我需要使用数据库,因为我使用 Google 应用引
我正在尝试制作一个宏,将下面的公式添加到单元格中,然后将其拖到整个列中并在 H 列中复制相同的公式 我想在 F 和 H 列中输入公式的数据 Range("F1").formula = "=IF(ISE
问题类似于this one ,但我想使用 OperatorPrecedenceParser 解析带有函数应用程序的表达式在 FParsec . 这是我的 AST: type Expression =
我想通过使用 sequelize 和 node.js 将这个查询更改为代码取决于在哪里 select COUNT(gender) as genderCount from customers where
我正在使用GNU bash,版本5.0.3(1)-发行版(x86_64-pc-linux-gnu),我想知道为什么简单的赋值语句会出现语法错误: #/bin/bash var1=/tmp
这里,为什么我的代码在 IE 中不起作用。我的代码适用于所有浏览器。没有问题。但是当我在 IE 上运行我的项目时,它发现错误。 而且我的 jquery 类和 insertadjacentHTMl 也不
我正在尝试更改标签的innerHTML。我无权访问该表单,因此无法编辑 HTML。标签具有的唯一标识符是“for”属性。 这是输入和标签的结构:
我有一个页面,我可以在其中返回用户帖子,可以使用一些 jquery 代码对这些帖子进行即时评论,在发布新评论后,我在帖子下插入新评论以及删除 按钮。问题是 Delete 按钮在新插入的元素上不起作用,
我有一个大约有 20 列的“管道分隔”文件。我只想使用 sha1sum 散列第一列,它是一个数字,如帐号,并按原样返回其余列。 使用 awk 或 sed 执行此操作的最佳方法是什么? Accounti
我需要将以下内容插入到我的表中...我的用户表有五列 id、用户名、密码、名称、条目。 (我还没有提交任何东西到条目中,我稍后会使用 php 来做)但由于某种原因我不断收到这个错误:#1054 - U
所以我试图有一个输入字段,我可以在其中输入任何字符,但然后将输入的值小写,删除任何非字母数字字符,留下“。”而不是空格。 例如,如果我输入: 地球的 70% 是水,-!*#$^^ & 30% 土地 输
我正在尝试做一些我认为非常简单的事情,但出于某种原因我没有得到想要的结果?我是 javascript 的新手,但对 java 有经验,所以我相信我没有使用某种正确的规则。 这是一个获取输入值、检查选择
我想使用 angularjs 从 mysql 数据库加载数据。 这就是应用程序的工作原理;用户登录,他们的用户名存储在 cookie 中。该用户名显示在主页上 我想获取这个值并通过 angularjs
我正在使用 autoLayout,我想在 UITableViewCell 上放置一个 UIlabel,它应该始终位于单元格的右侧和右侧的中心。 这就是我想要实现的目标 所以在这里你可以看到我正在谈论的
我需要与 MySql 等效的 elasticsearch 查询。我的 sql 查询: SELECT DISTINCT t.product_id AS id FROM tbl_sup_price t
我正在实现代码以使用 JSON。 func setup() { if let flickrURL = NSURL(string: "https://api.flickr.com/
我尝试使用for循环声明变量,然后测试cols和rols是否相同。如果是,它将运行递归函数。但是,我在 javascript 中执行 do 时遇到问题。有人可以帮忙吗? 现在,在比较 col.1 和
我举了一个我正在处理的问题的简短示例。 HTML代码: 1 2 3 CSS 代码: .BB a:hover{ color: #000; } .BB > li:after {
我是一名优秀的程序员,十分优秀!