预赛
首先，了解 3x3 矩阵求逆的算术很重要，请参阅here https://ardoris.wordpress.com/2008/07/18/general-formula-for-the-inverse-of-a-3x3-matrix/（及下文）。
用于内核设计的一般方法是为每个线程分配一个矩阵结果元素。因此，每个矩阵需要 9 个线程。最终，每个线程将负责计算每个矩阵的 9 个数值结果之一。为了计算两个矩阵，我们需要 18 个线程，3 个矩阵需要 27 个线程。
辅助任务是决定线程块/网格大小。这遵循典型的方法（总体问题大小决定所需的线程总数），但我们将为线程块大小选择 288，因为这是 9（每个矩阵的线程数）和 32（每个矩阵的线程数）的方便倍数。 CUDA 中每个扭曲的线程数），这为我们提供了一定的效率衡量标准（没有浪费的线程，没有数据存储的间隙）。
由于我们的线程策略是每个矩阵元素一个线程，因此我们必须使用 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 个线程用于此任务，矩阵有 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, and h）现在位于共享内存中。我们需要一种方法来按顺序加载它们，并执行所需的乘法和减法。
此时我们观察到两件事：
- each
M元素（辅因子）有一组不同的 4 个所需元素A
- each
M元素（辅因子）遵循相同的通用算术，给定四个任意元素A，我们将它们统称为X, Y, Z and W。算术是XY-ZW。我取第一个元素，将其乘以第二个元素，然后取第三个和第四个元素并将它们相乘，然后减去两个乘积。
由于所有 9 个辅助因子的一般操作顺序（上面的 2）都是相同的，因此我们只需要一种方法来安排 4 个所需矩阵元素的加载。此方法被编码到硬编码到示例中的负载模式中：
hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)
有9种负载模式，每种负载模式占用一个十六进制数，每个线程一种负载模式，即每个线程一种负载模式M矩阵元素（辅因子）。在特定的A矩阵，矩阵元素a, b, c等（已经）加载到共享内存中group偏移量为 0、1、2 等。给定线程的加载模式将允许我们生成组偏移量序列，需要检索以下矩阵元素A从它们在共享内存中的位置，按顺序使用来计算分配给该线程的辅因子。考虑线程 0 及其辅助因子ei-fh，负载模式如何0x7584对需要的模式进行编码以选择e, then i, then f, then h?
为此，我们有一个辅助函数getoff它采用加载模式，并连续（每次调用时）剥离索引。我第一次打电话getoff参数为0x7584，它“剥离”索引 4，返回该索引，并替换0x7584加载模式0x758以供下次使用。 4对应于e。下次我打电话的时候getoff with 0x758它“剥离”索引 8，返回该索引并替换0x758 with 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然后函数添加到该基地址以选择相关的输入元素。
行列式的计算
行列式的数值由下式给出：
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 个线程。因此，我们将在共享内存中选择 3 个（不再需要）位置来“发布”这些值。我们仍然需要位置 0、1、2 中的值，因为我们需要输入矩阵元素a, b, and 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;
