我想（超级）优化 Heaviside 函数的实现。
我正在研究一种速度特别重要的数值算法（在 Fortran 中）。这多次使用 Heaviside 函数，目前由 signum 内部函数实现，如下所示：
heaviside = 0.5*sign(1,x)+1
我主要对 x 是英特尔处理器上的双精度实数的情况感兴趣。
是否可以开发更有效的 Heaviside 功能实现？也许使用汇编语言、超级优化代码或调用现有的外部库？
问问题
4757 次
1 回答
9
你是故意的heaviside = 0.5*(sign(1,x)+1)吗？无论如何，使用 gcc 4.8.1 fortran 进行的测试表明高性能标记的想法应该是有益的。这里有3种可能性：
heaviside1 - 原始 heaviside2 - 高性能 Mark 的想法 heaviside3 - 另一种变体
function heaviside1 (x)
double precision heaviside1, x
heaviside1 = 0.5 * (sign(1d0,x) + 1)
end
function heaviside2 (x)
double precision heaviside2, x
heaviside2 = sign(0.5d0,x) + 0.5
end
function heaviside3 (x)
double precision heaviside3, x
heaviside3 = 0
if (x .ge. 0) heaviside3 = 1
end
program demo
double precision heaviside1, heaviside2, heaviside3, x, a, b, c
do
x = 0.5 - RAND(0)
a = heaviside1(x)
b = heaviside2(x)
c = heaviside3(x)
print *, "x=", x, "heaviside(x)=", a, b, c
enddo
end
编译后，gcc 会生成这 3 个独立的函数：
<heaviside1_>:
vmovsd xmm0,QWORD PTR [rcx]
vandpd xmm0,xmm0,XMMWORD PTR [rip+0x2d824]
vorpd xmm0,xmm0,XMMWORD PTR [rip+0x2d80c]
vaddsd xmm0,xmm0,QWORD PTR [rip+0x2d7f4]
vmulsd xmm0,xmm0,QWORD PTR [rip+0x2d81c]
ret
<heaviside2_>:
vmovsd xmm0,QWORD PTR [rcx]
vandpd xmm0,xmm0,XMMWORD PTR [rip+0x2d844]
vorpd xmm0,xmm0,XMMWORD PTR [rip+0x2d85c]
vaddsd xmm0,xmm0,QWORD PTR [rip+0x2d844]
ret
<heaviside3_>:
vxorpd xmm0,xmm0,xmm0
vmovsd xmm1,QWORD PTR [rip+0x2d844]
vcmplesd xmm0,xmm0,QWORD PTR [rcx]
vandpd xmm0,xmm1,xmm0
ret
当使用 gcc 编译时，heaviside1 会生成一个可能会减慢执行速度的乘法。heaviside2 消除了乘法。heaviside3 与 heaviside2 具有相同数量的指令，但使用的内存访问次数减少了 2 次。
对于独立功能：
instruction memory reference
count count
heaviside1 6 5
heaviside2 5 4
heaviside3 5 2
这些函数的内联代码避免了对返回指令的需要，并在理想情况下将参数传递到寄存器中，并用所需的常量预加载其他寄存器。确切的结果取决于使用的编译器和调用代码。内联代码的估计：
instruction memory reference
count count
heaviside1 4 0
heaviside2 3 0
heaviside3 2 0
看起来该函数可以由少至两个编译器生成的指令处理：vcmplesd+vandpd。如果参数为负，则第一条指令创建全零掩码，否则创建全零掩码。第二条指令将掩码应用于寄存器常量值 1，以产生结果值 0 或 1。
虽然我没有对这些函数进行基准测试，但看起来重载函数应该不会花费太多执行时间。
---09/23/2013: 添加 x86_64 汇编语言版本和 C 语言基准---
文件函数.s
//----------------------------------------------------------------------------
.intel_syntax noprefix
.text
//-----------------------------------------------------------------------------
// this heaviside function generates its own register constants
// double heaviside_a1 (double arg);
.globl heaviside_a1
heaviside_a1:
mov rax,0x3ff0000000000000
xorpd xmm1,xmm1 # xmm1: constant 0.0
cmplesd xmm1,xmm0 # xmm1: mask (all Fs or all 0s)
movq xmm0,rax # xmm0: constant 1.0
andpd xmm0,xmm1
retq
//-----------------------------------------------------------------------------
// this heaviside function uses register constants passed from caller
// double heaviside_a2 (double arg, double const0, double const1);
.globl heaviside_a2
heaviside_a2:
cmplesd xmm1,xmm0 # xmm1: mask (all Fs or all 0s)
movsd xmm0,xmm2 # xmm0: constant 1.0
andpd xmm0,xmm1
retq
//-----------------------------------------------------------------------------
文件 ctest.c
#define __USE_MINGW_ANSI_STDIO 1
#include <windows.h>
#include <stdio.h>
#include <stdint.h>
// functions.s
double heaviside_a1 (double x);
double heaviside_a2 (double arg, double const0, double const1);
//-----------------------------------------------------------------------------
static double heaviside_c1 (double x)
{
double result = 0;
if (x >= 0) result = 1;
return result;
}
//-----------------------------------------------------------------------------
//
// queryPerformanceCounter - similar to QueryPerformanceCounter, but returns
// count directly.
uint64_t queryPerformanceCounter (void)
{
LARGE_INTEGER int64;
QueryPerformanceCounter (&int64);
return int64.QuadPart;
}
//-----------------------------------------------------------------------------
//
// queryPerformanceFrequency - same as QueryPerformanceFrequency, but returns count direcly.
uint64_t queryPerformanceFrequency (void)
{
LARGE_INTEGER int64;
QueryPerformanceFrequency (&int64);
return int64.QuadPart;
}
//----------------------------------------------------------------------------
//
// lfsr64gpr - left shift galois type lfsr for 64-bit data, general purpose register implementation
//
static uint64_t lfsr64gpr (uint64_t data, uint64_t mask)
{
uint64_t carryOut = data >> 63;
uint64_t maskOrZ = -carryOut;
return (data << 1) ^ (maskOrZ & mask);
}
//---------------------------------------------------------------------------
int runtests (uint64_t pattern, uint64_t mask)
{
uint64_t startCount, elapsed, index, loops = 800000000;
double ns;
double total = 0;
startCount = queryPerformanceCounter ();
for (index = 0; index < loops; index++)
{
double x, result;
pattern = lfsr64gpr (pattern, mask);
x = (double) (int64_t) pattern;
result = heaviside_c1 (x);
total += result;
}
elapsed = queryPerformanceCounter () - startCount;
ns = (double) elapsed / queryPerformanceFrequency () * 1000000000 / loops;
printf ("heaviside_c1: %7.2f ns\n", ns);
startCount = queryPerformanceCounter ();
for (index = 0; index < loops; index++)
{
double x, result;
pattern = lfsr64gpr (pattern, mask);
x = (double) (int64_t) pattern;
result = heaviside_a1 (x);
//printf ("heaviside_a1 (%lf): %lf\n", x, result);
total += result;
}
elapsed = queryPerformanceCounter () - startCount;
ns = (double) elapsed / queryPerformanceFrequency () * 1000000000 / loops;
printf ("heaviside_a1: %7.2f ns\n", ns);
startCount = queryPerformanceCounter ();
for (index = 0; index < loops; index++)
{
double x, result;
const double const0 = 0.0;
const double const1 = 1.0;
pattern = lfsr64gpr (pattern, mask);
x = (double) (int64_t) pattern;
result = heaviside_a2 (x, const0, const1);
//printf ("heaviside_a2 (%lf): %lf\n", x, result);
total += result;
}
elapsed = queryPerformanceCounter () - startCount;
ns = (double) elapsed / queryPerformanceFrequency () * 1000000000 / loops;
printf ("heaviside_a2: %7.2f ns\n", ns);
return total;
}
//---------------------------------------------------------------------------
int main (void)
{
uint64_t mask;
mask = 0xBEFFFFFFFFFFFFFF;
// raise our priority to increase measurement accuracy
SetPriorityClass (GetCurrentProcess (), REALTIME_PRIORITY_CLASS);
printf ("using pseudo-random data\n");
runtests (1, mask);
return 0;
}
//---------------------------------------------------------------------------
mingw64 build command: gcc -Wall -Wextra -O3 -octest.exe ctest.c functions.s
Intel Core i7-2600K 4.0 GHz 的程序输出：
using pseudo-random data
heaviside_c1: 2.24 ns
heaviside_a1: 2.00 ns
heaviside_a2: 2.00 ns
这些时序结果包括执行伪随机参数生成和防止优化器消除其他未使用的 heaviside_c1 本地函数所需的结果汇总代码。
heaviside_c1 来自最初的 fortran 建议，移植到 C。heaviside_a1 是一种汇编语言实现。heaviside_a2 是对汇编语言版本的修改，它使用调用者传递的寄存器常量来避免生成它们的开销。对于我的处理器，基准测试显示传递常量没有任何优势。
汇编语言函数假定 xmm0 返回结果并且 xmm1 和 xmm2 可用作暂存寄存器。这对 Windows 使用的 x64 调用约定有效。对于其他调用约定，应确认此假设。
为了避免内存访问，汇编语言版本期望参数在寄存器 (XMM0) 中按值传递。因为这不是 fortran 的默认设置，所以需要一个特殊的声明。这似乎适用于 gfortran 64 位：
interface
real(c_double) function heaviside_a1(x)
use iso_c_binding, only: c_double
real(c_double), VALUE :: x
end function heaviside_a1
end interface
于 2013-09-16T05:18:15.240 回答