Как выполнить инверсию матрицы PyCUDA 4x4 с той же точностью, что и для функции «inv» или «pinv» ninaly linalg - PullRequest
Новая
       112

Как выполнить инверсию матрицы PyCUDA 4x4 с той же точностью, что и для функции «inv» или «pinv» ninaly linalg

0 голосов
/ 27 марта 2019

Я столкнулся с проблемой точности моего кода, который выполняет ряд (128, 256, 512) инверсий матрицы 4x4.Когда я использую оригинальную версию, то есть функцию numpy np.linalg.inv или np.linalg.pinv, все работает отлично.

К сожалению, с кодом CUDA ниже, я получаю значения nan и inf в инвертированной матрице.

Чтобы быть более точным, я возьму эту матрицу для инвертирования: 

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

Если я использую классический numpy "inv", я получу следующую инвертированную матрицу 3x3: 

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

Чтобы проверить достоверность этой обратной матрицы, я умножил ее на исходную матрицу, и в результате получилась единичная матрица.

Но с инверсией графического процессора CUDA я получаю после инверсии эту матрицу: 

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

Итак, я хотел бы повысить точность в моем коде ядра или Python CUDA, чтобы избежать этих значений nan и inf.

Вот код ядра CUDAи вызывая часть моего основного кода (я прокомментировал классический метод с помощью функции 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)
    """

, а здесь код ядра CUDA и функция 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 i.e. 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 < 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;
  }
}
""")

# 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

Как видите, я попытался увеличитьпростота операции инвертирования путем замены np.float32 на np.float64 или np.float128, но проблема остается.

Я также заменил typedef float T; на typedef double T;, но безуспешно.

Любой могпомогите мне выполнить правильную инверсию этих матриц и в основном избежать значений 'nan' и 'inf'?Я думаю, что это реальная проблема точности, но я не могу найти, как обойти эту проблему.

1 Ответ

5 голосов
/ 28 марта 2019

У этого вопроса есть предыдущие вопросы здесь и (в меньшей степени) здесь .Мне непонятно, почему название вопроса относится к 3x3, жирный текст в вопросе относится к 3x3, но представленная проблема является обратной матрицей 4x4 (и, как было указано ранее в OP, этот код может только используется для инвертирования матриц 4х4).Я буду исходить из предположения, что пример является желаемым.

В соответствии с моим тестированием, единственное, что необходимо, это взять предыдущий код ответа и преобразовать его в использование.double (или у pycuda float64) вместо float (или у pycuda float32).Я думаю, что это должно быть очевидно, потому что значения матрицы примера превышают диапазон типа float32.Вот рабочий пример: 

$ 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 i.e. 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]
$

Помимо обновления матрицы ввода для использования матрицы, представленной в этом вопросе, обратите внимание, что только 3 строки кода были изменены по сравнению с предыдущим ответом для преобразования из 32-битного плавающегоуказывают на 64-битную плавающую точку.Каждая из этих 3 строк помечена комментарием, содержащим тройную звездочку (***).

Кроме того, если вас беспокоит точность, превышающая первые 9 или около того цифр, отображаемых здесь, я не имею 'Я исследовал это.Может случиться так, что этот код численно не подходит для каждого случая или потребности.

Кроме того, код в вопросе также показывает тип float128 в разделе Python.Не существует собственного типа CUDA, который соответствует этому.

...