Cómo realizar la inversión de la matriz de PyCUDA 4×4 con la misma precisión que la función numpy linalg “inv” o “pinv”

Estoy enfrentando un problema de precisión sobre mi código que realiza un número (128, 256, 512) de inversiones de matriz 4×4. Cuando uso la versión original, es decir, la función np.linalg.inv o np.linalg.pinv , todo funciona bien.

Desafortunadamente, con el código CUDA a continuación, obtengo los valores nan e inf en la matriz invertida.

Para ser más explícito, aprovecho esta matriz para invertir:

 2.120771107884677649e+09 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 3.557266600921528288e+27 3.557266600921528041e+07 3.557266600921528320e+17 0.000000000000000000e+00 3.557266600921528041e+07 3.557266600921528288e+27 3.557266600921528041e+07 0.000000000000000000e+00 3.557266600921528320e+17 3.557266600921528041e+07 1.778633300460764144e+27 

Si uso el número clásico ” inv “, obtengo la siguiente matriz invertida de 3×3:

 4.715266047722758306e-10 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 2.811147187396482366e-28 -2.811147186834252285e-48 -5.622294374792964645e-38 0.000000000000000000e+00 -2.811147186834252285e-48 2.811147187396482366e-28 -5.622294374230735768e-48 0.000000000000000000e+00 -5.622294374792964645e-38 -5.622294374230735768e-48 5.622294374792964732e-28 

Para verificar la validez de esta matriz inversa, la he multiplicado por la matriz original y el resultado es la matriz de identidad.

    Pero con la inversión de CUDA GPU, obtengo después de la inversión esta matriz:

     0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 -inf -inf -9.373764907941219970e-01 -inf inf nan -inf nan 

    Por lo tanto, me gustaría boost la precisión en mi kernel CUDA o en el código python para evitar estos valores de nan e inf .

    Aquí está el código del kernel de CUDA y parte de mi código principal (he comentado el método clásico con la función numpy inv :

      # Create arrayFullCross_vec array arrayFullCross_vec = np.zeros((dimBlocks,dimBlocks,integ_prec,integ_prec)) # Create arrayFullCross_vec array invCrossMatrix_gpu = np.zeros((dimBlocks*dimBlocks*integ_prec**2)) # Create arrayFullCross_vec array invCrossMatrix = np.zeros((dimBlocks,dimBlocks,integ_prec,integ_prec)) # Build observables covariance matrix arrayFullCross_vec = buildObsCovarianceMatrix4_vec(k_ref, mu_ref, ir) """ # Compute integrand from covariance matrix for r_p in range(integ_prec): for s_p in range(integ_prec): # original version (without GPU) invCrossMatrix[:,:,r_p,s_p] = np.linalg.inv(arrayFullCross_vec[:,:,r_p,s_p]) """ # GPU version invCrossMatrix_gpu = gpuinv4x4(arrayFullCross_vec.flatten(),integ_prec**2) invCrossMatrix = invCrossMatrix_gpu.reshape(dimBlocks,dimBlocks,integ_prec,integ_prec) """ 

    y aquí el código del kernel CUDA y la función gpuinv4x4 :

     kernel = SourceModule(""" __device__ unsigned getoff(unsigned &off){ unsigned ret = off & 0x0F; off = off >> 4; return ret; } const int block_size = 256; const unsigned tmsk = 0xFFFFFFFF; // in-place is acceptable ie out == in) // T = double or double only typedef double T; __global__ void inv4x4(const T * __restrict__ in, T * __restrict__ out, const size_t n, const unsigned * __restrict__ pat){ __shared__ T si[block_size]; size_t idx = threadIdx.x+blockDim.x*blockIdx.x; if (idx >2)]*a; det += __shfl_down_sync(tmsk, det, 4, 16); // first add det += __shfl_down_sync(tmsk, det, 8, 16); // second add det = __shfl_sync(tmsk, det, 0, 16); // broadcast out[idx] = a / det; } } """) # python function for inverting 4x4 matrices # n should be an even number def gpuinv4x4(inp, n): # internal constants not to be modified hpat = ( 0x0EB51FA5, 0x1EB10FA1, 0x0E711F61, 0x1A710B61, 0x1EB40FA4, 0x0EB01FA0, 0x1E700F60, 0x0A701B60, 0x0DB41F94, 0x1DB00F90, 0x0D701F50, 0x19700B50, 0x1DA40E94, 0x0DA01E90, 0x1D600E50, 0x09601A50, 0x1E790F69, 0x0E391F29, 0x1E350F25, 0x0A351B25, 0x0E781F68, 0x1E380F28, 0x0E341F24, 0x1A340B24, 0x1D780F58, 0x0D381F18, 0x1D340F14, 0x09341B14, 0x0D681E58, 0x1D280E18, 0x0D241E14, 0x19240A14, 0x0A7D1B6D, 0x1A3D0B2D, 0x063D172D, 0x16390729, 0x1A7C0B6C, 0x0A3C1B2C, 0x163C072C, 0x06381728, 0x097C1B5C, 0x193C0B1C, 0x053C171C, 0x15380718, 0x196C0A5C, 0x092C1A1C, 0x152C061C, 0x05281618) # Convert parameters into numpy array # float32 """ inpd = np.array(inp, dtype=np.float32) hpatd = np.array(hpat, dtype=np.uint32) output = np.empty((n*16), dtype= np.float32) """ # float64 """ inpd = np.array(inp, dtype=np.float64) hpatd = np.array(hpat, dtype=np.uint32) output = np.empty((n*16), dtype= np.float64) """ # float128 inpd = np.array(inp, dtype=np.float128) hpatd = np.array(hpat, dtype=np.uint32) output = np.empty((n*16), dtype= np.float128) # Get kernel function matinv4x4 = kernel.get_function("inv4x4") # Define block, grid and compute blockDim = (256,1,1) # do not change gridDim = ((n/16)+1,1,1) # Kernel function matinv4x4 ( cuda.In(inpd), cuda.Out(output), np.uint64(n), cuda.In(hpatd), block=blockDim, grid=gridDim) return output 

    Como puede ver, traté de boost la precisión de la operación de inversión reemplazando np.float32 por np.float64 o np.float128 pero el problema persiste.

    También he sustituido typedef float T; por typedef double T; pero sin éxito.

    ¿Alguien podría ayudarme a realizar la inversión correcta de estas matrices y en su mayoría evitar los valores ‘ nan ‘ e ‘ inf ‘? Creo que este es un verdadero problema de precisión, pero no puedo encontrar una manera de evitar este problema.

    Esta pregunta tiene preguntas relacionadas anteriores aquí y (en menor medida) aquí . No me queda claro por qué el título de la pregunta se refiere a 3×3, el texto en negrita de la pregunta se refiere a 3×3, pero el problema presentado es una matriz inversa de 4×4 (y como se indicó anteriormente a OP, este código solo se puede usar para invertir Matrices 4×4). Continuaré suponiendo que el caso de ejemplo es el caso deseado.

    Según mis pruebas, lo único que es necesario es tomar el código de respuesta anterior y convertirlo para usar double (o, en pycuda, float64 ) en lugar de float (o en pycuda, float32 ). Creo que esto debería ser evidente porque los valores de la matriz de ejemplo exceden el rango de un tipo float32 . Aquí hay un ejemplo trabajado:

     $ cat t10.py import numpy as np # import matplotlib.pyplot as plt import pycuda.driver as cuda from pycuda.compiler import SourceModule import pycuda.autoinit # kernel kernel = SourceModule(""" __device__ unsigned getoff(unsigned &off){ unsigned ret = off & 0x0F; off = off >> 4; return ret; } const int block_size = 256; const unsigned tmsk = 0xFFFFFFFF; // in-place is acceptable ie out == in) // T = float or double only typedef double T; // *** change this typedef to convert between float and double __global__ void inv4x4(const T * __restrict__ in, T * __restrict__ out, const size_t n, const unsigned * __restrict__ pat){ __shared__ T si[block_size]; size_t idx = threadIdx.x+blockDim.x*blockIdx.x; if (idx < n*16){ si[threadIdx.x] = in[idx]; unsigned lane = threadIdx.x & 15; unsigned sibase = threadIdx.x & 0x03F0; __syncwarp(); unsigned off = pat[lane]; T a,b; a = si[sibase + getoff(off)]; a *= si[sibase + getoff(off)]; a *= si[sibase + getoff(off)]; if (!getoff(off)) a = -a; b = si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; if (getoff(off)) a += b; else a -=b; off = pat[lane+16]; b = si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; if (getoff(off)) a += b; else a -=b; b = si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; if (getoff(off)) a += b; else a -=b; off = pat[lane+32]; b = si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; if (getoff(off)) a += b; else a -=b; b = si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; b *= si[sibase + getoff(off)]; if (getoff(off)) a += b; else a -=b; T det = si[sibase + (lane>>2)]*a; det += __shfl_down_sync(tmsk, det, 4, 16); // first add det += __shfl_down_sync(tmsk, det, 8, 16); // second add det = __shfl_sync(tmsk, det, 0, 16); // broadcast out[idx] = a / det; } } """) # host code def gpuinv4x4(inp, n): # internal constants not to be modified hpat = ( 0x0EB51FA5, 0x1EB10FA1, 0x0E711F61, 0x1A710B61, 0x1EB40FA4, 0x0EB01FA0, 0x1E700F60, 0x0A701B60, 0x0DB41F94, 0x1DB00F90, 0x0D701F50, 0x19700B50, 0x1DA40E94, 0x0DA01E90, 0x1D600E50, 0x09601A50, 0x1E790F69, 0x0E391F29, 0x1E350F25, 0x0A351B25, 0x0E781F68, 0x1E380F28, 0x0E341F24, 0x1A340B24, 0x1D780F58, 0x0D381F18, 0x1D340F14, 0x09341B14, 0x0D681E58, 0x1D280E18, 0x0D241E14, 0x19240A14, 0x0A7D1B6D, 0x1A3D0B2D, 0x063D172D, 0x16390729, 0x1A7C0B6C, 0x0A3C1B2C, 0x163C072C, 0x06381728, 0x097C1B5C, 0x193C0B1C, 0x053C171C, 0x15380718, 0x196C0A5C, 0x092C1A1C, 0x152C061C, 0x05281618) # Convert parameters into numpy array # *** change next line between float32 and float64 to match float or double inpd = np.array(inp, dtype=np.float64) hpatd = np.array(hpat, dtype=np.uint32) # *** change next line between float32 and float64 to match float or double output = np.empty((n*16), dtype= np.float64) # Get kernel function matinv4x4 = kernel.get_function("inv4x4") # Define block, grid and compute blockDim = (256,1,1) # do not change gridDim = ((n/16)+1,1,1) # Kernel function matinv4x4 ( cuda.In(inpd), cuda.Out(output), np.uint64(n), cuda.In(hpatd), block=blockDim, grid=gridDim) return output inp = (1.0, 1.0, 1.0, 0.0, 0.0, 3.0, 1.0, 2.0, 2.0, 3.0, 1.0, 0.0, 1.0, 0.0, 2.0, 1.0, 2.120771107884677649e+09, 0.0, 0.0, 0.0, 0.0, 3.557266600921528288e+27, 3.557266600921528041e+07, 3.557266600921528320e+17, 0.0, 3.557266600921528041e+07, 3.557266600921528288e+27, 3.557266600921528041e+07, 0.0, 3.557266600921528320e+17, 3.557266600921528041e+07, 1.778633300460764144e+27) n = 2 result = gpuinv4x4(inp, n) print(result) $ python t10.py [ -3.00000000e+00 -5.00000000e-01 1.50000000e+00 1.00000000e+00 1.00000000e+00 2.50000000e-01 -2.50000000e-01 -5.00000000e-01 3.00000000e+00 2.50000000e-01 -1.25000000e+00 -5.00000000e-01 -3.00000000e+00 -0.00000000e+00 1.00000000e+00 1.00000000e+00 4.71526605e-10 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 2.81114719e-28 -2.81114719e-48 -5.62229437e-38 0.00000000e+00 -2.81114719e-48 2.81114719e-28 -5.62229437e-48 0.00000000e+00 -5.62229437e-38 -5.62229437e-48 5.62229437e-28] $ 

    Además de actualizar la matriz de entrada para utilizar la que se proporciona en esta pregunta, tenga en cuenta que solo se cambiaron 3 líneas de código de la respuesta anterior para convertir de punto flotante de 32 bits a punto flotante de 64 bits. Cada una de esas 3 líneas está marcada por un comentario que contiene triple asterisco (***).

    Además, si le preocupa la precisión más allá de los primeros 9 o más dígitos que se muestran aquí, no he investigado eso. Puede ser que este código no sea numéricamente adecuado para cada caso o necesidad.

    Además, el código de la pregunta también muestra un tipo float128 en la sección de python. No hay ningún tipo de CUDA nativo que corresponda a eso.