cuda - residuales - historia de la aritmetica modular
aritmética modular en la GPU (3)
Hace algún tiempo experimenté 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, la 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 solo funciona para los enteros de 31 bits (si 1 bit no es crítico para usted) y también necesita precomputar ''invk'' de antemano. Esto da el 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, puedes 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, comparé el rendimiento del algoritmo resultante usando su multiplicación modular:
unsigned r = (unsigned)(((u64)a * (u64)b) % m);
contra la versión optimizada con mul_m.
Aritmética modular con% de operación 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;
Entonces, 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 el rendimiento aritmético en bruto utilizando un kernel con un montón de operaciones mul_mod (como el ejemplo de saxpy ). Pero en aplicaciones reales con lógica de control, barreras de sincronización, etc., la aceleración es muy notable.
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 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 de 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 esto se implementa en CUDA, utilizo 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;
}
Se supone que este código no funciona, solo quería ver cómo se implementa la reducción modular en CUDA.
Cuando desarmo esto con cuobjdump --dump-sass (gracias njuffa 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 llamada al procedimiento 0xf8 que implementa algún algoritmo sofisticado (unas 50 instrucciones o incluso más). No sorprende que la operación mod (%) sea lenta
Una GPU Fermi de alta gama (por ejemplo, una GTX 580) probablemente le brinde 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 entre tres y cinco instrucciones de la 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.
[Agregado después, después de la aclaración de la pregunta:]
La división de 64 bits sin signo y el módulo con divisor variable por otro lado se llaman subrutinas de aproximadamente 65 instrucciones en 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;
Hay trucos para realizar eficientemente operaciones mod, pero si solo 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 de la implementación eficiente del módulo. Es posible que deba implementarlo usted mismo para aprovecharlo al máximo.