玩转 ArrayFire：12 GFOR — 并行 for 循环

---

前言

在《玩转 ArrayFire：11 计时函数》中，我们已经了解到 ArrayFire 的计时函数，在这一篇中，我们将继续学习 ArrayFire 的并行 for 循环：GFOR ，即在 GPU 或设备上同时运行多个独立的循环。

---

一、概述

     GFOR 可以在 GPU 或设备上同时启动 for 循环的所有迭代，只要迭代是独立的。标准的 for 循环按顺序执行每个迭代，而 ArrayFire 的 gfor 循环则同时（并行地）执行每个迭代。ArrayFire 通过平铺所有循环迭代，然后在这些平铺上执行一次计算来实现这一点。

    你可以把 GFOR 看作是对你的代码执行自动向量化，例如，你编写了一个 gfor 循环，对 vector 的每个元素递增，但在后台，ArrayFire 会重写它，以并行地对整个 vector 进行操作。

	for (int i = 0; i < n; ++i)
	   A(i) = A(i) + 1;
	gfor (seq i, n)
	   A(i) = A(i) + 1;

    在后台，ArrayFire 会把你的代码改写成这个等价的、更快的版本：

	A = A + 1;

    注意：最好尽可能地向量化计算，以避免for循环和gfor循环中的开销。
    再看另一个例子，你可以用 for 循环运行FFT，或者你也可以“向量化”，简单地在一个 gfor 循环操作中完成所有操作：

	for (int i = 0; i < N; ++i)
	   A(span,span,i) = fft2(A(span,span,i)); // runs each FFT in sequence
	gfor (seq i, N)
	   A(span,span,i) = fft2(A(span,span,i)); // runs N FFTs in parallel

    实例化 gfor 循环有三种格式：

• gfor(var,n) ：创建一个序列 {0,1，…, n - 1} ；
• gfor(var,first,last) ：创建一个序列 {first, first+1, …, last}；
• gfor(var,first,incr,last) ：创建一个序列{first, first + inc, first + 2 * inc, …, last}

    所以下面所有的代码都代表了等价序列｛0, 1, 2, 3, 4｝：

	gfor (seq i, 5)
	gfor (seq i, 0, 4)
	gfor (seq i, 0, 1, 4)

    有如下例子进行进一步说明：

	array A = constant(1, n, n);
	array B = constant(1, 1, n);
	gfor (seq k, 0, n-1) {
	   B(span, k) = sum(A(span, k) * A(span,k));  // inner product
	}

	array A = randu(n,m);
	array B = constant(0,n,m);
	gfor (seq k, 0, m-1) {
	   B(span,k) = fft(A(span,k));
	}

二、如何使用

1. 在 GFOR 中调用用户函数

    如果你定义了一个想在 GFOR 循环中调用的函数，那么这个函数必须满足本篇中描述的所有条件（见本篇第三节），以便能够像预期的那样工作。

2. 迭代器

    迭代器可以用在表达式中：

	A = constant(1,n,n,m);
	B = constant(1,n,n);
	gfor (seq k, m)
		A(span, span, k) = (array(k) + 1) * B + sin(array(k) + 1);  // expressions

3. 操作下标

    支持更复杂的下标：

	A = constant(1,n,n,m);
	B = constant(1,n,10);
	gfor (seq k, m)
	  A(span,seq(10),k) = k*B;  // subscripting, seq(10) generates index [0,9]

    迭代器可以与下标中的算术相结合：

	array A = randu(n,m);
	array B = constant(1,n,m);
	gfor (seq k, 1, m-1)
	  B(span,k) = A(span,k-1);
	A = randu(n,2*m);
	B = constant(1,n,m);
	gfor (seq k, m)
	  B(span,k) = A(span,2*(k+1)-1);
	A = randu(n,2*m);
	B = constant(1,n,m);
	gfor (seq k, m)
	  B(span,k) = A(span,floor(k+.2));

4. In-Place 计算

    在某些情况下，GFOR 的行为与典型的顺序for循环不同。例如，只要访问是独立的，您就可以在适当的地方读取和修改结果。

	A = constant(1,n,n);
	gfor (seq k, n)
	  A(span,k) = sin(k) + A(span,k);

5. 随机数据生成

    随机数据应该总是在 GFOR 循环之外生成。这是因为 GFOR 只经过循环体一次。因此，在循环体中对 randu() 的任何调用都将导致将相同的随机矩阵赋值给循环的每次迭代。

    gfor(seq ii, n) {
        array A = randu(3, 1);
        B(span, ii) = A;
    }
	af_print(B);
	//B
	//[3 3 1 1]
	//    0.6010     0.6010     0.6010
	//    0.0278     0.0278     0.0278
	//    0.9806     0.9806     0.9806

    这可以通过在循环外部引入随机数生成来纠正，如下所示：

	array A = randu(3,n);
	gfor (seq ii, n)
	  B(span,ii) = A(span,ii);
	af_print(B);

    这是一个简单的示例，但是演示了在大多数情况下应该在循环外部预先分配随机数的原则。

三、限制条件

     GFOR 的初步实现有以下限制：

1. 迭代独立

    循环体最重要的特性是每个迭代必须独立于其他迭代。请注意，访问单独迭代的结果会产生未定义的行为。

	array B = randu(3);
	gfor (seq k, n)
	  B = B + k; // bad

2. 条件语句限制

    循环体中没有条件语句（即没有分支）。然而，您通常可以找到克服这一限制的方法。考虑以下两个例子：
    示例 1 ：

	A = constant(1,n,m);
	gfor (seq k, n) {
	 if (k > 10) A(span,k) = k + 1; // bad
	}

    但是，您可以使用一些技巧来克服这个限制，方法是将条件语句表示为逻辑值的乘法。例如，上面的代码块可以转换为如下方式运行，不会出错：

	gfor (seq k, m) {
	  array condition = (k > 1); // good
	  A(span,k) = (!condition).as(f32) * A(span,k) + condition.as(f32) * (k + 1);
	}

    示例 2 ：

	array A = constant(1,n,n,m);
	array B = randu(n,n);
	gfor (seq k, 4) {
	  if ((k % 2) != 0)
	     A(span,span,k) = B + k;
	  else
	     A(span,span,k) = B * k;
	}

    相反，您可以对相同的数据进行两次传递，每次传递执行一个分支。

	A = constant(1,n,n,m);
	B = randu(n);
	gfor (seq k, 0, 2, 3)
	  A(span,span,k) = B + k;
	gfor (seq k, 1, 2, 3)
	  A(span,span,k) = B * k;

3. 嵌套循环

    不支持在 GFOR 循环中嵌套 GFOR 循环。你可以交错使用 for 循环，只要它们完全独立于 GFOR 迭代器。

	gfor (seq k, n) {
	  gfor (seq j, m) { // bad
	  // ...
	  }
	}

    支持在 GFOR 循环中嵌套FOR循环，只要 GFOR 迭代器不在 FOR 循环迭代器中使用，如下所示：

	gfor (seq k, n) {
	  for (int j = 0; j < (m+k); j++) { // bad
	  // ...
	  }
	}
	gfor (seq k, n) {
	  for (int j = 0; j < m; j++) { // good
	  //...
	  }
	}

    完全支持在for循环中嵌套gfor循环：

	for (int j = 0; j < m; j++) {
	  gfor (seq k, n) { // good
	  //  ...
	  }
	}

4. 逻辑索引限制

    不支持如下的逻辑索引：

	gfor (seq i, n) {
	  array B = A(span,i);
	  array tmp = B(B > .5); // bad
	  D(i) = sum(tmp);
	}

    问题是，每个 GFOR 都有不同数量的元素，这是GFOR无法处理的。类似于条件语句的解决方案，可以使用掩码算法：

	gfor (seq i, n) {
	  array B = A(span,i);
	  array mask = B > .5;
	  D(i) = sum(mask .* B);
	}

    支持使用标量和逻辑掩码的子赋值：

	gfor (seq i, n) {
	  a = A(span,i);
	  a(isnan(a)) = 0;
	  A(span,i) = a;
	}

四、内存考虑

    由于每个计算都是针对所有迭代器值并行进行的，因此需要有足够的内存来同时执行所有迭代。如果问题超出了内存，它将触发“内存不足”错误。
    你可以通过将 GFOR 循环分解成段来绕过 GPU 或设备的内存限制；然而，您可能需要考虑使用一个更大的内存GPU或设备。

	// BEFORE
	gfor (seq k, 400) {
	  array B = A(span,k);
	  C(span,span,k) = matmulNT(B * B);  // outer product expansion runs out of memory
	}
	// AFTER
	for (int kk = 0; kk < 400; kk += 100) {
	  gfor (seq k, kk, kk+99) { // four batches of 100
	     array B = A(span,k);
	     C(span,span,k) = matmulNT(B, B); // now several smaller problems fit in card memory
	  }
	}

---