2012-10-18 19 views
11

Opracowuję algorytm kryptograficzny na GPU i obecnie utknąłem z algorytmem do wykonywania dużych dodatków całkowitych. Duże liczby całkowite są reprezentowane w zwykły sposób jako grupa 32-bitowych słów.duża liczba całkowita dodana z CUDA

Na przykład możemy użyć jednego wątku, aby dodać dwa słowa 32-bitowe. Dla uproszczenia przyjmijmy , że liczby do dodania mają tę samą długość i liczbę wątków na blok == liczba słów. Następnie:

__global__ void add_kernel(int *C, const int *A, const int *B) { 
    int x = A[threadIdx.x]; 
    int y = B[threadIdx.x]; 
    int z = x + y; 
    int carry = (z < x); 
    /** do carry propagation in parallel somehow ? */ 
    ............ 

    z = z + newcarry; // update the resulting words after carry propagation 
    C[threadIdx.x] = z; 
} 

Jestem pewien, że istnieje sposób, aby zrobić propagację carry przez jakiś skomplikowanej procedury redukcji, ale nie mógł zrozumieć to ..

musiałem spojrzeć na CUDA thrust extensions ale duży pakiet całkowitą wydaje jeszcze nie wdrożone. Być może ktoś może dać mi wskazówkę, jak to zrobić w CUDA?

+2

GPU może obsługiwać do 64 bitów (Dawno) bezpośrednio. Jedno podejście do 128-bitowego jest opisane w [to SO pytanie/odpowiedź] (http://stackoverflow.com/questions/6162140/128-bit-integer-on-cuda). –

+0

Myślę, że to, czego oczekujesz od CUDA, można osiągnąć dzięki technikom C. Dlatego też powtórzyłem pytanie w "C". Mam nadzieję uzyskać dobrą odpowiedź od ekspertów C. – ahmad

+0

Tak, można również zaprogramować długie dodawanie całkowite przy użyciu tylko konstruktów wysokiego poziomu C (w przeciwieństwie do zespołu linii PXT w CUDA), ale wymagałoby to znacznie więcej instrukcji, jak wskazałem w tej odpowiedzi: http: // stackoverflow .pl/questions/12448549/is-inline-ptx-assembly-code-powerful/12453534 # 12453534 – njuffa

Odpowiedz

8

Masz rację, przeniesienia propagacji można dokonać za pomocą obliczenia sumy prefiksów, ale trochę trudniej jest zdefiniować funkcję binarną dla tej operacji i udowodnić, że jest ona asocjacyjna (potrzebna dla sumy przedrostka równoległego). W rzeczywistości ten algorytm jest używany (teoretycznie) w Carry-lookahead adder.

Załóżmy, że mamy dwie duże liczby całkowite a [0..n-1] i b [0..n-1]. Następnie obliczyć (i = 0..n-1):

s[i] = a[i] + b[i]l; 
carryin[i] = (s[i] < a[i]); 

określamy dwie funkcje:

generate[i] = carryin[i]; 
propagate[i] = (s[i] == 0xffffffff); 

z bardzo intuicyjny sposób: generowanie [i] == 1 oznacza, że ​​przeniesienia jest generowany na pozycji podczas propagacji [i] == 1 oznacza, że ​​przeniesienie będzie propagowane z pozycji (i - 1) do (i + 1). Naszym celem jest obliczenie funkcji [0..n-1] używanej do aktualizacji sumy wynikowej [0..n-1]. Carryout mogą być obliczane rekurencyjnie w następujący sposób:

carryout[i] = generate[i] OR (propagate[i] AND carryout[i-1]) 
carryout[0] = 0 

tutaj Carryout [i] == 1, w przypadku przenoszenia jest generowany w pozycji I lub generowany jest niekiedy wcześniej, propagowana do położenia I. Wreszcie, musimy zaktualizować wynikowy Suma:

s[i] = s[i] + carryout[i-1]; for i = 1..n-1 
carry = carryout[n-1]; 

Teraz jest to dość łatwe do udowodnienia, że ​​rzeczywiście jest funkcja Carryout binarny asocjacyjne i stąd równoległe obliczenie sumy przedrostek odnosi. Aby zaimplementować to na CUDA, możemy połączyć „generować” obie flagi i „propagować” w pojedynczej zmiennej, gdyż wykluczają się wzajemnie, tj:

cy[i] = (s[i] == -1u ? -1u : 0) | carryin[i]; 

Innymi słowy,

cy[i] = 0xffffffff if propagate[i] 
cy[i] = 1   if generate[i] 
cy[u] = 0   otherwise 

Następnie można sprawdzić, że następująca formuła oblicza sumę prefiks dla funkcji carryout:

cy[i] = max((int)cy[i], (int)cy[k]) & cy[i]; 

dla wszystkich k < i. Poniższy przykładowy kod pokazuje duży dodatek dla liczb całkowitych z 2048 liczb. Użyłem tu bloki CUDA 512 wątków:

// add & output carry flag 
#define UADDO(c, a, b) \ 
    asm volatile("add.cc.u32 %0, %1, %2;" : "=r"(c) : "r"(a) , "r"(b)); 
// add with carry & output carry flag 
#define UADDC(c, a, b) \ 
    asm volatile("addc.cc.u32 %0, %1, %2;" : "=r"(c) : "r"(a) , "r"(b)); 

#define WS 32 

__global__ void bignum_add(unsigned *g_R, const unsigned *g_A,const unsigned *g_B) { 

extern __shared__ unsigned shared[]; 
unsigned *r = shared; 

const unsigned N_THIDS = 512; 
unsigned thid = threadIdx.x, thid_in_warp = thid & WS-1; 
unsigned ofs, cf; 

uint4 a = ((const uint4 *)g_A)[thid], 
     b = ((const uint4 *)g_B)[thid]; 

UADDO(a.x, a.x, b.x) // adding 128-bit chunks with carry flag 
UADDC(a.y, a.y, b.y) 
UADDC(a.z, a.z, b.z) 
UADDC(a.w, a.w, b.w) 
UADDC(cf, 0, 0) // save carry-out 

// memory consumption: 49 * N_THIDS/64 
// use "alternating" data layout for each pair of warps 
volatile short *scan = (volatile short *)(r + 16 + thid_in_warp + 
     49 * (thid/64)) + ((thid/32) & 1); 

scan[-32] = -1; // put identity element 
if(a.x == -1u && a.x == a.y && a.x == a.z && a.x == a.w) 
    // this indicates that carry will propagate through the number 
    cf = -1u; 

// "Hillis-and-Steele-style" reduction 
scan[0] = cf; 
cf = max((int)cf, (int)scan[-2]) & cf; 
scan[0] = cf; 
cf = max((int)cf, (int)scan[-4]) & cf; 
scan[0] = cf; 
cf = max((int)cf, (int)scan[-8]) & cf; 
scan[0] = cf; 
cf = max((int)cf, (int)scan[-16]) & cf; 
scan[0] = cf; 
cf = max((int)cf, (int)scan[-32]) & cf; 
scan[0] = cf; 

int *postscan = (int *)r + 16 + 49 * (N_THIDS/64); 
if(thid_in_warp == WS - 1) // scan leading carry-outs once again 
    postscan[thid >> 5] = cf; 

__syncthreads(); 

if(thid < N_THIDS/32) { 
    volatile int *t = (volatile int *)postscan + thid; 
    t[-8] = -1; // load identity symbol 
    cf = t[0]; 
    cf = max((int)cf, (int)t[-1]) & cf; 
    t[0] = cf; 
    cf = max((int)cf, (int)t[-2]) & cf; 
    t[0] = cf; 
    cf = max((int)cf, (int)t[-4]) & cf; 
    t[0] = cf; 
} 
__syncthreads(); 

cf = scan[0]; 
int ps = postscan[(int)((thid >> 5) - 1)]; // postscan[-1] equals to -1 
scan[0] = max((int)cf, ps) & cf; // update carry flags within warps 
cf = scan[-2]; 

if(thid_in_warp == 0) 
    cf = ps; 
if((int)cf < 0) 
    cf = 0; 

UADDO(a.x, a.x, cf) // propagate carry flag if needed 
UADDC(a.y, a.y, 0) 
UADDC(a.z, a.z, 0) 
UADDC(a.w, a.w, 0) 
((uint4 *)g_R)[thid] = a; 
} 

Zauważ, że makr UADDO/UADDC może nie być konieczne już od CUDA 4.0 posiada odpowiednie intrinsics (jednak nie jestem do końca pewien).

Należy również zauważyć, że choć równoległe zmniejszenie jest dość szybkie, jeśli trzeba dodać kilka dużych liczb całkowitych z rzędu, lepiej byłoby użyć trochę redundantnej reprezentacji (co zasugerowano w komentarzach powyżej), tj. Najpierw zebrać wyniki uzupełnień w słowach 64-bitowych, a następnie wykonać jedną propagację nośną na samym końcu w "jednym przeciągnięciu".

+0

Próbowałem skompilować to, ale mam błąd na ta linia: lotny krótki * skan = (lotny krótki *) (r + 16 + thid_in_warp + (49 * (thid/64)) + ((thid/32) & 1); Wydaje się, że brakuje bliskiego nawiasu. dodając jeden na końcu przed średnikiem.Czy możesz to sprawdzić? (Po ustaleniu, że dostałem awarię uruchamiania próbując użyć go do dodania dwóch 2048x32bit unsigned ints .Może mam błąd w moim własnym kodzie.) –

+0

oops, jesteś w prawo, powinno być (r + 16 + thid_in_warp + 49 * (thid/64)) + ((thid/32) i 1). Naprawiłem to. Proszę również upewnić się, że alokowałeś wystarczającą ilość pamięci dla uruchomienia algorytmu poprawnie, czyli około (49 * 512/64) + 32 o rds. –

+0

Przydzielam 4096 bajtów, co wydaje się być wystarczające. Zaktualizowałem test za pomocą zmodyfikowanej linii. Nadal dostaję nieokreśloną awarię uruchamiania. Być może to jest coś, co robię. –

4

Pomyślałem, że dodam moją odpowiedź, oprócz @asm, więc to pytanie SO może być rodzajem repozytorium pomysłów. Podobnie do @asm, wykrywam i przechowuję warunki przenoszenia, jak również warunek "carry-through", tj. gdy pośrednim wynikiem słowa są wszystkie 1 (0xF ... FFF), więc jeśli przeniesienie zostanie przeniesione do tego słowa, to "przenosi się" do następnego słowa.

Nie użyłem żadnego PTX lub asm w moim kodzie, więc zdecydowałem się użyć 64-bitowych unsigned intów zamiast 32-bitowych, aby osiągnąć zdolność 2048x32bit, używając 1024 wątków.

Większa różnica od kodu @ asm jest w moim schemacie propagacji przenoszenia równoległego. Konstruuję bitową tablicę ("carry"), gdzie każdy bit reprezentuje warunek przenoszenia wygenerowany z niezależnego pośredniego 64-bitowego dodania z każdego z 1024 wątków. Konstruuję również tablicę bitów ("carry_through"), gdzie każdy bit reprezentuje warunek carry_through poszczególnych 64-bitowych wyników pośrednich. W przypadku 1024 wątków odpowiada to 1024/64 = 16x64 bitowym słowom pamięci współużytkowanej dla każdej tablicy bitowej, więc całkowite współużytkowanie pamięci jest 64 + 3 32bitowe. Z tych tablic bitowych zapakowany, to należy wykonać następujące czynności w celu wygenerowania łącznej propagowane wskaźnik posiadające:

carry = carry | (carry_through^((carry & carry_through) + carry_through); 

(uwagę, że przeniesienia jest przesunięta w lewo o jeden: przeprowadzenie [i] wskazuje, że wynik A [I-1] + B [j -1] generowane niosą) wytłumaczeniem jest następujący:

  1. to bitowe i przenoszenia i carry_through generuje kandydatami, gdy transportowa będą interakcję z sekwencją jednego lub więcej przenoszenia pracy w trudnych warunkach
  2. dodanie wyniku z pierwszego kroku do przeniesienia genu stopy którego wynik nie zmienione bity reprezentujące wszystkie słowa, które są dotknięte propagację przeniesienia do carry_through sekwencji
  3. biorąc wyłączne albo z carry_through oraz wyniki z etapu 2 przedstawia wyniki dotkniętych oznaczonych 1-bitowy
  4. przyjmowanie bitowego wyniku lub wyniku z etapu 3 i zwykłych wskaźników przenoszenia daje łączny stan przenoszenia, który jest następnie używany do aktualizacji wszystkich wyników pośrednich.

Należy pamiętać, że dodanie w kroku 2 wymaga kolejnego dodania wielu słów (dla dużych intów złożonych z więcej niż 64 słów). Uważam, że ten algorytm działa i przeszedł pomyślnie przypadki testowe, które zostały przeze mnie rzucone.

jest mój przykład kod, który realizuje to:

// parallel add of large integers 
// requires CC 2.0 or higher 
// compile with: 
// nvcc -O3 -arch=sm_20 -o paradd2 paradd2.cu 
#include <stdio.h> 
#include <stdlib.h> 

#define MAXSIZE 1024 // the number of 64 bit quantities that can be added 
#define LLBITS 64 // the number of bits in a long long 
#define BSIZE ((MAXSIZE + LLBITS -1)/LLBITS) // MAXSIZE when packed into bits 
#define nTPB MAXSIZE 

// define either GPU or GPUCOPY, not both -- for timing 
#define GPU 
//#define GPUCOPY 

#define LOOPCNT 1000 

#define cudaCheckErrors(msg) \ 
    do { \ 
     cudaError_t __err = cudaGetLastError(); \ 
     if (__err != cudaSuccess) { \ 
      fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \ 
       msg, cudaGetErrorString(__err), \ 
       __FILE__, __LINE__); \ 
      fprintf(stderr, "*** FAILED - ABORTING\n"); \ 
      exit(1); \ 
     } \ 
    } while (0) 

// perform c = a + b, for unsigned integers of psize*64 bits. 
// all work done in a single threadblock. 
// multiple threadblocks are handling multiple separate addition problems 
// least significant word is at a[0], etc. 

__global__ void paradd(const unsigned size, const unsigned psize, unsigned long long *c, const unsigned long long *a, const unsigned long long *b){ 

    __shared__ unsigned long long carry_through[BSIZE]; 
    __shared__ unsigned long long carry[BSIZE+1]; 
    __shared__ volatile unsigned mcarry; 
    __shared__ volatile unsigned mcarry_through; 

    unsigned idx = threadIdx.x + (psize * blockIdx.x); 
    if ((threadIdx.x < psize) && (idx < size)){ 
    // handle 64 bit unsigned add first 
    unsigned long long cr1 = a[idx]; 
    unsigned long long lc = cr1 + b[idx]; 
    // handle carry 
    if (threadIdx.x < BSIZE){ 
     carry[threadIdx.x] = 0; 
     carry_through[threadIdx.x] = 0; 
     } 
    if (threadIdx.x == 0){ 
     mcarry = 0; 
     mcarry_through = 0; 
     } 
    __syncthreads(); 
    if (lc < cr1){ 
     if ((threadIdx.x%LLBITS) != (LLBITS-1)) 
     atomicAdd(&(carry[threadIdx.x/LLBITS]), (2ull<<(threadIdx.x%LLBITS))); 
     else atomicAdd(&(carry[(threadIdx.x/LLBITS)+1]), 1); 
     } 
    // handle carry-through 
    if (lc == 0xFFFFFFFFFFFFFFFFull) 
     atomicAdd(&(carry_through[threadIdx.x/LLBITS]), (1ull<<(threadIdx.x%LLBITS))); 
    __syncthreads(); 
    if (threadIdx.x < ((psize + LLBITS-1)/LLBITS)){ 
     // only 1 warp executing within this if statement 
     unsigned long long cr3 = carry_through[threadIdx.x]; 
     cr1 = carry[threadIdx.x] & cr3; 
     // start of sub-add 
     unsigned long long cr2 = cr3 + cr1; 
     if (cr2 < cr1) atomicAdd((unsigned *)&mcarry, (2u<<(threadIdx.x))); 
     if (cr2 == 0xFFFFFFFFFFFFFFFFull) atomicAdd((unsigned *)&mcarry_through, (1u<<threadIdx.x)); 
     if (threadIdx.x == 0) { 
     unsigned cr4 = mcarry & mcarry_through; 
     cr4 += mcarry_through; 
     mcarry |= (mcarry_through^cr4); 
     } 
     if (mcarry & (1u<<threadIdx.x)) cr2++; 
     // end of sub-add 
     carry[threadIdx.x] |= (cr2^cr3); 
     } 
    __syncthreads(); 
    if (carry[threadIdx.x/LLBITS] & (1ull<<(threadIdx.x%LLBITS))) lc++; 
    c[idx] = lc; 
    } 
} 

int main() { 

    unsigned long long *h_a, *h_b, *h_c, *d_a, *d_b, *d_c, *c; 
    unsigned at_once = 256; // valid range = 1 .. 65535 
    unsigned prob_size = MAXSIZE ; // valid range = 1 .. MAXSIZE 
    unsigned dsize = at_once * prob_size; 
    cudaEvent_t t_start_gpu, t_start_cpu, t_end_gpu, t_end_cpu; 
    float et_gpu, et_cpu, tot_gpu, tot_cpu; 
    tot_gpu = 0; 
    tot_cpu = 0; 


    if (sizeof(unsigned long long) != (LLBITS/8)) {printf("Word Size Error\n"); return 1;} 
    if ((c = (unsigned long long *)malloc(dsize * sizeof(unsigned long long))) == 0) {printf("Malloc Fail\n"); return 1;} 

    cudaHostAlloc((void **)&h_a, dsize * sizeof(unsigned long long), cudaHostAllocDefault); 
    cudaCheckErrors("cudaHostAlloc1 fail"); 
    cudaHostAlloc((void **)&h_b, dsize * sizeof(unsigned long long), cudaHostAllocDefault); 
    cudaCheckErrors("cudaHostAlloc2 fail"); 
    cudaHostAlloc((void **)&h_c, dsize * sizeof(unsigned long long), cudaHostAllocDefault); 
    cudaCheckErrors("cudaHostAlloc3 fail"); 

    cudaMalloc((void **)&d_a, dsize * sizeof(unsigned long long)); 
    cudaCheckErrors("cudaMalloc1 fail"); 
    cudaMalloc((void **)&d_b, dsize * sizeof(unsigned long long)); 
    cudaCheckErrors("cudaMalloc2 fail"); 
    cudaMalloc((void **)&d_c, dsize * sizeof(unsigned long long)); 
    cudaCheckErrors("cudaMalloc3 fail"); 
    cudaMemset(d_c, 0, dsize*sizeof(unsigned long long)); 

    cudaEventCreate(&t_start_gpu); 
    cudaEventCreate(&t_end_gpu); 
    cudaEventCreate(&t_start_cpu); 
    cudaEventCreate(&t_end_cpu); 

    for (unsigned loops = 0; loops <LOOPCNT; loops++){ 
    //create some test cases 
    if (loops == 0){ 
    for (int j=0; j<at_once; j++) 
    for (int k=0; k<prob_size; k++){ 
    int i= (j*prob_size) + k; 
    h_a[i] = 0xFFFFFFFFFFFFFFFFull; 
    h_b[i] = 0; 
    } 
    h_a[prob_size-1] = 0; 
    h_b[prob_size-1] = 1; 
    h_b[0] = 1; 
    } 
    else if (loops == 1){ 
    for (int i=0; i<dsize; i++){ 
    h_a[i] = 0xFFFFFFFFFFFFFFFFull; 
    h_b[i] = 0; 
    } 
    h_b[0] = 1; 
    } 
    else if (loops == 2){ 
    for (int i=0; i<dsize; i++){ 
    h_a[i] = 0xFFFFFFFFFFFFFFFEull; 
    h_b[i] = 2; 
    } 
    h_b[0] = 1; 
    } 
    else { 
    for (int i = 0; i<dsize; i++){ 
    h_a[i] = (((unsigned long long)lrand48())<<33) + (unsigned long long)lrand48(); 
    h_b[i] = (((unsigned long long)lrand48())<<33) + (unsigned long long)lrand48(); 
    } 
    } 
#ifdef GPUCOPY 
    cudaEventRecord(t_start_gpu, 0); 
#endif 
    cudaMemcpy(d_a, h_a, dsize*sizeof(unsigned long long), cudaMemcpyHostToDevice); 
    cudaCheckErrors("cudaMemcpy1 fail"); 
    cudaMemcpy(d_b, h_b, dsize*sizeof(unsigned long long), cudaMemcpyHostToDevice); 
    cudaCheckErrors("cudaMemcpy2 fail"); 
#ifdef GPU 
    cudaEventRecord(t_start_gpu, 0); 
#endif 
    paradd<<<at_once, nTPB>>>(dsize, prob_size, d_c, d_a, d_b); 
    cudaCheckErrors("Kernel Fail"); 
#ifdef GPU 
    cudaEventRecord(t_end_gpu, 0); 
#endif 
    cudaMemcpy(h_c, d_c, dsize*sizeof(unsigned long long), cudaMemcpyDeviceToHost); 
    cudaCheckErrors("cudaMemcpy3 fail"); 
#ifdef GPUCOPY 
    cudaEventRecord(t_end_gpu, 0); 
#endif 
    cudaEventSynchronize(t_end_gpu); 
    cudaEventElapsedTime(&et_gpu, t_start_gpu, t_end_gpu); 
    tot_gpu += et_gpu; 
    cudaEventRecord(t_start_cpu, 0); 
    //also compute result on CPU for comparison 
    for (int j=0; j<at_once; j++) { 
    unsigned rc=0; 
    for (int n=0; n<prob_size; n++){ 
     unsigned i = (j*prob_size) + n; 
     c[i] = h_a[i] + h_b[i]; 
     if (c[i] < h_a[i]) { 
     c[i] += rc; 
     rc=1;} 
     else { 
     if ((c[i] += rc) != 0) rc=0; 
     } 
     if (c[i] != h_c[i]) {printf("Results mismatch at offset %d, GPU = 0x%lX, CPU = 0x%lX\n", i, h_c[i], c[i]); return 1;} 
     } 
    } 
    cudaEventRecord(t_end_cpu, 0); 
    cudaEventSynchronize(t_end_cpu); 
    cudaEventElapsedTime(&et_cpu, t_start_cpu, t_end_cpu); 
    tot_cpu += et_cpu; 
    if ((loops%(LOOPCNT/10)) == 0) printf("*\n"); 
    } 
    printf("\nResults Match!\n"); 
    printf("Average GPU time = %fms\n", (tot_gpu/LOOPCNT)); 
    printf("Average CPU time = %fms\n", (tot_cpu/LOOPCNT)); 

    return 0; 
} 
+0

faktycznie wierzę, że moje propagowanie przenoszenia może być dalej zredukowane do: carry = carry | (carry_through^(carry + carry_through)); –

+0

huh ciekawe pomysły)) Sprawdzę twój kod później dzisiaj –

+0

Jest to bardzo przydatne. Czy możesz podać dane dla średniego czasu na procesor vs GPU dla twojego komputera (z określeniem CPu, GPU, OS itd.)? – user3891236