2

SIMD（SSE2など）を使用して次の関数を最適化したい：

int64_t fun(int64_t N, int size, int* p)
{
    int64_t sum = 0;
    for(int i=1; i<size; i++)
       sum += (N/i)*p[i];
    return sum;
}

これは、必要な命令がそこにないことを除いて、非常にベクトル化可能なタスクのように思えます...

N は非常に大きく (10^12 から 10^18)、size~sqrt(N) であると仮定できます。また、p は -1、0、および 1 の値しかとれないと仮定することもできます。したがって、実際の乗算は必要ありません。(N/i)*p[i] は、何らかの方法で N/i を計算できれば、4 つの命令 (pcmpgt、pxor、psub、pand) で実行できます。

4

4 に答える 4

2

の導関数は1/xです。-1/x^2つまり、 asxが大きくなると、 になりN/x==N/(x + 1)ます。

N/x(その値と呼びましょう) の既知の値についてr、次の値を決定できます(その値を次のようxに呼びましょう:x'N/x'<r

x'= N/(r - 1)

そして、整数を扱っているので：

x'= ceiling(N/(r - 1))

したがって、ループは次のようになります。

int64_t sum = 0;   
int i=1; 
int r= N;
while (i<size)
{
    int s= (N + r - 1 - 1)/(r - 1);

    while (i<s && i<size)
    {
        sum += (r)*p[i];
        ++i;
    }

    r= N/s;
}
return sum;

N が十分に大きい場合、 に対して同一の値が多数実行されN/iます。確かに、注意しないと、0 で除算することになります。

于 2010-11-05T00:06:58.547 に答える
2

これは、そのコードのベクトル化に最も近いものです。私はそれが速くなることを本当に期待していません。私はちょうど SIMD コードを書こうとしていました。

#include <stdint.h>

int64_t fun(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    for(i=1; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

typedef int64_t v2sl __attribute__ ((vector_size (2*sizeof(int64_t))));

int64_t fun_simd(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    v2sl v_2 = { 2, 2 };
    v2sl v_N = { N, N };
    v2sl v_i = { 1, 2 };
    union { v2sl v; int64_t a[2]; } v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) {
        v2sl v_p = { p[i], p[i+1] };
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    }
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

typedef double v2df __attribute__ ((vector_size (2*sizeof(double))));

int64_t fun_simd_double(int64_t N, int size, const int* p)
{
    int64_t sum = 0;
    int i;
    v2df v_2 = { 2, 2 };
    v2df v_N = { N, N };
    v2df v_i = { 1, 2 };
    union { v2df v; double a[2]; } v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) {
        v2df v_p = { p[i], p[i+1] };
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    }
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) {
        sum += (N/i)*p[i];
    }
    return sum;
}

#include <stdio.h>

static const int test_array[] = {
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0
};
#define test_array_len (sizeof(test_array)/sizeof(int))

#define big_N (1024 * 1024 * 1024)

int main(int argc, char *argv[]) {
    int64_t res1;
    int64_t res2;
    int64_t res3;
    v2sl a = { 123, 456 };
    v2sl b = { 100, 200 };
    union { v2sl v; int64_t a[2]; } tmp;

    a = a + b;
    tmp.v = a;
    printf("a = { %ld, %ld }\n", tmp.a[0], tmp.a[1]);

    printf("test_array size = %zd\n", test_array_len);

    res1 = fun(big_N, test_array_len, test_array);
    printf("fun() = %ld\n", res1);

    res2 = fun_simd(big_N, test_array_len, test_array);
    printf("fun_simd() = %ld\n", res2);

    res3 = fun_simd_double(big_N, test_array_len, test_array);
    printf("fun_simd_double() = %ld\n", res3);

    return 0;
}
于 2010-11-01T13:18:35.170 に答える
1

精度要件に応じて、単精度または倍精度の浮動小数点 SIMD 演算でこれを行うことをお勧めします。int から float または double への変換は、SSE を使用すると比較的高速です。

于 2010-10-29T07:09:02.550 に答える
1

除算の計算にコストが集中します。SSE2 には整数除算用のオペコードがないため、除算アルゴリズムを自分で少しずつ実装する必要があります。努力する価値があるとは思いません.SSE2では2つのインスタンスを並行して実行できます（64ビットの数値を使用し、SSE2レジスタは128ビットです）が、手作りの除算アルゴリズムは少なくともCPUidivオペコードの 2 倍遅い。

(ちなみに、32 ビット モードと 64 ビット モードのどちらでコンパイルしますか?後者は 64 ビット整数でより快適です。)

部門の全体数を減らすことは、より有望な方法のように見えます。正の整数xおよびyの場合、floor(x/(2y)) = floor(floor(x/y)/2)であることに注意してください。C 用語では、N/i(切り捨て除算) を計算したら、それを 1 ビット右にシフトするだけで を取得できN/(2*i)ます。適切に使用すると、これにより分割の半分がほぼ無料になります (「適切に」にはp[i]、キャッシュに大混乱をもたらさない方法で数十億の値にアクセスすることも含まれるため、非常に簡単ではないようです)。

于 2010-11-01T12:19:42.807 に答える