速度的需要 - C++ 与汇编语言
4.76/5 (28投票s)
无法击败 C++ 编译器？
- 下载 CalcDetStackMem.zip - 9.3 KB
- 下载 CalcDet.zip - 9 KB
- 下载 det32asm.zip - 2.9 KB
- 下载 det64asm.zip - 2.9 KB
- 下载 det64masm.zip - 2 KB
- 下载 CSharpDet.zip - 104.5 KB
- 下载 LazarusDet.zip - 1.6 MB
引言
乍一看，我们倾向于相信汇编代码，即通过汇编器从ASM（汇编语言）源代码生成的机器代码，应该比编译代码运行得更快，即通过编译器从高级语言源代码生成的机器代码。原因是ASM与机器代码几乎是一对一（同构）映射（ASM可以称为符号机器代码），而编译器至少需要经过一个两阶段的过程，最后一个阶段是翻译成ASM。
过去确实如此，但编译器，尤其是某些C++编译器，在优化方面变得越来越巧妙，在大多数情况下，只有非常有才华的ASM程序员才能与之媲美或做得更好。一如既往，这条规则有很多例外，但我们在此不予详述。举一个例子，涉及向量化SIMD指令的任务目前最好使用ASM来处理，尽管内在函数（实际上是伪装的ASM）几乎可以同样好地完成工作。
尽管如此，ASM仍然是一门重要的编程语言，不仅因为存在需要尽可能接近硬件的任务，而且最重要的是，它提供了通过高级语言学习和实践不可能获得的见解。
为了说明起见，我将介绍一个程序，该程序将首先通过调用C++函数，然后通过调用链接的汇编ASM对象文件来执行一些耗时的计算。耗时的计算意味着几秒钟，尽管程序可以轻松修改以花费更少或更多的时间。我没有采取任何步骤来改进C++代码；编译器必须尽力处理现有的代码。另一方面，我对ASM源代码进行了一些优化，并将汇编性能提高了30%以上。
背景
我们的程序将使用克莱姆法则计算行列式，这是线性方程组的显式公式，其方程数与未知数相同。
只需考虑行列式的计算是克莱姆法则的核心——这足以满足我们的目的，我们将不再花费更多时间在背后的数学上。
克莱姆法则非常适合计算小矩阵的行列式；然而，对于较大的矩阵尺寸，它在计算上变得非常昂贵，因为其时间复杂度为O(n.n！)——计算13x13矩阵的行列式需要几秒钟，14x14矩阵需要几分钟，而计算15x15矩阵的行列式则需要等待一个小时。
使用代码
1）它是如何编译的？
C++源代码是用Visual Studio 2015编译的。编译的Release版本具有以下设置：
最大化速度，使用流式SIMD扩展2（/arch:SSE2），快速（/fp:fast）浮点模型，禁用安全检查（/GS-），不启用C++异常，Win32禁用（/SAFESEH:NO）。
2）它是如何工作的？
该程序将使用克莱姆法则计算13x13矩阵的行列式。大部分时间花在递归和循环上。之所以需要递归，是因为矩阵的行列式是从低阶矩阵（称为“子式”）的行列式计算的，而子式是通过一次消除顶行和一列来获得的。低阶矩阵的行列式以同样的方式计算，依此类推。在我们的程序中，当达到4x4矩阵时，递归结束，其行列式直接计算。4x4、3x3和2x2阶的矩阵的行列式是直接计算的，不进行任何递归。该程序通过三个级别进行循环以准备子矩阵。对于13x13矩阵，外层循环将执行321,408,685次，中间循环执行1,359,245,472次，最内层循环将执行7,191,246,668次。
4x4子式行列式的计算次数相对少得多，仅为259,459,200次。
3）准备矩阵
矩阵用范围在-50.0到+50.0之间的随机浮点数填充，四舍五入到小数点后1位。
int i, j;
int numRows, numCols;
float *inMat;
srand(time(NULL));
numRows = 13; /*Change as you like for other matrix orders */
numCols = numRows;
inMat = (float*)malloc(numRows*numCols * sizeof(float));
for (i = 0;i < numRows;i++)
for (j = 0;j < numCols;j++)
/* Pointer to matrix, filled as one dimensional array */
inMat[i*numCols + j] = roundf((((float)rand() / (float)(RAND_MAX)) * 100 - 50) * 10) / 10;
printf("%dx%d Matrix\n\n", numRows, numCols);
printf("Input Matrix\n");
for (int i = 0; i < numRows; i++)
{
for (int j = 0; j < numCols; j++)
printf("%.1f\t", inMat[i*numCols + j]);
printf("\n");
}
printf("\n");
4）测量C++执行时间
void C_example(float* inMat, int rows)
{
LARGE_INTEGER frequency;
LARGE_INTEGER t1, t2;
double elapsedTime;
QueryPerformanceFrequency(&frequency);
QueryPerformanceCounter(&t1);
float det = determinant(inMat, rows);
QueryPerformanceCounter(&t2);
elapsedTime = (t2.QuadPart - t1.QuadPart) * 1000.0 / frequency.QuadPart;
printf("\nDeterminant of Matrix using C++ is %.2f.\n", det);
printf("Elapsed time: %f miliseconds\n", elapsedTime);
}
5）在C++中计算行列式
float determinant(float*inMat, int rows)
{
HANDLE heapPtr;
float* minorMat;
int i, j, k, s;
float det = 0.0f, sign = 1.0f;
if (4 == rows)
{
/*
float A11 = inMat[0];
float A12 = inMat[1];
float A13 = inMat[2];
float A14 = inMat[3];
float A21 = inMat[rows];
float A22 = inMat[rows + 1];
float A23 = inMat[rows + 2];
float A24 = inMat[rows + 3];
float A31 = inMat[rows * 2];
float A32 = inMat[rows * 2 + 1];
float A33 = inMat[rows * 2 + 2];
float A34 = inMat[rows * 2 + 3];
float A41 = inMat[rows * 3];
float A42 = inMat[rows * 3 + 1];
float A43 = inMat[rows * 3 + 2];
float A44 = inMat[rows * 3 + 3];
(A13*A24-A14*A23) * ((A31*A42)-(A32*A41)) +
(A14*A22-A12*A24) * ((A31*A43)-(A33*A41)) +
(A13*A22-A12*A23) * ((A34*A41)-(A31*A44)) +
(A11*A24-A14*A21) * ((A32*A43)-(A33*A42)) +
(A11*A23-A13*A21) * ((A34*A42)-(A32*A44)) +
(A11*A22-A12*A21) * ((A33*A44)-(A34*A43))
*/
det = (inMat[2] * inMat[rows + 3] - inMat[3] * inMat[rows + 2]) * ((inMat[rows * 2] * inMat[rows * 3 + 1]) - (inMat[rows * 2 + 1] * inMat[rows * 3])) +
(inMat[3] * inMat[rows + 1] - inMat[1] * inMat[rows + 3]) * ((inMat[rows * 2] * inMat[rows * 3 + 2]) - (inMat[rows * 2 + 2] * inMat[rows * 3])) +
(inMat[2] * inMat[rows + 1] - inMat[1] * inMat[rows + 2]) * ((inMat[rows * 2 + 3] * inMat[rows * 3]) - (inMat[rows * 2] *inMat[rows * 3 + 3])) +
(inMat[0] * inMat[rows + 3] - inMat[3] * inMat[rows]) * ((inMat[rows * 2 + 1] * inMat[rows * 3 + 2]) - (inMat[rows * 2 + 2] * inMat[rows * 3 + 1])) +
(inMat[0] * inMat[rows + 2] - inMat[2] * inMat[rows]) * ((inMat[rows * 2 + 3] * inMat[rows * 3 + 1]) - (inMat[rows * 2 + 1] * inMat[rows * 3 + 3])) +
(inMat[0] * inMat[rows + 1] - inMat[1] * inMat[rows]) * ((inMat[rows * 2 + 2] * inMat[rows * 3 + 3]) - (inMat[rows * 2 + 3] * inMat[rows * 3 + 2]));
}
else if (3 == rows)
{
det = inMat[0] * (inMat[rows + 1] * inMat[rows * 2 + 2] - inMat[rows + 2] * inMat[rows * 2 + 1]) -
inMat[1] * (inMat[rows] * inMat[rows * 2 + 2] - inMat[rows + 2] * inMat[rows * 2]) +
inMat[2] * (inMat[rows] * inMat[rows * 2 + 1] - inMat[rows + 1] * inMat[rows * 2]);
}
else if (2 == rows)
{
det = inMat[0] * inMat[rows + 1] - inMat[1] * inMat[rows];
}
else
{
heapPtr = GetProcessHeap();
minorMat = (float*)HeapAlloc(heapPtr, 0, (rows - 1)*(rows - 1) * sizeof(float));
for (i = 0;i < rows;i++, sign *= -1)
{
for (j = 1;j < rows;j++)
{
for (s = k = 0;k < rows;k++)
{
if (i != k)
{
minorMat[((j - 1)*(rows - 1) + s)] = inMat[(j*rows + k)];
s++;
}
}
}
det += (sign * inMat[i] * determinant(minorMat, rows - 1));
}
HeapFree(heapPtr, 0, minorMat);
}
return det;
}
6）测量汇编ASM的执行时间
void Asm_Example(float* inMat, int rows)
{
LARGE_INTEGER frequency;
LARGE_INTEGER t1, t2;
double elapsedTime;
QueryPerformanceFrequency(&frequency);
QueryPerformanceCounter(&t1);
float det = Determinant(inMat, rows);
QueryPerformanceCounter(&t2);
elapsedTime = (t2.QuadPart - t1.QuadPart) * 1000.0 / frequency.QuadPart;
printf("\nDeterminant of Matrix using Asm is %.2f\n", det);
printf("Elapsed time: %f miliseconds\n", elapsedTime);
}
请注意，Determinant（大写D）函数是一个外部函数。
extern "C"
{
float Determinant(void* intMat, size_t rows);
}
7）ASM源代码。
我们有3种源代码变体：1个用于32位，2个用于64位。
32位版本可以使用MASM直接汇编。
对于64位，我们有一个版本准备好使用JWASM（或HJWASM，它是前者的最新演进）进行汇编。JWASM在64位下的功能至少与MASM在32位下的功能相同。换句话说，MASM 64位缺少其自身的32位功能，但它运行良好，并且有其粉丝不断开发宏来弥补不足（或者说，实际上MASM的意思是“宏汇编器”，所以它在你期望做什么方面是名副其实的！）。无论哪种方式，我还制作了一个可以使用MASM 64位进行汇编的版本。
以下是JWASM/HJWASM的64位版本（32位和MASM 64位版本可以从上面下载）
option casemap:none
option frame:auto
_MM_SHUFFLE MACRO fp3, fp2, fp1, fp0
exitm <( ( ( fp3 ) shl 6) or ( ( fp2 ) shl 4 ) or ( ( fp1 ) shl 2 ) or ( ( fp0 ) ) ) >
ENDM
_PERMUTE_PS MACRO v, c
shufps v, v, c
ENDM
_VECTORSWIZZLE MACRO v, fp0, fp1, fp2, fp3
_PERMUTE_PS v, _MM_SHUFFLE (fp3, fp2, fp1, fp0)
ENDM
_REALSTOXMM MACRO par1, par2, par3, par4
Local xmmValue
.const
align 16
xmmValue dd par1, par2, par3, par4
.code
exitm <xmmValue>
ENDM
T4x4MATRIX struct 16
r0 XMMWORD ?
r1 XMMWORD ?
r2 XMMWORD ?
r3 XMMWORD ?
T4x4MATRIX ends
TFLOAT3X3 struct 4
r0 REAL4 3 dup (?)
r1 REAL4 3 dup (?)
r2 REAL4 3 dup (?)
TFLOAT3X3 ends
HANDLE TYPEDEF PTR
HeapAlloc proto :HANDLE, : DWORD, :QWORD
HeapFree proto :HANDLE, : DWORD, : PTR
GetProcessHeap proto
.const
_MINUSONE REAL4 -1.0
.code
RecursiveDetCalc proc private FRAME uses rsi rdi r12 r13 r14 inMat : ptr, rows : qword
LOCAL heapPtr : ptr
LOCAL det : REAL4
LOCAL sign : REAL4
.if rdx==4
ASSUME rcx : ptr T4x4MATRIX
movups xmm6, [rcx].r2
movups xmm1, [rcx].r3
movaps xmm2, xmm6
_VECTORSWIZZLE xmm2, 1, 0, 0, 0
movaps xmm3, xmm1
_VECTORSWIZZLE xmm3, 2, 2, 1, 1
movaps xmm4, xmm6
_VECTORSWIZZLE xmm4, 1, 0, 0, 0
movaps xmm5, xmm1
_VECTORSWIZZLE xmm5, 3, 3, 3, 2
movaps xmm0, xmm6
_VECTORSWIZZLE xmm0, 2, 2, 1, 1
movaps xmm7, xmm1
_VECTORSWIZZLE xmm7, 3, 3, 3, 2
mulps xmm2, xmm3
movaps xmm8, xmm2
mulps xmm4, xmm5
movaps xmm9, xmm4
mulps xmm0, xmm7
movaps xmm10, xmm0
movaps xmm2, xmm6
_VECTORSWIZZLE xmm2, 2, 2, 1, 1
movaps xmm3, xmm1
_VECTORSWIZZLE xmm3, 1, 0, 0, 0
movaps xmm4, xmm6
_VECTORSWIZZLE xmm4, 3, 3, 3, 2
movaps xmm5, xmm1
_VECTORSWIZZLE xmm5, 1, 0, 0, 0
movaps xmm0, xmm6
_VECTORSWIZZLE xmm0, 3, 3, 3, 2
movaps xmm7, xmm1
_VECTORSWIZZLE xmm7, 2, 2, 1, 1
mulps xmm2, xmm3
movaps xmm3, xmm8
subps xmm3, xmm2
movaps xmm8, xmm3
mulps xmm4, xmm5
subps xmm9, xmm4
mulps xmm0, xmm7
subps xmm10, xmm0
movups xmm6, [rcx].r1
movaps xmm2, xmm6
_VECTORSWIZZLE xmm2, 3, 3, 3, 2
movaps xmm3, xmm6
_VECTORSWIZZLE xmm3, 2, 2, 1, 1
movaps xmm4, xmm6
_VECTORSWIZZLE xmm4, 1, 0, 0, 0
movups xmm0, [rcx].r0
mulps xmm0, XMMWORD ptr _REALSTOXMM(1.0, -1.0, 1.0, -1.0)
mulps xmm2, xmm8
mulps xmm3, xmm9
subps xmm2, xmm3
mulps xmm4, xmm10
addps xmm4, xmm2
mulps xmm0, xmm4
shufps xmm4, xmm0, _MM_SHUFFLE(1,0,0,0)
addps xmm4, xmm0
shufps xmm0, xmm4, _MM_SHUFFLE(0,3,0,0)
addps xmm0, xmm4
_PERMUTE_PS xmm0, _MM_SHUFFLE(2,2,2,2)
ASSUME rcx : nothing
.elseif rdx==3 ; Only used for 3x3 matrixes
ASSUME rcx : ptr TFLOAT3X3
movups xmm3, XMMWORD ptr [rcx].r0
mulps xmm3, XMMWORD ptr _REALSTOXMM(1.0,-1.0,1.0,1.0)
movups xmm0, XMMWORD ptr [rcx].r1+4
movq xmm2, qword ptr [rcx].r2+4
shufps xmm0, xmm2, _MM_SHUFFLE(1,0,1,0)
;calculate first determinant
shufps xmm1, xmm0, _MM_SHUFFLE(0,1,0,0)
mulps xmm0, xmm1
shufps xmm1, xmm0, _MM_SHUFFLE(2,0,0,0)
shufps xmm0, xmm0, _MM_SHUFFLE(3,3,3,3)
shufps xmm1, xmm1, _MM_SHUFFLE(3,3,3,3)
subps xmm0, xmm1
movaps xmm4, xmm0
movups xmm0, XMMWORD ptr [rcx].r1
movups xmm2, XMMWORD ptr [rcx].r1+8
shufps xmm0, xmm2, _MM_SHUFFLE(3,1,2,0)
;calculate second determinant
shufps xmm1, xmm0, _MM_SHUFFLE(0,1,0,0)
mulps xmm0, xmm1
shufps xmm1, xmm0, _MM_SHUFFLE(2,0,0,0)
shufps xmm0, xmm0, _MM_SHUFFLE(3,3,3,3)
shufps xmm1, xmm1, _MM_SHUFFLE(3,3,3,3)
subps xmm0, xmm1
movaps xmm5, xmm0
movups xmm0, XMMWORD ptr [rcx].r1
movq xmm2, qword ptr [rcx].r2
shufps xmm0, xmm2, _MM_SHUFFLE(1,0,1,0)
;calculate third determinant
shufps xmm1, xmm0, _MM_SHUFFLE(0,1,0,0)
mulps xmm0, xmm1
shufps xmm1, xmm0, _MM_SHUFFLE(2,0,0,0)
shufps xmm0, xmm0, _MM_SHUFFLE(3,3,3,3)
shufps xmm1, xmm1, _MM_SHUFFLE(3,3,3,3)
subps xmm0, xmm1
unpckhps xmm4, xmm5
shufps xmm4, xmm0, _MM_SHUFFLE(0,0,1,0)
mulps xmm3, xmm4
movaps xmm0, xmm3
movhlps xmm0, xmm0
addps xmm0, xmm3
shufps xmm3, xmm3, 1
addss xmm0, xmm3
shufps xmm0, xmm0,0
ASSUME rcx : nothing
.elseif rdx==2 ; Only used for 2x2 matrixes
movss xmm0, dword ptr [rcx+12]
movss xmm1, dword ptr [rcx+8]
mulss xmm0, dword ptr [rcx]
mulss xmm1, dword ptr [rcx+4]
subss xmm0, xmm1
.else
mov rsi, rcx
mov rows, rdx
mov r13, rdx
mov det, 0.0
mov sign, 1.0; sign
INVOKE GetProcessHeap
mov heapPtr, rax
mov rax, rows
dec rax
mul rax
shl rax, 2
INVOKE HeapAlloc, heapPtr, 0, rax ; space for new (n-1)x(n-1) matrix
mov rdi, rax
mov r12,0
.while r12<r13
mov r10, 1
.while r10<r13
mov r11,0
mov r9,0
lea rax, [4*r13]
mul r10
mov r14, rax
lea rax, [4*r13-4]
lea rcx, [r10-1]
mul rcx
mov r8, rax
.while r11<r13
.if r12 != r11
lea rcx, [r14 + 4*r11]
mov eax, REAL4 ptr [rsi+rcx]
lea rcx, [r8+4*r9]
mov REAL4 ptr [rdi+rcx], eax
add r9, 1
.endif
add r11, 1
.endw
add r10, 1
.endw
lea rdx, [r13-1]
INVOKE RecursiveDetCalc, rdi, rdx
movss xmm1, dword ptr [rsi+4*r12]
mulss xmm1, sign
mulss xmm0, xmm1
addss xmm0, det
movss xmm1, sign
movss det, xmm0
mulss xmm1, _MINUSONE
movss sign, xmm1
add r12, 1
.endw
INVOKE HeapFree, heapPtr, 0, rdi
movss xmm0, det
.endif
ret
RecursiveDetCalc endp
Determinant proc public FRAME uses xmm6 xmm7 xmm8 xmm9 xmm10 inMat : ptr, rows : qword
; Save xmm registers here and do the calculation in the recursive function.
INVOKE RecursiveDetCalc, rcx, rdx
ret
Determinant endp
end
8）编译64位Release版本并运行后，您将获得类似下面图像的结果。
忽略行列式值的差异，这是由于不同计算变体下的舍入误差造成的。
为什么ASM更慢？
如前所述，我对ASM代码进行了许多优化，但仍不足以击败C++编译。SIMD部分，特别是4x4子式行列式计算，在ASM代码中更快，但它们在整个过程的总耗时中所占比例不大。
我看到了C++编译器生成的ASM列表，有些部分简直令人难以置信——没有人会相信人类会那样编码（如果他这样做了，代码将几乎无法维护）。编译器以自动化方式使用了各种技巧——难以超越。它了解管道和预测分支的所有工作原理，它重新排序指令，进行循环切片并使用缓存无关算法。
更新1
在撰写本文时，与其它编程语言或编译器相比，ASM编程仍然能产生更快お代码。现在断言ASM编程已经过时还为时过早。我选择了两种流行的语言和环境，并复制了C++代码。
Lazarus/Free Pascal
Lazarus是一个IDE，具有Delphi（旧版本）的大部分RAD功能。其集成的Free Pascal编译器似乎能生成比Delphi更快的64位可执行文件（尽管32位稍慢）。
与Visual Studio C++编译器相比，它在我们测试中落后了，花费了将近30%的时间。
正如您所知，对于大多数应用程序来说，这不成问题。
您可以从上面下载测试的完整Free Pascal源代码。有几点需要提及。
- 对于32位，您需要使用JWASM/HJWASM（使用-zt0开关）而不是MASM来汇编ASM，以删除API调用的STDCALL名称修饰（例如_HeapAlloc@12）。Free Pascal（以及Delphi）无法处理这一点。
- 与Delphi不同，Free Pascal在链接时无法在ASM代码中处理API引用，而无需包装函数。这是一个不便之处，尽管对整体性能没有实质性影响。
Visual Studio/C#
C#是微软发明的一个不错的编程语言。在其当前实现中，它与庞大的.Net Framework和庞大的公共语言运行时（Common Language Runtime）紧密集成。
尽管C#的性能多年来一直在不断提高，但在我们的测试中，JIT编译器在速度方面并未产生惊人的结果。我确实相信C#源代码可以针对速度进行改进，这一点毫无疑问，但C++和Free Pascal也可以。因此，不要急于得出结论。
您可以从上面下载测试的完整C#源代码。照常，ASM通过C++ CLI类调用，否则我们需要另一个DLL。
- 我们使用了VS 2015。您必须为x64或x86 CPU构建，而不是为AnyCpu构建。
- 您也可以通过简单的批处理文件将其构建为单个x86或x64可执行文件（完全没有DLL），如您的下载中的CSharpSingleExe目录所示。这是基于我们的文章在独立的64位exe中混合.NET和汇编语言。
更新2
我决定在不使用堆内存的情况下，仅使用堆栈内存来复制C++与ASM的测试。
堆分配的影响是真实的，但并非首要（约16-18%），并且对C++和ASM的影响程度（几乎）相同。结论仍然是相同的——正如预期的那样——但测试现在更加透明。更改后的代码可以从上面下载（CalcDetStackMem.zip）。
更新历史
- 2017年4月18日 - 初始文章
- 2017年4月28日 - 包含Free Pascal和C#的等效测试。x86和x64的ASM进行了小型修复。
- 2017年5月5日 - C++与ASM测试使用堆栈而不是堆内存分配。