aritmética modular en la gpu

Question

Estoy trabajando en el algoritmo de GPU que se supone que debe hacer muchos cálculos modulares. Particularmente, varias operaciones en matrices en un campo finito que a la larga se reducen a operaciones primitivas como: (a * b - c * d) mod m o (a * b + c) mod m donde a, b, c y d son residuos módulo m y m es un primo de 32 bits.

A través de la experimentación, aprendí que el rendimiento del algoritmo está limitado en gran medida por una aritmética modular lenta porque las operaciones de módulo entero (%) y división no son compatibles con la GPU en el hardware.

Agradezco si alguien me puede dar una idea de cómo realizar cálculos modulares eficientes con CUDA?

Para ver cómo se implementa en CUDA, yo uso el siguiente fragmento de código:

__global__ void mod_kernel(unsigned *gout, const unsigned *gin) { 

unsigned tid = threadIdx.x; 
unsigned a = gin[tid], b = gin[tid * 2], m = gin[tid * 3]; 

typedef unsigned long long u64; 

__syncthreads(); 
unsigned r = (unsigned)(((u64)a * (u64)b) % m); 
__syncthreads(); 
gout[tid] = r; 
} 

no se supone que este código funcione, sólo quería ver cómo la reducción modular está implementado en CUDA.

Cuando desmonte este con cuobjdump --dump-Sass (! Njuffa gracias por el consejo), veo lo siguiente:

/*0098*/  /*0xffffdc0450ee0000*/  BAR.RED.POPC RZ, RZ; 
/*00a0*/  /*0x1c315c4350000000*/  IMUL.U32.U32.HI R5, R3, R7; 
/*00a8*/  /*0x1c311c0350000000*/  IMUL.U32.U32 R4, R3, R7; 
/*00b0*/  /*0xfc01dde428000000*/  MOV R7, RZ; 
/*00b8*/  /*0xe001000750000000*/  CAL 0xf8; 
/*00c0*/  /*0x00000007d0000000*/  BPT.DRAIN 0x0; 
/*00c8*/  /*0xffffdc0450ee0000*/  BAR.RED.POPC RZ, RZ; 

Tenga en cuenta que entre las dos llamadas a bar.red.popc hay una llame al procedimiento 0xf8 que implementa algún algoritmo sofisticado (alrededor de 50 instrucciones o incluso más). No sorprende que la operación de mod (%) sea lenta

Answer 1

Hace algún tiempo he experimentado mucho con la aritmética modular en la GPU. En las GPU de Fermi, puede utilizar la aritmética de precisión doble para evitar costosas operaciones de div y mod. Por ejemplo, multiplicación modular se puede hacer de la siguiente manera:

// fast truncation of double-precision to integers 
#define CUMP_D2I_TRUNC (double)(3ll << 51) 
// computes r = a + b subop c unsigned using extended precision 
#define VADDx(r, a, b, c, subop) \ 
    asm volatile("vadd.u32.u32.u32." subop " %0, %1, %2, %3;" : \ 
      "=r"(r) : "r"(a) , "r"(b), "r"(c)); 

// computes a * b mod m; invk = (double)(1<<30)/m 
__device__ __forceinline__ 
unsigned mul_m(unsigned a, unsigned b, volatile unsigned m, 
    volatile double invk) { 

    unsigned hi = __umulhi(a*2, b*2); // 3 flops 
    // 2 double instructions 
    double rf = __uint2double_rn(hi) * invk + CUMP_D2I_TRUNC; 
    unsigned r = (unsigned)__double2loint(rf); 
    r = a * b - r * m; // 2 flops 

    // can also be replaced by: VADDx(r, r, m, r, "min") // == umin(r, r + m); 
    if((int)r < 0) 
     r += m; 
    return r; 
} 

Sin embargo esto sólo funciona para modulos enteros de 31 bits (si 1 bit no es crítica para usted) y también tiene que calcular previamente 'invk' de antemano. Esto da mínimo absoluto de instrucciones que puedo lograr, es decir .:

SHL.W R2, R4, 0x1; 
SHL.W R8, R6, 0x1; 
IMUL.U32.U32 R4, R4, R6; 
IMUL.U32.U32.HI R8, R2, R8; 
I2F.F64.U32 R8, R8; 
DFMA R2, R2, R8, R10; 
IMAD.U32.U32 R4, -R12, R2, R4; 
ISETP.GE.AND P0, pt, R4, RZ, pt; 
@!P0 IADD R4, R12, R4; 

Para la descripción del algoritmo, se puede echar un vistazo a mi artículo: gpu_resultants. Otras operaciones como (x y - z w) mod m también se explican allí.

Por curiosidad, he comparado el rendimiento del algoritmo resultante la utilización de su multiplicación modular:

unsigned r = (unsigned)(((u64)a * (u64)b) % m); 

contra la versión optimizada con mul_m.

aritmética modular con funcionamiento por defecto%:

low_deg: 11; high_deg: 2481; bits: 10227 
nmods: 330; n_real_pts: 2482; npts: 2495 

res time: 5755.357910 ms; mod_inv time: 0.907008 ms; interp time: 856.015015 ms; CRA time: 44.065857 ms 
GPU time elapsed: 6659.405273 ms; 

aritmética modular con mul_m:

low_deg: 11; high_deg: 2481; bits: 10227 
nmods: 330; n_real_pts: 2482; npts: 2495 

res time: 1100.124756 ms; mod_inv time: 0.192608 ms; interp time: 220.615143 ms; CRA time: 10.376352 ms 
GPU time elapsed: 1334.742310 ms; 

Así que en promedio es aproximadamente 5 veces más rápido. Tenga en cuenta también que, es posible que no vea una aceleración si solo evalúa rendimiento aritmético bruto usando un kernel con un montón de operaciones mul_mod (como saxpy ejemplo). Pero en aplicaciones reales con lógica de control, barreras de sincronización, etc., la aceleración es muy notable.

Answer 2

Una GPU Fermi de alta gama (por ejemplo, una GTX 580) probablemente le proporcione el mejor rendimiento entre las tarjetas de envío para esto. Querría que todos los operandos de 32 bits fueran del tipo "unsigned int" para obtener el mejor rendimiento, ya que hay una sobrecarga adicional para el manejo de divisiones y módulos firmados.

El compilador genera un código muy eficiente para la división y el módulo con un divisor fijo. Como recuerdo, por lo general hay alrededor de tres a cinco instrucciones de instrucciones de máquina en Fermi y Kepler. Puede verificar el SASS generado (código de máquina) con cuobjdump --dump-sass. Es posible que pueda utilizar funciones con plantilla con divisores constantes si solo usa algunos divisores diferentes.

Debería ver en el orden de dieciséis instrucciones SASS en línea que se generan para las operaciones sin signo de 32 bits con divisor variable, a través de Fermi y Kepler. El código está limitado por el rendimiento de las multiplicaciones enteras y para las GPU de clase Fermi es competitivo con las soluciones de hardware. Se observa un rendimiento algo reducido en las GPU de clase Kepler actualmente enviadas debido a su rendimiento de multiplicación de enteros reducido.

[añadió más tarde, después de la clarificación de la cuestión:]

sin signo de división de 64 bits y modulo con divisor variable sobre la otra parte se llaman subrutinas de unos 65 instrucciones sobre Fermi y Kepler. Se ven cercanos a lo óptimo. En Fermi, esto todavía es razonablemente competitivo con las implementaciones de hardware (tenga en cuenta que las divisiones enteras de 64 bits no son exactamente súper rápidas en las CPU que proporcionan esto como una instrucción incorporada). A continuación se muestra un código que publiqué en los foros de NVIDIA hace algún tiempo para el tipo de tarea descrita en la aclaración. Evita la costosa división, pero supone que lotes bastante grandes de operandos comparten la misma división. Utiliza una aritmética de doble precisión, que es especialmente rápida en las GPU de clase Tesla (a diferencia de las tarjetas de consumidor). Solo hice una prueba superficial del código, es posible que desee probar esto con más cuidado antes de implementarlo.

// Let b, p, and A[i] be integers < 2^51 
// Let N be a integer on the order of 10000 
// for i from 1 to N 
// A[i] <-- A[i] * b mod p 

/*---- kernel arguments ----*/ 
unsigned long long *A; 
double b, p; /* convert from unsigned long long to double before passing to kernel */ 
double oop; /* pass precomputed 1.0/p to kernel */ 

/*---- code inside kernel -----*/ 
double a, q, h, l, rem; 
const double int_cvt_magic = 6755399441055744.0; /* 2^52+2^51 */ 

a = (double)A[i]; 

/* approximate quotient and round it to the nearest integer */ 
q = __fma_rn (a * b, oop, int_cvt_magic); 
q = q - int_cvt_magic; 

/* back-multiply, representing p*q as a double-double h:l exactly */ 
h = p * q; 
l = __fma_rn (p, q, -h); 

/* remainder is double-width product a*b minus double-double h:l */ 
rem = __fma_rn (a, b, -h); 
rem = rem - l; 

/* remainder may be negative as quotient rounded; fix if necessary */ 
if (rem < 0.0) rem += p; 

A[i] = (unsigned long long)rem; 

Answer 3

Hay trucos para llevar a cabo eficientemente las operaciones mod, pero sólo si m es radix 2.

Por ejemplo, x mod y == x & (y-1), donde y es 2^n. La operación bit a bit es la más rápida.

De lo contrario, probablemente una tabla de búsqueda? A continuación se muestra un enlace sobre la discusión sobre la implementación eficiente del módulo. Es posible que deba implementarlo usted mismo para aprovecharlo al máximo.

Efficient computation of mod