如何以与numpy linalg"inv"相同的精度执行PyCUDA 4x4矩阵求逆.或"pinv"功能 [英] How to perform PyCUDA 4x4 matrix inversion with same accuracy than numpy linalg "inv" or "pinv" function
问题描述
我遇到了我的代码的准确性问题，该代码执行4x4矩阵求逆(128、256、512)个数.当我使用原始版本(即numpy函数np.linalg.inv或np.linalg.pinv)时，一切正常.
I am facing an issue of accuracy about my code which performs a number (128, 256, 512) of 4x4 matrix inversions. When I use the original version, i.e the numpy function np.linalg.inv or np.linalg.pinv, everything works fine.
不幸的是，使用下面的CUDA代码，我将nan和inf的值转换为倒置矩阵.
Unfortunately, with the CUDA code below, I get nan and inf values into inverted matrix.
更明确地说，我将这个矩阵求反:
To be more explicit, I take this matrix to invert :
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倒置矩阵:
If I use classical numpy "inv", I get for the following inverted 3x3 matrix :
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
要检查此逆矩阵的有效性，我将其乘以原始矩阵，得出的结果是单位矩阵.
To check the validity of this inverse matrix, I have multiplied it by original matrix and the result is the identity matrix.
但是使用CUDA GPU反转后，我得到了这个矩阵:
But with CUDA GPU inversion, I get after the inversion this matrix :
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
因此，我想提高CUDA内核或python代码的精度，以避免出现这些nan和inf值.
So, I woul like to increase the precision into my CUDA kernel or python code to avoid these nanand inf values.
这是CUDA内核代码，并调用了我的主要代码的一部分(我已经用numpy inv函数注释了经典方法:
Here is the CUDA kernel code and calling part of my main code (I have commented the classical method with numpy inv function :
# 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函数:
and here the CUDA kernel code and gpuinv4x4 function :
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来提高反相操作的精度，但是问题仍然存在.
As you can see, I tried to increase accuracy of inverting operation by replacing np.float32 by np.float64 or np.float128 but problem remains.
我也用typedef double T;替换了typedef float T;，但没有成功.
I have also replaced typedef float T; by typedef double T;but without success.
任何人都可以帮助我对这些矩阵进行正确的求逆，并且大多避免使用'nan'和'inf'值吗?我认为这是一个真正的精度问题，但我找不到解决该问题的方法.
Anyone could help me to perform the right inversion of these matrices and mostly avoid the 'nan' and 'inf' values ? I think this is a real issue of precision but I can't find how to circumvent this problem.
推荐答案
此问题以前有相关问题此处.我不清楚问题的标题为何指向3x3，问题中的粗体文本指向3x3，但是出现的问题是4x4矩阵逆的(如先前OP所述，此代码只能 用于反转4x4矩阵).我将假设示例案例是理想案例.
This question has previous related questions here and (to a lesser extent) here. It's unclear to me why the title of the question refers to 3x3, the bolded text in the question refers to to 3x3, but the presented problem is a 4x4 matrix inverse (and as stated to OP previously, this code can only be used to invert 4x4 matrices). I will proceed under the assumption that the example case is the desired case.
根据我的测试，唯一需要做的就是获取先前的答案代码并将其转换为double(或者在pycuda中，是float64)而不是float(或者在pycuda中，float32).我认为这应该是显而易见的，因为示例矩阵值超出了float32类型的范围.这是一个可行的示例:
According to my testing, the only thing that is necessary is to take the previous answer code and convert it to use double (or, in pycuda, float64) rather than float (or in pycuda, float32). I think this should be evident because the example matrix values exceed the range of a float32 type. Here is a worked example:
$ 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行中的每行都用一个注释标记，该注释中包含三星号(***).
Apart from updating the input matrix to use the one supplied in this question, note that only 3 lines of code were changed from the previous answer to convert from 32-bit floating point to 64-bit floating point. Each of those 3 lines are marked by a comment that contains triple-asterisk (***) in it.
此外，如果您担心此处显示的前9个数字以外的准确性，我没有进行调查.可能该代码在数字上并不适合所有情况或需求.
Also, if you are concerned about accuracy beyond the first 9 or so digits displayed here, I haven't investigated that. It may be that this code is not numerically suitable for every case or need.
顺便说一句，问题中的代码在python部分中还显示了float128类型.没有与之对应的本机CUDA类型.
As an aside, the code in the question also shows a float128 type in the python section. There is no native CUDA type that corresponds to that.
这篇关于如何以与numpy linalg"inv"相同的精度执行PyCUDA 4x4矩阵求逆.或"pinv"功能的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
