基于WebGL的GPGPU指南 - 矩阵乘法

[llwanghong]

整理译文的目录章节如下：

• 概述
• 各阶段实现
• 初始化矩阵
• 减速带
• 矩阵乘法

矩阵乘法

我们已经完成了矩阵的初始化。进一步可以做的最简单但有趣的操作是矩阵平方。

$$\textbf{R}_{ij} = \textbf{M}_{ij}^2 = \sum_{k=1}^{n}\textbf{M}_{ik}\textbf{M}_{kj}$$

其中 $\textbf{M}$ 是一个 $\mathbf{n{\times}n}$ 矩阵. 计算这个总和将展示如何在纹理中处理数据的技巧。

在深入研究代码之前，我们再来看一下预期的GPGPU应用程序的结构。

WebGL GPGPU实现的基本结构

我们将从画布（canvas）开始自底向上地构建整个程序。要计算结果矩阵 $\textbf R$ 的每个元素 $\textbf{R}_{ij}$，其中 $\textbf R = \textbf M^2$，所以矩阵 $\textbf{R}$ 和 $\textbf{M}$ 一样为 $\mathbf{n{\times}n}$ 矩阵。这样画布大小也同样是 $\mathbf{n{\times}n}$.

  var gpgpUtility;
  var matrixColumns;
  var matrixRows;

  matrixColumns = 128;
  matrixRows = 128;
  gpgpUtility = new vizit.utility.GPGPUtility(matrixColumns, matrixRows);

纹理的 $x$ 坐标 $s$，从 $0$ 变化到 $1$，而画布从 $0$ 变化到 $canvas.width$。类似的，纹理的 $y$ 坐标从 $0$ 变化到 $1$，而画布从 $0$ 变化到 $canvas.height$。我们因此可得到纹理的一些重要结论信息。

纹理元素之间的 $x$ 间距 $\delta_s$ 和 $y$ 间距 $\delta_t$，分别为：

$$\delta_s = \frac{1}{canvas.width}$$

$$\delta_t = \frac{1}{canvas.height}$$

在片段着色器中，我们可以访问到纹理坐标，从这个坐标可以准确地知道正在写入哪个纹理元素。具体来说，

$$i = floor\left( canvas.width \times s \right)$$

$$j = floor\left( canvas.height \times t \right)$$

几何形状再次是完全覆盖画布的标准三角形。由这个画布加几何形状组合生成的片段为矩阵 $\textbf R$ 每个元素创建了一个片段。

我们将初始化步骤中的渲染目标 $m$ 作为矩阵乘法阶段的输入纹理。将上一步骤中的目标纹理作为下一个步骤的输入纹理是一种常见的做法。

 /**
   * Runs the program to do the actual work. On exit the framebuffer &
   * texture are populated with the square of the input matrix, m. Use
   * gl.readPixels to retrieve texture values.
   *
   * @param m        {WebGLTexture} A texture containing the elements of m.
   * @param mSquared {WebGLTexture} A texture to be incorporated into a fbo,
   *                                the target for our operations.
   */
  this.square = function(m, mSquared)
  {
    var m2FrameBuffer;

    // Create and bind a framebuffer
    m2FrameBuffer = gpgpUtility.attachFrameBuffer(mSquared);

    gl.useProgram(program);

    gpgpUtility.getStandardVertices();

    gl.vertexAttribPointer(positionHandle,     3, gl.FLOAT, gl.FALSE, 20, 0);
    gl.vertexAttribPointer(textureCoordHandle, 2, gl.FLOAT, gl.FALSE, 20, 12);

    gl.activeTexture(gl.TEXTURE0);
    gl.bindTexture(gl.TEXTURE_2D, m);
    gl.uniform1i(textureHandle, 0);

    gl.drawArrays(gl.TRIANGLE_STRIP, 0, 4);
  };

$$\textbf{R}_{ij} = \textbf{M}_{ij}^2 = \sum_{k=1}^{n}\textbf{M}_{ik}\textbf{M}_{kj}$$

我们设置几何形状以便每个 $\textbf{R}_{ij}$ 元素都精确对应一个片段。在着色器中，对 $\textbf{R}_{ij}$ 的赋值对应到对 $\textbf{gl\_FragColor}$ 的赋值，相当于设置了 $\textbf{m2}$ 纹理中 $(s, t)$ 位置的纹素。

方程右侧 $\textbf{M}$ 相关元素将从 $\textbf m$ 纹理来读取。其中 $\textbf i$ 和 $\textbf j$ 即为计算得到纹素坐标 $\textbf s$ 和 $\textbf t$， $\textbf k$ 如何算呢？

$\textbf k$ 从 $\textbf 0$ 变化到 $\textbf n$，其中 $\textbf n$ 为矩阵尺寸，本例中为128。转换到纹理坐标体系，则从 $\textbf 0$ 变化到 $\textbf 1$，步长为 $\mathbf{\frac{1}{n}}$ 。

  #ifdef GL_FRAGMENT_PRECISION_HIGH
  precision highp float;
  #else
  precision mediump float;
  #endif

  uniform sampler2D m;

  varying vec2 vTextureCoord;

  void main()
  {
    float i, j;
    float value = 0.0;

    i = vTextureCoord.s;
    j = vTextureCoord.t;

    for(float k=0.0; k<128.0; ++k)
    {
      value += texture2D(m, vec2(i, k/128.0)).r * texture2D(m, vec2(k/128.0, j)).r;
    }
    gl_FragColor.r = value;
  }

GLSL（特别是WebGL使用的GLSL版本）有一个有趣的要求，即循环的界限必须是固定的。根据GLSL规范附录A控制流（Control Flow）章节，循环的最大迭代次数需要在编译时确定。这意味着我们需要为不同大小的矩阵编写不同的着色器。

另一种方法是设置一个非常大的循环索引，并在完成时跳出循环，类似于下面的方法：

  uniform sampler2D m;
  uniform float mSize;

  varying vec2 vTextureCoord;

  void main()
  {
    float i, j;
    float value = 0.0;

    i = vTextureCoord.s;
    j = vTextureCoord.t;

    for(float k=0.0; k<2048.0; ++k)
    {
      if (k >= mSize)
      {
        break;
      }
      value += texture2D(m, vec2(i, k/mSize)).r * texture2D(m, vec2(k/mSize, j)).r;
    }
    gl_FragColor.r = value;

虽然这种方法可以让你拥有一个适用于多种情况的单一着色器，但仍存在一些问题。许多GPU，尤其是老款GPU，会展开循环，并执行所有真假条件路径，然后丢弃不需要的结果。这会导致意料之外的大内存消耗以及糟糕的执行时间。

问题： 你将如何修改这个算法来实现两个不同尺寸矩阵相乘？

下面表格是这段代码的运行结果。同时将GPU计算结果与JavaScript的64位计算结果进行了比较。通常情况下这些结果不会匹配。

GPGPU结果与JavaScript 64位结果相比较

坐标	GPU结果	正确期望值	相对误差
(1, 1)	8819015680	8819016128	5.1e-8
(10, 12)	81986330624	81986337536	8.4e-8
(20, 30)	163327934464	163327923840	6.5e-8
(100, 100)	814771535872	814771692800	1.9e-7
(101, 101)	822925459456	822925428928	3.7e-8
(127, 127)	1035012341760	1035012424256	8.0e-8

这就引出了关于GPU计算的一个重要问题。WebGL中你不会看到64位（双精度）浮点数。最多只能看到32位（单精度）甚至24位或16位浮点数。针对问题和示例选择足够的精度非常重要，就像许多教学模拟的情况一样。