2012-10-18 4 views
11

Ho sviluppato un algoritmo di crittografia sulla GPU e attualmente sono bloccato con un algoritmo per eseguire l'aggiunta di interi di grandi dimensioni. I grandi numeri interi sono rappresentati in modo normale come un gruppo di parole a 32 bit.aggiunta intera di grandi dimensioni con CUDA

Ad esempio, è possibile utilizzare un thread per aggiungere due parole a 32 bit. Per semplicità, supponiamo che i numeri da aggiungere siano della stessa lunghezza e numero di thread per blocco == numero di parole. Poi:

__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; 
} 

Sono abbastanza sicuro che ci sia un modo per farlo trasportare propagazione attraverso qualche procedura di riduzione difficile, ma non riusciva a capirlo ..

ho dato un'occhiata a CUDA thrust extensions ma grande pacchetto intero sembra non essere ancora implementato. Forse qualcuno può darmi un suggerimento su come farlo su CUDA?

+2

La GPU in grado di gestire fino a 64 bit (long long) direttamente. Un approccio per 128 bit è delineato in [questa domanda/risposta SO] (http://stackoverflow.com/questions/6162140/128-bit-integer-on-cuda). –

+0

Penso che quello che vuoi da CUDA possa essere ottenuto dalle tecniche C. Pertanto, ho ritirato la domanda anche in 'C'. Spero di ottenere una buona risposta dagli esperti di C. – ahmad

+0

Sì, puoi anche programmare un'aggiunta intera lunga usando solo costrutti C di alto livello (in contrapposizione all'assembly linea PXT in CUDA), ma richiederebbe molte più istruzioni, come ho sottolineato in questa risposta: http: // stackoverflow .com/questions/12448549/is-inline-ptx-assembly-code-powerful/12453534 # 12453534 – njuffa

risposta

8

Hai ragione, la propagazione può essere effettuata tramite calcolo della somma del prefisso, ma è un po 'complicato definire la funzione binaria per questa operazione e dimostrare che è associativa (necessaria per la somma del prefisso parallelo). In realtà, questo algoritmo viene utilizzato (in teoria) in Carry-lookahead adder.

Supponiamo di avere due interi grandi a [0..n-1] e b [0..n-1]. Poi calcoliamo (i = 0..n-1):

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

Definiamo due funzioni:

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

con significato molto intuitivo: generiamo [i] == 1 significa che il riporto viene generato in posizione mentre propagando [i] == 1 significa che il trasporto verrà propagato dalla posizione (i - 1) a (i + 1). Il nostro obiettivo è calcolare il carryout della funzione [0..n-1] usato per aggiornare la somma risultante s [0..n-1]. carryout può essere calcolata in modo ricorsivo come segue:

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

Qui asporto [i] == 1 se riporto è generato alla posizione i O viene generato volte prima e propagato per posizionare i. Infine, si aggiorna la somma risultante:

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

Ora è abbastanza semplice per dimostrare che la funzione asporto è effettivamente associativo binario e quindi parallelo computazione prefisso somma applica. Per implementare questo su CUDA, possiamo unire 'generare' entrambe le bandierine e 'propagare' in una singola variabile dal momento che sono mutuamente esclusivi, ossia:

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

In altre parole,

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

Poi, si può verificare che la seguente formula calcola somma prefisso per funzione asporto:

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

per ogni k < i. Il codice di esempio seguente mostra una grande aggiunta per interi di 2048 parole. Qui ho usato blocchi CUDA con 512 discussioni:

// 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; 
} 

Nota che le macro UADDO/UADDC potrebbe non essere più necessario dal momento CUDA 4.0 ha intrinseche corrispondenti (ma non sono del tutto sicuro).

Inoltre, sebbene la riduzione parallela sia abbastanza veloce, se è necessario aggiungere diversi numeri interi di grandi dimensioni in una riga, potrebbe essere preferibile utilizzare una rappresentazione ridondante (che è stata suggerita nei commenti sopra), vale a dire, prima accumulare il risultati di aggiunte in parole a 64 bit e quindi eseguire una propagazione di carry alla fine in "one sweep".

+0

Ho provato a compilarlo, ma ho ricevuto un errore su questa linea: volatile short * scan = (volatile breve *) (r + 16 + thid_in_warp + (49 * (th/64)) + ((thid/32) & 1); Sembra mancare una parentesi chiusa. aggiungendo uno alla fine prima del punto e virgola.Potresti controllare? (Dopo aver risolto che ho avuto un errore di avvio cercando di usarlo per aggiungere due pollici non firmati 2048x32 bit. Potrei avere un errore nel mio codice.) –

+0

oops, sei a destra, dovrebbe essere (r + 16 + thid_in_warp + 49 * (th/64)) + ((thid/32) & 1). L'ho corretto. Assicurati inoltre di allocare abbastanza mem di parti per l'algoritmo da eseguire correttamente che è circa (49 * 512/64) + 32 wo RDS. –

+0

Sto allocando 4096 byte che sembra dovrebbe essere più che sufficiente. Ho aggiornato il mio test con la tua linea modificata. Sto ancora ricevendo l'errore di avvio non specificato. Forse è qualcosa che sto facendo. –

4

Ho pensato di pubblicare anche la mia risposta, oltre a @asm, quindi questa domanda SO può essere una sorta di repository di idee. Analogamente a @asm, rilevo e memorizzo la condizione carry e anche la condizione "carry-through", ovvero. quando il risultato della parola intermedia è tutto 1 (0xF ... FFF) in modo che se un carry si propaga in questa parola, esso "trasporta" alla parola successiva.

Non ho utilizzato alcun PTX o asm nel mio codice, quindi ho scelto di utilizzare 64 bit senza segno anziché 32 bit, per ottenere la capacità di 2048 x 32 bit, utilizzando 1024 thread.

Una differenza maggiore rispetto al codice di @ asm si trova nel mio schema di trasmissione di trasporto parallelo. Costruisco un array bit-packed ("carry") in cui ogni bit rappresenta la condizione di carry generata dagli addizionali intermedi indipendenti a 64 bit da ciascuno dei 1024 thread. Costruisco anche un array bit-packed ("carry_through") in cui ogni bit rappresenta la condizione carry_through dei singoli risultati intermedi a 64 bit. Per 1024 thread, ciò equivale a 1024/64 = 16x64 bit di memoria condivisa per ogni array con bit, quindi l'utilizzo di mem condiviso totale è 64 + 3 quantiti a 32 bit. Con queste matrici bit imballato, effettuo quanto segue per generare un combinato propagato indicatore riporto:

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

(si noti che riporto viene spostato a sinistra di uno: portare [i] indica che il risultato di un [i-1] + b [i-1] generato un riporto) la spiegazione è la seguente:

  1. il bit e di riporto e carry_through genera i candidati cui un carry interagire con una sequenza di uno o più di riporto se le condizioni
  2. aggiungendo il risultato del passaggio uno al gene carry_through tassi un risultato che è cambiato bit che rappresentano tutte le parole che saranno interessate dal propagazione del riporto nella sequenza carry_through
  3. prendendo la OR esclusivo del carry_through più il risultato del punto 2 mostra i risultati affetti indicati con a 1 bit
  4. prendendo il bit a bit o il risultato dal punto 3 e gli indicatori di trasporto ordinari danno una condizione di trasporto combinata, che è quindi utilizzata per aggiornare tutti i risultati intermedi.

Si noti che l'aggiunta nel passaggio 2 richiede un'altra aggiunta di più parole (per i big ints composti da più di 64 parole). Credo che questo algoritmo funzioni e ha superato i test che ho lanciato.

Ecco il mio codice di esempio che implementa questo:

// 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

in realtà credo che la mia propagazione del carry possa essere ulteriormente ridotta a: carry = carry | (carry_through^(carry + carry_through)); –

+0

eh idee interessanti)) Controllerò il tuo codice più tardi oggi –

+0

Questo è molto utile. Puoi dare le cifre per tempo medio su CPU vs GPU per la tua macchina (affermando CPu, GPU, OS ecc.)? – user3891236