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.

    Related of "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”"

    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.