详细了解大量 3x3 矩阵的求逆算法

Posted

技术标签:

【中文标题】详细了解大量 3x3 矩阵的求逆算法【英文标题】:Understanding in details the algorithm for inversion of a high number of 3x3 matrixes 【发布时间】:2020-03-22 04:23:22 【问题描述】:

我遵循这个原始帖子:PyCuda code to invert a high number of 3x3 matrixes。 建议作为答案的代码是:

$ 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();
    

确实,sibase定义的索引unsigned sibase = (threadIdx.x / 9)*9;的作用是什么

还有,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 的初学者,但我不够专业,无法完全理解这个内核代码。

【问题讨论】:

【参考方案1】:

预赛

首先,了解 3x3 矩阵求逆的算法很重要,请参阅here(及以下)。

用于内核设计的一般方法是为每个线程分配一个矩阵结果元素。因此,每个矩阵需要 9 个线程。最终,每个线程将负责为每个矩阵计算 9 个数值结果之一。为了计算两个矩阵,我们需要 18 个线程,3 个矩阵需要 27 个线程。

辅助任务是决定线程块/网格的大小。这遵循典型的方法(总体问题大小决定了所需的线程总数),但我们将为线程块大小做出特定选择 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 的 9 个元素要计算,所以我们将为 M 的每个元素分配一个线程。 M 的每个元素都依赖于多个输入值(abc 等),所以我们首先将每个输入值(有 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 矩阵元素(abc、...)都已加载到共享内存中,我们可以开始计算M 中的辅因子。让我们关注一个特定的线程 (0) 及其辅因子 (ei-fh)。计算这个辅因子所需的所有矩阵元素(eifh)现在都在共享内存中。我们需要一种方法来按顺序加载它们,并执行所需的乘法和减法。

此时我们观察到两件事:

    每个M 元素(辅因子)都有一组不同的A 所需的4 个元素 每个M 元素(辅因子)遵循相同的通用算​​法,给定A 的四个任意元素,让我们将它们统称为XYZW。算术是XY-ZW。我取第一个元素,乘以第二个元素,然后取第三个和第四个元素并将它们相乘,然后减去两个乘积。

由于所有 9 个辅因子的一般操作顺序(上面的 2)是相同的,我们只需要一个方法来安排 4 个所需矩阵元素的加载。这种方法被编码到硬编码到示例中的负载模式中:

 hpat = (0x07584, 0x08172, 0x04251, 0x08365, 0x06280, 0x05032, 0x06473, 0x07061, 0x03140)

有 9 种负载模式,每种占用一个十六进制数量,每个线程一个负载模式,即每个M 矩阵元素(辅因子)一个负载模式。在特定的A 矩阵中,矩阵元素abc 等(已经)在 group 偏移量为 0、1、2 等处加载到共享内存中. 给定线程的加载模式将允许我们生成组偏移序列,需要从它们在共享内存中的位置检索A 的矩阵元素,以便按顺序使用来计算分配给该线程的辅因子。考虑线程 0 及其辅因子 ei-fh,负载模式 0x7584 如何编码所需的模式以选择 e,然后是 i,然后是 f,然后是 h

为此,我们有一个辅助函数getoff,它采用加载模式,并连续(每次调用它)剥离索引。我第一次使用0x7584 的参数调用getoff 时,它“剥离”索引4,返回它,并用0x758 替换0x7584 加载模式以供下次使用。 4 对应于e。下次我用0x758 调用getoff 时,它会“剥离”索引8,返回它,并用0x75 替换0x758。 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 中的值,因为我们需要输入矩阵元素 abc 来计算行列式。但是我们不再需要位置 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;

【讨论】:

非常感谢!你帮助我掌握了内核代码的不同过程和结构。你的解释很清楚!问候

以上是关于详细了解大量 3x3 矩阵的求逆算法的主要内容,如果未能解决你的问题,请参考以下文章

是否可以使用 SIMD 指令进行 3x3 矩阵求逆?

矩阵求特征值 和矩阵求逆 计算复杂度分析, 继续求助

线性的求逆元

python,c矩阵求逆问题记录及解决方案

python,c矩阵求逆问题记录及解决方案

python,c矩阵求逆问题记录及解决方案