详细了解大量 3x3 矩阵的求逆算法
Understanding in details the algorithm for inversion of a high number of 3x3 matrixes
我按照这个原文 post : 。
建议作为答案的代码是:
$ cat t14.py
import numpy as np
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 >>= 4;
return ret;
}
// in-place is acceptable i.e. out == in)
// T = float or double only
const int block_size = 288;
typedef double T; // *** can set to float or double
__global__ void inv3x3(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;
T det = 1;
if (idx < n*9)
det = in[idx];
unsigned sibase = (threadIdx.x / 9)*9;
unsigned lane = threadIdx.x - sibase; // cheaper modulo
si[threadIdx.x] = det;
__syncthreads();
unsigned off = pat[lane];
T a = si[sibase + getoff(off)];
a *= si[sibase + getoff(off)];
T b = si[sibase + getoff(off)];
b *= si[sibase + getoff(off)];
a -= b;
__syncthreads();
if (lane == 0) si[sibase+3] = a;
if (lane == 3) si[sibase+4] = a;
if (lane == 6) si[sibase+5] = a;
__syncthreads();
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
if (idx < n*9)
out[idx] = a / det;
}
""")
# host code
def gpuinv3x3(inp, n):
# internal constants not to be modified
hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)
# 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*9), dtype= np.float64)
# Get kernel function
matinv3x3 = kernel.get_function("inv3x3")
# Define block, grid and compute
blockDim = (288,1,1) # do not change
gridDim = ((n/32)+1,1,1)
# Kernel function
matinv3x3 (
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, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
n = 2
result = gpuinv3x3(inp, n)
print(result.reshape(2,3,3))
在包含 18 个值(因此 2 个矩阵 3x3)的初始一维数组上,结果给出了正确的倒矩阵,即:
[[[ 2. -0. -1. ]
[-1. -0.33333333 1. ]
[-0. 0.33333333 -0. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]]
主要问题:我想详细了解该算法的工作原理,尤其是当我在大量 3x3矩阵。
我理解这一行:size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
它给出了由当前工作组块的本地 threadIdx 和 blockIdx 标识的当前工作项的全局索引。
我知道 __shared__ T si[block_size];
代表一个共享数组,即关联到工作组块:这就是我们所说的 Local Memory
.
另一方面,我不明白内核代码的以下部分:
__shared__ T si[block_size];
size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
T det = 1;
if (idx < n*9)
det = in[idx];
unsigned sibase = (threadIdx.x / 9)*9;
unsigned lane = threadIdx.x - sibase; // cheaper modulo
si[threadIdx.x] = det;
__syncthreads();
unsigned off = pat[lane];
c
__syncthreads();
if (lane == 0) si[sibase+3] = a;
if (lane == 3) si[sibase+4] = a;
if (lane == 6) si[sibase+5] = a;
__syncthreads();
的确,unsigned sibase = (threadIdx.x / 9)*9;
定义的sibase
索引有什么作用
还有,参数 lane
定义的效用是什么:unsigned lane = threadIdx.x - sibase; // cheaper modulo
最后,移位应用:
T a = si[sibase + getoff(off)];
a *= si[sibase + getoff(off)];
T b = si[sibase + getoff(off)];
b *= si[sibase + getoff(off)];
a -= b;
但是功能没看清楚
关于这部分我也有同样的问题:
if (lane == 0) si[sibase+3] = a;
if (lane == 3) si[sibase+4] = a;
if (lane == 6) si[sibase+5] = a;
行列式的计算方式很奇怪,我无法理解,即:
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
我不是 OpenCL 的初学者,但我还不够专业,无法完全理解这个内核代码。
预赛
首先,了解 3x3 矩阵求逆的算法很重要,请参阅 here(及下文)。
用于内核设计的一般方法是为每个线程分配一个矩阵结果元素。因此每个矩阵需要 9 个线程。最终每个线程将负责计算每个矩阵的 9 个数值结果之一。为了计算两个矩阵,我们需要18个线程,3个矩阵需要27个线程。
一项辅助任务是决定 threadblock/grid 尺码。这遵循典型的方法(总体问题大小决定了所需的线程总数),但我们将具体选择线程块大小 288,因为这是 9(每个矩阵的线程数)和 32(每个矩阵的线程数)的方便倍数CUDA 中每个 warp 的线程数),这为我们提供了一定的效率衡量标准(没有浪费的线程,数据存储中没有间隙)。
由于我们的线程策略是每个矩阵元素一个线程,我们必须使用 9 个线程共同解决矩阵求逆算法。主要任务是计算余因子的转置矩阵,然后计算行列式,然后做最后的算术(除以行列式)计算每个结果元素。
计算余因子
第一个任务是计算 A
的余因子转置矩阵,称为 M
:
|a b c|
let A = |d e f|
|g h i|
|ei-fh ch-bi bf-ce|
M = |fg-di ai-cg cd-af|
|dh-eg bg-ah ae-bd|
这项任务有 9 个线程,矩阵 M
的九个元素要计算,因此我们将为 M
的每个元素分配一个线程。 M
的每个元素都依赖于多个输入值(a
、b
、c
等)所以我们首先加载每个输入值(有 9 个,每个线程), 进入共享内存:
// allocate enough shared memory for one element per thread in the block:
__shared__ T si[block_size];
// compute a globally unique thread index, so each thread has a unique number 0,1,2,etc.
size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
// establish a temporary variable that will use and reuse during thread processing
T det = 1;
// do a thread check to make sure that our next load will be in-bounds for the input array in
if (idx < n*9)
// load one element per thread, 9 threads per matrix will load an entire matrix
det = in[idx];
// for a given matrix (9 threads) compute the base offset into shared memory, where this matrix data (9 elements) will be stored. All 9 threads have the same base offset
unsigned sibase = (threadIdx.x / 9)*9;
// for each group of 9 threads handling a matrix, compute for each thread in that group, a group offset or "lane" from 0..8, so each thread in the group has a unique identifier/assignment in the group
unsigned lane = threadIdx.x - sibase; // cheaper modulo
// let each thread place its matrix element a,b,c, etc. into shared memory
si[threadIdx.x] = det;
// shared memory is now loaded, make sure all threads have loaded before any calculations begin
__syncthreads();
现在每个 A
矩阵元素(a
、b
、c
、...)都已加载到共享内存中,我们可以开始计算M
中的辅因子。让我们关注一个特定线程 (0) 及其辅助因子 (ei-fh
)。计算此余因子所需的所有矩阵元素(e
、i
、f
和 h
)现在都在共享内存中。我们需要一种方法来按顺序加载它们,并执行所需的乘法和减法。
此时我们观察到两件事:
- 每个
M
元素(辅助因子)都有一组不同的 4 个所需元素 A
- 每个
M
元素(辅助因子)遵循相同的通用算法,给定 A
的四个任意元素,让我们将它们统称为 X
、Y
, Z
和 W
。算术是XY-ZW。我取第一个元素,乘以第二个,然后取第三个和第四个元素,将它们相乘,然后减去两个乘积。
由于所有 9 个辅助因子的一般操作顺序(上面的 2)相同,我们只需要一种方法来安排 4 个所需矩阵元素的加载。这种方法被编码到硬编码到示例中的负载模式中:
hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)
有9种负载模式,每种占用一个十六进制数,每个线程一个负载模式,即每个M
矩阵元素(余因子)一个负载模式。在特定的 A
矩阵中,矩阵元素 a
、b
、c
等(已经)加载到 group[=132 的共享内存中=] 0、1、2 等的偏移量。给定线程的加载模式将允许我们生成组偏移量序列,需要从它们在共享内存中的位置检索 A
的矩阵元素,以按顺序用于计算分配给该线程的辅助因子。考虑线程 0 及其辅助因子 ei-fh
,负载模式 0x7584
如何将所需模式编码为 select e
,然后 i
,然后 [=34] =],然后 h
?
为此,我们有一个辅助函数 getoff
,它采用加载模式,并连续(每次调用)剥离索引。我第一次使用参数 0x7584
调用 getoff
时,它“剥离”了索引 4,returns,并将 0x7584
负载模式替换为 0x758
用于下一次使用。 4对应e
。下次我用 0x758
调用 getoff
时,它会“剥离”索引 8,returns,并将 0x758
替换为 0x75
。 8对应i
。下一次产生索引5,对应f
,上一次产生索引7,对应h
.
有了那个描述,然后我们将遍历代码,假装我们是线程 0,并描述计算 ei-fh
:
的过程
// get the load pattern for my matrix "lane"
unsigned off = pat[lane];
//load my temporary variable `a` with the first item indexed in the load pattern:
T a = si[sibase + getoff(off)];
// multiply my temporary variable `a` with the second item indexed in the load pattern
a *= si[sibase + getoff(off)];
//load my temporary variable `b` with the third item indexed in the load pattern
T b = si[sibase + getoff(off)];
// multiply my temporary variable `b` with the fourth item indexed in the load pattern
b *= si[sibase + getoff(off)];
// compute the cofactor by subtracting the 2 products
a -= b;
sibase
,如第一个注释代码部分中所示,是共享内存中存储 A
矩阵元素的基本偏移量。 getoff
函数然后将此基地址添加到 select 相关输入元素。
计算行列式
行列式的数值由下式给出:
det(A) = det = a(ei-fh) - b(di-fg) + c(dh-eg)
如果我们分解它,我们会看到所有的项实际上都已经计算好了:
a,b,c: these are input matrix elements, in shared locations (group offsets) 0, 1, 2
ei-fh: cofactor computed by thread 0
di-fg: cofactor computed by thread 3 (with sign reversed)
dh-eg: cofactor computed by thread 6
现在,每个线程都需要行列式的值,因为每个线程在计算其最终(结果)元素时都会使用它。因此,我们将让矩阵中的每个线程冗余地计算相同的值(这比在一个线程中计算它,然后将该值广播给其他线程更有效)。为了促进这一点,我们需要将 3 个已计算的辅助因子提供给所有 9 个线程。所以我们将 select 共享内存中的 3 个(不再需要)位置来“发布”这些值。我们仍然需要位置 0、1、2 中的值,因为我们需要输入矩阵元素 a
、b
和 c
来计算行列式。但是对于我们剩余的工作,我们不再需要位置 3、4 或 5 中的输入元素,因此我们将重用这些元素:
// we are about to change shared values, so wait until all previous usage is complete
__syncthreads();
// load cofactor computed by thread 0 into group offset 3 in shared
if (lane == 0) si[sibase+3] = a;
// load cofactor computed by thread 3 into group offset 4 in shared
if (lane == 3) si[sibase+4] = a;
// load cofactor computed by thread 6 into group offset 5 in shared
if (lane == 6) si[sibase+5] = a;
// make sure shared memory loads are complete
__syncthreads();
// let every thread compute the determinant (same for all threads)
// a * (ei-fh) + b * -(fg-di) + c * (dh-eg)
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
最终结果的计算
这仅涉及(对于每个线程)将先前为该线程计算的余因子除以刚刚计算的行列式,并存储该结果:
// another thread check: make sure this thread is actually doing useful work
if (idx < n*9)
// take previously computed cofactor, divide by determinant, store result
out[idx] = a / det;
我按照这个原文 post :
$ cat t14.py
import numpy as np
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 >>= 4;
return ret;
}
// in-place is acceptable i.e. out == in)
// T = float or double only
const int block_size = 288;
typedef double T; // *** can set to float or double
__global__ void inv3x3(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;
T det = 1;
if (idx < n*9)
det = in[idx];
unsigned sibase = (threadIdx.x / 9)*9;
unsigned lane = threadIdx.x - sibase; // cheaper modulo
si[threadIdx.x] = det;
__syncthreads();
unsigned off = pat[lane];
T a = si[sibase + getoff(off)];
a *= si[sibase + getoff(off)];
T b = si[sibase + getoff(off)];
b *= si[sibase + getoff(off)];
a -= b;
__syncthreads();
if (lane == 0) si[sibase+3] = a;
if (lane == 3) si[sibase+4] = a;
if (lane == 6) si[sibase+5] = a;
__syncthreads();
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
if (idx < n*9)
out[idx] = a / det;
}
""")
# host code
def gpuinv3x3(inp, n):
# internal constants not to be modified
hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)
# 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*9), dtype= np.float64)
# Get kernel function
matinv3x3 = kernel.get_function("inv3x3")
# Define block, grid and compute
blockDim = (288,1,1) # do not change
gridDim = ((n/32)+1,1,1)
# Kernel function
matinv3x3 (
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, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
n = 2
result = gpuinv3x3(inp, n)
print(result.reshape(2,3,3))
在包含 18 个值(因此 2 个矩阵 3x3)的初始一维数组上,结果给出了正确的倒矩阵,即:
[[[ 2. -0. -1. ]
[-1. -0.33333333 1. ]
[-0. 0.33333333 -0. ]]
[[ 1. 0. 0. ]
[ 0. 1. 0. ]
[ 0. 0. 1. ]]]
主要问题:我想详细了解该算法的工作原理,尤其是当我在大量 3x3矩阵。
我理解这一行:
size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
它给出了由当前工作组块的本地 threadIdx 和 blockIdx 标识的当前工作项的全局索引。我知道
__shared__ T si[block_size];
代表一个共享数组,即关联到工作组块:这就是我们所说的Local Memory
.另一方面,我不明白内核代码的以下部分:
__shared__ T si[block_size]; size_t idx = threadIdx.x+blockDim.x*blockIdx.x; T det = 1; if (idx < n*9) det = in[idx]; unsigned sibase = (threadIdx.x / 9)*9; unsigned lane = threadIdx.x - sibase; // cheaper modulo si[threadIdx.x] = det; __syncthreads(); unsigned off = pat[lane]; c __syncthreads(); if (lane == 0) si[sibase+3] = a; if (lane == 3) si[sibase+4] = a; if (lane == 6) si[sibase+5] = a; __syncthreads();
的确,unsigned sibase = (threadIdx.x / 9)*9;
sibase
索引有什么作用
还有,参数 lane
定义的效用是什么:unsigned lane = threadIdx.x - sibase; // cheaper modulo
最后,移位应用:
T a = si[sibase + getoff(off)];
a *= si[sibase + getoff(off)];
T b = si[sibase + getoff(off)];
b *= si[sibase + getoff(off)];
a -= b;
但是功能没看清楚
关于这部分我也有同样的问题:
if (lane == 0) si[sibase+3] = a; if (lane == 3) si[sibase+4] = a; if (lane == 6) si[sibase+5] = a;
行列式的计算方式很奇怪,我无法理解,即:
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
我不是 OpenCL 的初学者,但我还不够专业,无法完全理解这个内核代码。
预赛
首先,了解 3x3 矩阵求逆的算法很重要,请参阅 here(及下文)。
用于内核设计的一般方法是为每个线程分配一个矩阵结果元素。因此每个矩阵需要 9 个线程。最终每个线程将负责计算每个矩阵的 9 个数值结果之一。为了计算两个矩阵,我们需要18个线程,3个矩阵需要27个线程。
一项辅助任务是决定 threadblock/grid 尺码。这遵循典型的方法(总体问题大小决定了所需的线程总数),但我们将具体选择线程块大小 288,因为这是 9(每个矩阵的线程数)和 32(每个矩阵的线程数)的方便倍数CUDA 中每个 warp 的线程数),这为我们提供了一定的效率衡量标准(没有浪费的线程,数据存储中没有间隙)。
由于我们的线程策略是每个矩阵元素一个线程,我们必须使用 9 个线程共同解决矩阵求逆算法。主要任务是计算余因子的转置矩阵,然后计算行列式,然后做最后的算术(除以行列式)计算每个结果元素。
计算余因子
第一个任务是计算 A
的余因子转置矩阵,称为 M
:
|a b c|
let A = |d e f|
|g h i|
|ei-fh ch-bi bf-ce|
M = |fg-di ai-cg cd-af|
|dh-eg bg-ah ae-bd|
这项任务有 9 个线程,矩阵 M
的九个元素要计算,因此我们将为 M
的每个元素分配一个线程。 M
的每个元素都依赖于多个输入值(a
、b
、c
等)所以我们首先加载每个输入值(有 9 个,每个线程), 进入共享内存:
// allocate enough shared memory for one element per thread in the block:
__shared__ T si[block_size];
// compute a globally unique thread index, so each thread has a unique number 0,1,2,etc.
size_t idx = threadIdx.x+blockDim.x*blockIdx.x;
// establish a temporary variable that will use and reuse during thread processing
T det = 1;
// do a thread check to make sure that our next load will be in-bounds for the input array in
if (idx < n*9)
// load one element per thread, 9 threads per matrix will load an entire matrix
det = in[idx];
// for a given matrix (9 threads) compute the base offset into shared memory, where this matrix data (9 elements) will be stored. All 9 threads have the same base offset
unsigned sibase = (threadIdx.x / 9)*9;
// for each group of 9 threads handling a matrix, compute for each thread in that group, a group offset or "lane" from 0..8, so each thread in the group has a unique identifier/assignment in the group
unsigned lane = threadIdx.x - sibase; // cheaper modulo
// let each thread place its matrix element a,b,c, etc. into shared memory
si[threadIdx.x] = det;
// shared memory is now loaded, make sure all threads have loaded before any calculations begin
__syncthreads();
现在每个 A
矩阵元素(a
、b
、c
、...)都已加载到共享内存中,我们可以开始计算M
中的辅因子。让我们关注一个特定线程 (0) 及其辅助因子 (ei-fh
)。计算此余因子所需的所有矩阵元素(e
、i
、f
和 h
)现在都在共享内存中。我们需要一种方法来按顺序加载它们,并执行所需的乘法和减法。
此时我们观察到两件事:
- 每个
M
元素(辅助因子)都有一组不同的 4 个所需元素A
- 每个
M
元素(辅助因子)遵循相同的通用算法,给定A
的四个任意元素,让我们将它们统称为X
、Y
,Z
和W
。算术是XY-ZW。我取第一个元素,乘以第二个,然后取第三个和第四个元素,将它们相乘,然后减去两个乘积。
由于所有 9 个辅助因子的一般操作顺序(上面的 2)相同,我们只需要一种方法来安排 4 个所需矩阵元素的加载。这种方法被编码到硬编码到示例中的负载模式中:
hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)
有9种负载模式,每种占用一个十六进制数,每个线程一个负载模式,即每个M
矩阵元素(余因子)一个负载模式。在特定的 A
矩阵中,矩阵元素 a
、b
、c
等(已经)加载到 group[=132 的共享内存中=] 0、1、2 等的偏移量。给定线程的加载模式将允许我们生成组偏移量序列,需要从它们在共享内存中的位置检索 A
的矩阵元素,以按顺序用于计算分配给该线程的辅助因子。考虑线程 0 及其辅助因子 ei-fh
,负载模式 0x7584
如何将所需模式编码为 select e
,然后 i
,然后 [=34] =],然后 h
?
为此,我们有一个辅助函数 getoff
,它采用加载模式,并连续(每次调用)剥离索引。我第一次使用参数 0x7584
调用 getoff
时,它“剥离”了索引 4,returns,并将 0x7584
负载模式替换为 0x758
用于下一次使用。 4对应e
。下次我用 0x758
调用 getoff
时,它会“剥离”索引 8,returns,并将 0x758
替换为 0x75
。 8对应i
。下一次产生索引5,对应f
,上一次产生索引7,对应h
.
有了那个描述,然后我们将遍历代码,假装我们是线程 0,并描述计算 ei-fh
:
// get the load pattern for my matrix "lane"
unsigned off = pat[lane];
//load my temporary variable `a` with the first item indexed in the load pattern:
T a = si[sibase + getoff(off)];
// multiply my temporary variable `a` with the second item indexed in the load pattern
a *= si[sibase + getoff(off)];
//load my temporary variable `b` with the third item indexed in the load pattern
T b = si[sibase + getoff(off)];
// multiply my temporary variable `b` with the fourth item indexed in the load pattern
b *= si[sibase + getoff(off)];
// compute the cofactor by subtracting the 2 products
a -= b;
sibase
,如第一个注释代码部分中所示,是共享内存中存储 A
矩阵元素的基本偏移量。 getoff
函数然后将此基地址添加到 select 相关输入元素。
计算行列式
行列式的数值由下式给出:
det(A) = det = a(ei-fh) - b(di-fg) + c(dh-eg)
如果我们分解它,我们会看到所有的项实际上都已经计算好了:
a,b,c: these are input matrix elements, in shared locations (group offsets) 0, 1, 2
ei-fh: cofactor computed by thread 0
di-fg: cofactor computed by thread 3 (with sign reversed)
dh-eg: cofactor computed by thread 6
现在,每个线程都需要行列式的值,因为每个线程在计算其最终(结果)元素时都会使用它。因此,我们将让矩阵中的每个线程冗余地计算相同的值(这比在一个线程中计算它,然后将该值广播给其他线程更有效)。为了促进这一点,我们需要将 3 个已计算的辅助因子提供给所有 9 个线程。所以我们将 select 共享内存中的 3 个(不再需要)位置来“发布”这些值。我们仍然需要位置 0、1、2 中的值,因为我们需要输入矩阵元素 a
、b
和 c
来计算行列式。但是对于我们剩余的工作,我们不再需要位置 3、4 或 5 中的输入元素,因此我们将重用这些元素:
// we are about to change shared values, so wait until all previous usage is complete
__syncthreads();
// load cofactor computed by thread 0 into group offset 3 in shared
if (lane == 0) si[sibase+3] = a;
// load cofactor computed by thread 3 into group offset 4 in shared
if (lane == 3) si[sibase+4] = a;
// load cofactor computed by thread 6 into group offset 5 in shared
if (lane == 6) si[sibase+5] = a;
// make sure shared memory loads are complete
__syncthreads();
// let every thread compute the determinant (same for all threads)
// a * (ei-fh) + b * -(fg-di) + c * (dh-eg)
det = si[sibase]*si[sibase+3]+si[sibase+1]*si[sibase+4]+si[sibase+2]*si[sibase+5];
最终结果的计算
这仅涉及(对于每个线程)将先前为该线程计算的余因子除以刚刚计算的行列式,并存储该结果:
// another thread check: make sure this thread is actually doing useful work
if (idx < n*9)
// take previously computed cofactor, divide by determinant, store result
out[idx] = a / det;