算法导论26 Problems 答案

发表于 分类于 阅读次数： 本文字数： 12k 阅读时长 ≈ 11 分钟

26-1

a

使用得到P-MAT-VEC-RECURSIVE类似的方法对SUM-ARRAYS进行修改，同样可以得到其基于递归的并行版本SUM-ARRAYS-RECURSIVE：

1
2
3
4
5
6
7
8
SUM-ARRAYS-RECURSIVE(A, B, C, i, i')
  if i == i'
    C[i] = A[i] + B[i]
  else
    mid = ⌊(i + i') / 2⌋
    spawn SUM-ARRAYS-RECURSIVE(A, B, C, i, mid)
    spawn SUM-ARRAYS-RECURSIVE(A, B, C, mid + 1, i')
    sync

可见这个算法的工作量为\(T_1(n)=\Theta(n)\)。由于其递归深度达到\(\Theta(\lg n)\)，并且只需要花费\(\Theta(1)\)的时间即可完成，因此其持续时间为\(T_{\infty}(n)=\Theta(\lg n)\cdot\Theta(1)=\Theta(\lg n)\)。因此的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n/\lg n)\)。

b

如果\(grain\text{-}size=1\)，那么\(r=n\)，SUM-ARRAYS'中的for循环将会串行地执行\(n\)次，每次调用ADD-SUBARRAY都只花费\(\Theta(1)\)的时间。因此持续时间\(T_{\infty}(n)=n\cdot \Theta(1)=\Theta(n)\)。可以知道其工作量为\(T_1(n)=\Theta(n)\)，因此其并行度为\(T_1(n)/T_{\infty}(n)=\Theta(1)\)。

c

令\(g=grain\text{-}size\)。执行单次ADD-SUBARRAY所需要的时间为\(O(g)\)。在SUM-ARRAYS'中，一共需要执行\(n/g\)次ADD-SUBARRAY，但是for循环并没有带有parrllel关键字，也就是说，这些调用ADD-SUBARRAY的线程是按顺序启动的，它们并非同时启动的，因此这里仍然需要花费\(O(n/g)\)的时间完成这个for循环。

因此，SUM-ARRAYS'需要花费\(O(g+n/g)\)的时间完成。令\(f(g)=g+n/g\)，那么有\(f'(g)=1-n/g^2\)。令\(f'(g)=00\)，即\(g=\sqrt{n}\)时，SUM-ARRAYS'只需要花费\(O(\sqrt{n})\)的时间就可以完成。

26-2

a

修改后的代码由P-MATRIX-MULTIPLY-RECURSIVE'所示，它消去了临时矩阵\(D\)的存在。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
P-MATRIX-MULTIPLY-RECURSIVE'(A, B, C, n)
  if n == 1
    c11 = c11 + a11 · b11
    return
  partition A, B, C, and D into n/2 × n/2 submatrices A11, A12, A21, A22; B11, B12, B21, B22; and C11, C12, C21, C22; respectively
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A11, B11, C11, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A11, B12, C12, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A21, B11, C21, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A21, B12, C22, n / 2)
  sync
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A12, B21, C11, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A12, B22, C12, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A22, B21, C21, n / 2)
  spawn P-MATRIX-MULTIPLY-RECURSIVE'(A22, B22, C22, n / 2)
  sync

可见，它的串行投影是MATRIX-MULTIPLY-RECURSIVE，因此其工作量为\(\Theta(n^3)\)。

b

对于其工作量\(T_1(n)\)，可以给出其递推式\(T_1(n)=8T_1(n/2)+\Theta(n^2)\)，因此得到\(T_1(n)=\Theta(n^3)\)。

对于其持续时间，除了递归调用，P-MATRIX-MULTIPLY-RECURSIVE'没有再做其它工作。因此可以给出其递推式\(T_{\infty}(n)=2T_{\infty}(n/2)+\Theta(1)\)，最终通过主定理可以得知\(T_{\infty}(n)=\Theta(n)\)。

也就是说，算法P-MATRIX-MULTIPLY-RECURSIVE'的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^2)\)。

c

忽略掉\(\Theta\)符号后，对于\(1000\times1000\)的矩阵而言，调用P-MATRIX-MULTIPLY-RECURSIVE'的并行度为\(10^6\)。相比于调用P-MATRIX-MULTIPLY-RECURSIVE得到\(10^7\)的并行度，计算机依旧不会有这么多的处理器，但是在算法的运行时间的常数却明显减少了，因此这种权衡是值得的。

26-3

a

对LU-DECOMPOSITION修改后得到的并行化版本为P-LU-DECOMPOSITION：

1
2
3
4
5
6
7
8
9
10
11
12
13
P-LU-DECOMPOSITION(A, n)
  let L and U be new n × n matrices 
  initialize U with 0s below the diagonal
  initialize L with 1s on the diagonal and 0s above the diagonal
  for k = 1 to n
    u_{kk} = a_{kk}
    parallel for i = k + 1 to n
      l_{ik} = a_{ik} / a_{kk}
      u_{ki} = a_{ki}
    parallel for i = k + 1 to n
      parallel for j = k + 1 to n
        a_{ij} = a_{ij} − l_{ik} * u_{kj}
  return L and U

也就是说，所有内循环都可以进行并行，但是最外层的那一层循环不允许并行。第一个内循环能够并行的原因是：它只读取\(a\)矩阵的元素，并修改\(L,U\)对应行和对应列的元素，并不会进行多次存取。第二个双重内循环能够并行的原因是第\(i\)轮二层循环，第\(j\)轮三次循环中，他只会读取\(L,U\)矩阵的值并且修改\(a_{ij}\)的值，不会造成冲突。外层循环不能是并行的原因是按照LU分解的定义，后一次的运行依赖于前一次的结果。

P-LU-DECOMPOSITION的串行投影是LU-DECOMPOSITION，因此其工作量为\(T_1(n)=\Theta(n^3)\)。

对于P-LU-DECOMPOSITION的持续时间，第一层的第一次内循环花费的时间是\(\Theta(\lg(n-k))=\Theta(\lg n)\)，因为转化成spawn ... sync结构后就变成了分支法，然后进行递归。第二次的两个双重嵌套寻呼按的原因相同，其花费的时间为\(\Theta(\lg(n-k))+\Theta(\lg(n-k))=\Theta(\lg n)\)。因此，外层for循环需要\(\Theta(\lg n)\)的时间完成一次。最终，P-LU-DECOMPOSITION的持续时间为\(T_{\infty}(n)=n\cdot \Theta(\lg n)=\Theta(n\lg n)\)。

最终我们可以计算出P-LU-DECOMPOSITION的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^2/\lg n)\)。

b

对LUP-DECOMPOSITION修改后得到的并行化版本为P-LUP-DECOMPOSITION：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
P-LUP-DECOMPOSITION(A, n)
  let π[1 : n] be a new array
  parallel for i = 1 to n
    π[i] = i
  for k = 1 to n
    p = 0
    parallel for i = k to n
      if |a_{ik}| > p
        p = |a_{ik}|
        k' = i
    if p == 0
      error "singular matrix"
    exchange π[k] with π[k']
    parallel for i = 1 to n
      exchange a_{ki} with a_{k'i} 
    parallel for i = k + 1 to n
      a_{ik} = a_{ik} / a_{kk}
      parallel for j = k + 1 to n
        a_{ij} = a_{ij} − a_{ik} * a_{kj}

第2行的for循环可以并行，这是显而易见的。第6-9行的for循环是可以并行的，因为这是寻找一个最大值，在这个过程中，\(n\)个元素构成一棵\(n\)个叶节点的二叉树，对于内部节点，它是两个子节点中的最大值，由于这棵树的深度为\(\Theta(\lg n)\)，因此这段循环的持续时间为\(\Theta(\lg n)\)。第13行的for循环用于交换两行元素，显而易见是可以并行的。对于第15行和17行的for循环，第16行仅仅是对\(a_{ik}\)进行更新，此后都是对\(a_{ij}(j>k)\)中的所有元素进行更新。至于先对哪一行操作都是没有关系的，因此第15行和17行的for循环它们都可以并行。外面针对\(k\)的for循环不能并行，因为会导致数据读取冲突。

P-LUP-DECOMPOSITION的串行投影是LUP-DECOMPOSITION，因此其工作量为\(T_1(n)=\Theta(n^3)\)。

对于P-LUP-DECOMPOSITION的持续时间，分析和P-LU-DECOMPOSITION非常相似。第6，13，15，17行的for循环都可以进行，因此这部分只使用\(\Theta(\lg n)\)的时间即可完成。最终，P-LUP-DECOMPOSITION的持续时间为\(T_{\infty}(n)=n\cdot \Theta(\lg n)=\Theta(n\lg n)\)。

最终我们可以计算出P-LUP-DECOMPOSITION的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^2/\lg n)\)。

c

对LUP-SOLVE修改后得到的并行化版本为P-LUP-SOLVE：

1
2
3
4
5
6
7
8
9
10
11
12
13
P-LUP-SOLVE(L, U, π, b, n)
  let x and y be new vectors of length n
  for i = 1 to n
    val = 0
    parallel for j = 1 to i - 1
      val = val + l_{ij} * y_j
    y_i = b_{π[i]} - val
  for i = n downto 1
    val = 0
    parallel for j = i + 1 to n
      val = val + u_{ij} * x_j
    x_i = (y_u - val) / u_{ii}
  return x

第4和9行的for循环可以并行，只需要将长度为\(n\)的序列构造出一棵\(n\)个叶子节点的完全二叉树即可，对于内部节点，它是两个子节点中的值之和。因此，这两个内部for循环都能够在\(\Theta(\lg n)\)的时间内完成。对于第2和7行的for循环，因为后面计算\(x,y\)值需要依赖前面已经计算出的\(x,y\)值。

最终，P-LUP-SOLVE的串行投影即为LUP-SOLVE，因此其工作量为\(T_1(n)=\Theta(n^2)\)。按照上面的结论，可以得到持续时间\(T_{\infty}(n)=n\cdot \Theta(\lg n)=\Theta(n\lg n)\)。最终我们可以计算出P-LUP-SOLVE的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n/\lg n)\)。

d

一个对正定矩阵求逆的并行算法由P-MATRIX-INVERSION-POSITIVE-DEFINITE给出。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
P-MATRIX-INVERSION-POSITIVE-DEFINITE(A, n)
  partition A into n/2 × n/2 submatrices B, CT, C, D; respectively
  let W[1 : n, 1 : n], X[1 : n, 1 : n], Y[1 : n, 1 : n], Z[1 : n, 1 : n], S[1 : n, 1 : n] be new matrices
  B' = P-MATRIX-INVERSION-POSITIVE-DEFINITE(B, n / 2)
  W = P-MATRIX-MULTIPLY(C, B', W, n / 2)
  WT = W
  P-TRANSPOSE(WT, n / 2)
  P-MATRIX-MULTIPLY(W, CT, X, n / 2)
  S = D - X
  S' = P-MATRIX-INVERSION-POSITIVE-DEFINITE(S, n / 2)
  P-MATRIX-MULTIPLY(S', W, Y, n / 2)
  YT = Y
  P-TRANSPOSE(YT, n / 2)
  P-MATRIX-MULTIPLY(WT, y, Z)
  R = B' + Z
  A' = [[R, -YT], [-Y, S']]
  return A'

按照第28.2章的结论以及不等式28.15，可以知道其工作量\(I_1(n)\le 2I_1(n/2)+4M_1(n/2)+O(n^2)\)，最终得到\(I_1(n)=O(M(n))\)。

在并行版本中，由题目26.1-8，26.2-3和题26-1-a可知，矩阵的转置、加法和乘法都可以在\(\Theta(\lg n)\)时间内完成。在P-MATRIX-INVERSION-POSITIVE-DEFINITE中，一共进行了\(2\)次矩阵转置，\(4\)次矩阵加减法，\(4\)次矩阵乘法。同样的，它们还进行了\(2\)次大小为\(n/2\times n/2\)矩阵的求逆。可见，这些步骤都是依赖于前一步的结果，因此步骤之间不能直接并行。我们可以写出它的持续时间\(I_{\infty}(n)\)满足\(I_{\infty}(n)=2I_{\infty}(n/2)+10\cdot\Theta(\lg n)\)，因此最终得到\(I_{\infty}(n)=\Theta(n)\)。

最终我们可以得到它的并行度为\(I_1(n)/I_{\infty}(n)=O(M(n)/n)\)。

26-4

a

设计的P-REDUCE如下所示：

1
2
3
4
5
6
7
8
9
P-REDUCE(x, i, j)
  if i == j
    return x[i]
  else
    mid = ⌊(i + j) / 2⌋
    l = spawn P-REDUCE(x, i, mid)
    r = spawn P-REDUCE(x, mid + 1, j)
    sync
    return l ⊗ r

不难发现它的串行投影即为REDUCE，因此它的工作量为\(T_1(n)=\Theta(n)\)。由于这棵进行搜索的二叉树的深度为\(\Theta(\lg n)\)，因此这个算法的持续时间为\(T_{\infty}(n)=\Theta(\lg n)\)。

b

这个算法相当于是以\(k=1,2,\dots,n\)对P-REDUCE都进行了一次调用。因此P-SCAN-1的工作量为\(\displaystyle{T_1(n)=\sum_{k=1}^n\Theta(k)=\Theta(n^2)}\)。

接下来求解P-SCAN-1的持续时间。由于对P-SCAN-1对P-REDUCE的间接调用呈树形，假设伪代码中，第\(i\)次对P-REDUCE的调用的持续时间为\(iter(i)\)，那么有\(T_{\infty}(n)=\Theta(\lg n)+\max\{iter(i):1\le i\le n\}\)，按照题目26-4-a的结论，有\(T_{\infty}(n)=\Theta(\lg n)\)。

因此P-SCAN-1的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^2/\lg n)\)。

c

我们可以使用归纳法进行证明其正确性。即证明，调用P-SCAN-2-AUX(i, j)，对于\(\forall k\in[i,j]\)，它都正确计算出了\(y[k]=x[i]\otimes x[i+1]\otimes\dots\otimes x[k]\)。

对于基本情况：\(i=j\)，P-SCAN-2-AUX都正确计算出了\(y\)数组，这由P-SCAN-2-AUX的前两行明显可知。

当\(i<j\)时，假设对于区间\([i,j]\)内的所有真子区间\([i',j']\)，即\(i\le i'\le j'\le j\)，且\(i=i',j=j'\)不同时成立，P-SCAN-2-AUX都计算出了正确的结果。那么第三行得到了一个\(k=\lfloor(i+j)/2\rfloor\)。对于第4和5行的代码，它们分别对\(y\)数组的区间\([i,k]\)和\([k+1,j]\)进行写入，并且读取的内容也不相交，因此这两行代码不会构成竞争。在第6行结束后，\(y\)数组满足：如果\(p\le k\)，那么\(y[p]=x[i]\otimes x[i+1]\otimes\dots\otimes x[p]\)，否则\(y[p]=x[k+1]\otimes x[k+2]\otimes\dots\otimes x[k]\)。对于\(p>k\)的情况，第8行将\(y[p]\)重新赋值成\(y[k]\otimes y[p]\)。这个步骤完成后，\(\forall p\in[i,j]\)，都有\(y[p]=x[i]\otimes x[i+1]\otimes\dots\otimes x[p]\)。因此P-SCAN-2-AUX是正确的。

接下来考虑P-SCAN-2-AUX的工作量\(T_1(n)\)，消去最后的for循环中的parallel关键字后，那么除去递归部分，它的运行时间是\(\Theta(n)\)。因此可以写出\(T_1(n)=2T_1(n/2)+\Theta(n)\)，从而得到\(T_1(n)=\Theta(n\lg n)\)。

接下来考虑其持续时间\(T_{\infty}(n)\)。可见出了递归部分，其余部分仍然需要\(\Theta(\lg n)\)进行。因此有\(T_{\infty}(n)=T_{\infty}(n/2)+\Theta(\lg n)\)，最终得到\(T_{\infty}(n)=\Theta(\lg^2n)\)。

因此P-SCAN-2的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n/\lg n)\)。

d

P-SCAN-UP第8行填入的是$right ⊗ t[k]；P-SCAN-DOWN的第5行填入v，第6行填入v ⊗ t[k]。

需要注意的是，\(t[k]\)存储的是当一个区间\([i,j]\)被尽量均匀地划分成两个区间\([i,k],[k+1,j]\)后，\(t[k]\)就是记录区间\([i,j]\)前面一半元素（即\([i,k]\)）的元素之和\(t[k]=x[i]\otimes x[i+1]\otimes\dots\otimes x[k]\)。

首先证明，调用时每个\(t[k]\)最多只会被赋值一次（因此，在执行P-SCAN-DOWN时不会造成任何冲突）。只有调用P-SCAN-UP，当\(i<j\)时，才会对\(t[k]\)进行赋值，在此之后更深的递归讲不会再对\(t[k]\)进行赋值，原因如下。如果递归的区间是左半子区间\([i,k]\)，那么对于所有\([i,k]\)的长度大于等于\(2\)的（等于\(1\)则不会进入到这个分支）真子区间\([i',k']\)，都有\(\lfloor(i'+k')/2\rfloor<k\)，因此\(t[k]\)不会被重复赋值；如果递归的是右子区间\([k+1,j]\)，那么它的访问和读写只会在这个区间内进行，更不会对\(t[k]\)进行访问。因此，\(t[k]\)确实能够正确记录区间\([i,j]\)的信息。根据P-SCAN-DOWB的第5行代码，就可以知道\(t[k]\)记录的是区间i,j的左半区间的元素之和。并且，P-SCAN-UP的返回值是\([i,j]\)这个区间的所有元素之和，因此结论成立。

接下来证明每次调用P-SCAN-DOWN(v, x, t, y, i, j)时，总满足\(v=x[1]\otimes x[2]\otimes\dots\otimes x[i-1]\)。同样使用归纳法来证明。在P-SCAN-3调用P-SCAN-DOWN时，有\(v=x[1],i=2\)，因此基本情况下是成立的。P-SCAN-DOWN首先调用P-SCAN-DOWN(v, x, t, y, i, k)，由于参数\(i\)没有变化，因此\(v\)仍然使用原来的\(v\)，原结论成立；然后调用P-SCAN-DOWN(v ⊗ t[k], x, t, y, k + 1, j)，由于此时\(t[k]=x[i]\otimes x[i+1]\otimes\dots\otimes x[k]\)，因此有\(v\otimes t[k]=x[1]\otimes x[2]\otimes\dots\otimes x[k]\)。因此第二次调用时同样满足题目的条件。因此当P-SCAN-DOWN进入第2行后，y[i] = v ⊗ x[i]则是\(y[i]=x[1]\otimes x[2]\otimes\dots\otimes x[i]\)，\(y[i]\)被正确地计算出来。

因此，算法P-SCAN-3是正确的。

e

可以发现，这棵树的节点数仍然是\(\Theta(n)\)，因此P-SCAN-3的工作量\(T_1(n)=\Theta(n)\)。

接下来首先考虑P-SCAN-UP。由于其每次划分出了两个长度尽量均等的区间，并且其余处理部分都只要\(\Theta(1)\)的时间，因此P-SCAN-UP这段代码的持续时间满足\(T_{\infty}(n)=T_{\infty}(n/2)+\Theta(1)\)，从而得到\(T_{\infty}(n)=\Theta(\lg n)\)。P-SCAN-DOWN和P-SCAN-UP的结构基本相同，因此对其分析也一样。

最终，P-SCAN-3的持续时间为\(T_{\infty}(n)=\Theta(\lg n)\)。

因此，P-SCAN-3的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n/\lg n)\)。

f

修改后的P-SCAN-3由P-SCAN-3'给出。代价是并发度将会降低。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
P-SCAN-3'(x, n)
  let y[1 : n] be a new array
  y[1] = x[1]
  if n > 1
    P-SCAN-3-AUX(x[1], x, y, 2, n)
  return y

P-SCAN-SUM-ARRAYS(x, i, j)
  if i == j
    return x[i]
  else
    k = ⌊(i + j) / 2⌋
    l = spawn P-SCAN-SUM-ARRAYS(x, t, i, k)
    r = P-SCAN-SUM-ARRAYS(x, t, k + 1, j)
    return l ⊗ r

P-SCAN-3-AUX(v, x, y, i, j)
  if i == j
    y[i] = v ⊗ x[i]
  else
    k = ⌊(i + j) / 2⌋
    spawn P-SCAN-3-AUX(v, x, y, i, k)
    t = P-SCAN-SUM-ARRAYS(x, i, k)
    P-SCAN-3-AUX(v ⊗ t, x, y, k + 1, j)
    sync

\(\star\) g

1
2
3
4
5
6
7
8
P-SCAN-4(x, n)
  l = 1
  while l <= n
    parallel for i = 1 to n by l * 2
      parallel j = 0 to l - 1
        if i + l + j <= n
          x[i + l + j] = x[i + l + j] + x[i + l - 1]
    l = l * 2

这个算法的基本思想是，将\(x\)分成\(\lceil n/(2l)\rceil\)块，每一块的长度为\(2l\)（最后一块不足\(2l\)也以一块记），并且\(l\)是\(2\)的幂。一开始\(l=1\)，将前半块的最后一个元素\(x[i+l-1]\)添加到后半块\(x[i+l+j](0\le j<l)\)中的每一个元素。从而最终完成后，每一块内部都是对应值的前缀和。可见两个parallel for循环不会导致冲突，因此这个算法的持续时间为\(T_{\infty}(n)=\Theta(\lg^2n)\)。

h

将字符串中的(视为\(+1\)，')'视为\(-1\)，得到一个数组\(x\)，求出其前缀和\(y\)。那么一个括号字符串是合法的，当且仅当\(\forall i\in[1,n],y_i\ge 0\)，并且有\(y_n=0\)。通过对P-SCAN-3进行改造，我们可以得到一个在\(\Theta(\lg n)\)时间内判断一个括号字符串是否合法的程序PARENTHESES-IS-WELL-FORMED：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
PARENTHESES-IS-WELL-FORMED(s, n)
  let y[0 : n] and t[1 : n] be new arrays
  y[0] = 0
  P-SCAN-UP'(s, t, 1, n)
  mn = P-SCAN-DOWN'(0, s, t, y, 1, n)
  return mn >= 0 and y[n] == 0

P-SCAN-UP'(s, t, i, j)
  if i == j
    if s[i] == '('
      return 1
    else
      return -1
  else
    k = ⌊(i + j) / 2⌋
    t[k] = spawn P-SCAN-UP'(s, t, i, k)
    right = P-SCAN-UP'(s, t, k + 1, j)
    sync
    return t[k] + right

P-SCAN-DOWN'(v, s, t, y, i, j)
  if i == j
    y[i] = v + x[i]
    return y[i]
  else
    k = ⌊(i + j) / 2⌋
    l = spawn P-SCAN-DOWN (v, x, t, y, i, k)
    r = P-SCAN-DOWN(v + t[k], x, t, y, k + 1, j)
    sync
    return min{l, r}

改造后的程序还返回了前缀和数组中的最小值。

26-5

a

不失一般性，这里假设矩阵\(A\)的大小\(n\)是\(2\)的幂次。基于等式26.9，那么对矩阵\(A\)的填充并行分治算法由SIMPLE-STENCIL给出。

1
2
3
4
5
6
7
8
9
10
SIMPLE-STENCIL(A, n)
  if n == 1
    generate the value of a11
    return
  partition A into n/2 × n/2 submatrices A11, A12, A21, A22; respectively
  SIMPLE-STENCIL(A11, n / 2)
  spawn SIMPLE-STENCIL(A12, n / 2)
  spawn SIMPLE-STENCIL(A21, n / 2)
  sync
  SIMPLE-STENCIL(A22, n / 2)

可见这个算法只是为了填充原来矩阵的所有元素，由其串行投影可以知道其工作量满足\(T_1(n)=4T_1(n/2)+\Theta(1)\)，可以得到\(T_1(n)=\Theta(n^2)\)。

这个算法的持续时间满足\(T_{\infty}(n)=3T_\infty(n/2)+\Theta(1)\)，可以得到\(T_{\infty}(n)=\Theta(n^{\lg 3})\)。

最终可以得到并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^{2-\lg 3})\)。

b

不失一般性，这里假设矩阵\(A\)的大小\(n\)是\(3\)的幂次，那么按照题目26-5-a的结果，我们可以对SIMPLE-STENCIL修改成SIMPLE-STENCIL3。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
SIMPLE-STENCIL3(A, n)
  if n == 1
    generate the value of a11
    return
  partition A into n/3 × n/3 submatrices A11, A12, A13, A21, A22, A23, A31, A32, A33; respectively
  SIMPLE-STENCIL3(A11, n / 3)
  spawn SIMPLE-STENCIL3(A12, n / 3)
  spawn SIMPLE-STENCIL3(A21, n / 3)
  sync
  spawn SIMPLE-STENCIL3(A13, n / 3)
  spawn SIMPLE-STENCIL3(A22, n / 3)
  spawn SIMPLE-STENCIL3(A31, n / 3)
  sync
  spawn SIMPLE-STENCIL3(A23, n / 3)
  spawn SIMPLE-STENCIL3(A32, n / 3)
  sync
  SIMPLE-STENCIL3(A33, n / 3)

和题目26-3-a分析的过程类似，可以知道其工作量满足\(T_1(n)=9T_1(n/3)+\Theta(1)\)，可以得到\(T_1(n)=\Theta(n^2)\)。

这个算法的持续时间满足\(T_{\infty}(n)=5T_\infty(n/3)+\Theta(1)\)，可以得到\(T_{\infty}(n)=\Theta(n^{\log_3 5})\)。

最终可以得到并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^{2-\log_3 5})\)。

c

利用类似的方式修改题目26-5-a和26-5-b的代码，可以给出SIMPLE-STENCIL-B代码为：

1
2
3
4
5
6
7
8
9
10
11
SIMPLE-STENCIL-B(A, n)
  if n == 1
    generate the value of a11
    return
  partition A into n/b × n/b submatrices A11, A12, ..., A1b, A21, A22, ..., A2b, ...,  Ab1, Ab2, ..., Abb; respectively
  for k = 2 to b * b
    //注意这里的parallel for 循环将是以硬编码的形式构造出来，因此后面的分析不会对Θ(lg n)这个因子进行考虑。
    parallel for i = 1 to b
      j = k - i
      if 1 <= j and j <= b
      SIMPLE-STENCIL3(Aij, n / b)

和题目26-3-a和26-3-b分析的过程类似，可以知道其工作量满足\(T_1(n)=b^2T_1(n/b)+\Theta(1)\)，可以得到\(T_1(n)=\Theta(n^2)\)。

这个算法的持续时间满足\(T_{\infty}(n)=(2b-1)T_\infty(n/b)+\Theta(1)\)，可以得到\(T_{\infty}(n)=\Theta(n^{\log_b (2b-1)})\)。

最终可以得到并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^{2-\log_b (2b-1)})\)。

为了证明\(T_1(n)/T_{\infty}(n)=\Theta(n^{2-\log_b (2b-1)})=o(n)\)，那么需要证明\(f\forall b\ge 2\)，都有\(2-\log_b(2b-1)<1\)，即证明\(\log_b(2b-1)>1\)。由于\(\forall b\ge 2\)，都有\(2b-1>b\)，因此\(\log_b(2b-1)>1\)成立是很显然的。因此有\(T_1(n)/T_{\infty}(n)=o(n)\)。

d

这个算法由STENCIL给出，并且它还能对更一般形式的矩阵进行填充。

1
2
3
4
5
6
7
STENCIL(A, n)
  for k = 2 to n + n - 1
    up = max{1, k - n}
    down = min{n, k - 1}
    parallel for i = up to down
      j = k - i
      generate the value of aij

可以知道STENCIL的工作量是\(T_1(n)=\Theta(n^2)\)，因为它仅仅是从小到大枚举矩阵的每条副对角线，并按顺序填入一个个元素。

可以知道STENCIL的持续时间满足\(T_\infty(n)=n\cdot\Theta(\lg n)=\Theta(n\lg n)\)。因为第\(k\)轮迭代最多也只会进行\(n\)次操作。如果去掉关键字parallel关键字并转化为spawn ... sync结构，那么这一部分需要花费\(\Theta(\lg n)\)的时间。

因此，STENCIL的并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n/\lg n)\)。

如果需要达到\(\Theta(n)\)的并行度，那么可以将parallel for循环转化成硬编码形式， 从而去掉这个parallel for循环，并且一次性在\(\Theta(1)\)的时间内生产\(\Theta(n)\)个元素，从而使并行度提升到\(\Theta(n)\)。但是这里使用了parallel for循环，它是基于分治法实现的，因此达不到这个最大并行度。

26-6

a

工作量定律将改写成\(E[T_P]\ge E[T_1]/P\)；持续时间定律将改写成\(E[T_P]\ge E[T_{\infty}]\)；贪心调度界限改写成\(E[T_P]\le E[T_1]/P+E[T_{\infty}]\)。

b

按照题目给出的数据，令\(P=10000\)，那么可以写出以下三条式子：

\(\begin{aligned} E[T_1]&=10^4\cdot 0.01+10^9\cdot 0.99\\ E[T_P]&=1\cdot 0.01+10^9\cdot0.99\\ E[T_1/T_P]&=\dfrac{10^4}{1}\cdot0.01+\dfrac{10^9}{10^9}\cdot0.99\\ \end{aligned}\)

从而得到\(E[T_1/T_P]\approx100,\dfrac{E[T_1]}{E[T_P]}\approx 1\)。

选用\(\dfrac{E[T_1]}{E[T_P]}\)的原因有如下两个：

  1. 由于在绝大多数时间中，无论是\(P=1\)还是\(P=10000\)，其运行时间大多数都在\(10^9\)，因此说明这个添加到\(P=10000\)的操作应该对加速作用不大，因此选用\(\dfrac{E[T_1]}{E[T_P]}\)是一个比较恰当的值。

  2. \(T_1/T_P\)这个随机变量并不成立。题目中没有提到\(T_1=10^4\)和\(T_1=10^9\)的时机是否和\(T_P=1\)和\(T_P=10^9\)的时机相同。也就是说，它们不一定是相关的。因此采用独立计算更加合适，即选用\(\dfrac{E[T_1]}{E[T_P]}\)。

c

当\(P\)趋于无穷时，期望的加速比应该和并行度相等。这和题目26-6-b使用的定义是一致的。

d

我们可以对第7.3章提到的RANDOMIZED-QUICKSORT算法提出其并行化版本P-RANDOMIZED-QUICKSORT（按照题目要求，不对RANDOMIZED-PARTITION并行化）：

1
2
3
4
5
6
P-RANDOMIZED-QUICKSORT(A, p, r)
  if p < r:
    q = RANDOMIZED-PARTITION(A, p, r)
    spawn P-RANDOMIZED-QUICKSORT(A, p, q - 1)
    spawn P-RANDOMIZED-QUICKSORT(A, q + 1, r)
    sync

e

可见，P-RANDOMIZED-QUICKSORT的串行投影为RANDOMIZED-QUICKSORT，因此其期望工作量\(E[T_1(n)]=O(n\lg n)\)。

对于P-RANDOMIZED-QUICKSORT持续时间，考虑将P-RANDOMIZED-QUICKSORT和RANDOMIZED-SELECT的行为进行对比。P-RANDOMIZED-QUICKSORT的RANDOMIZED-PARTITION和RANDOMIZED-SELECT中的一样，这部分都是占据了主导的时间\(\Theta(n)\)。因此，对其划分的推导也和第9.2章给的完全一样。接下来考虑其递归部分。可以发现，P-RANDOMIZED-QUICKSORT的阶段划分和RANDOMIZED-SELECT也相同，并且证明过程同样也考虑了那个执行时间比较长的执行分支。因此对P-RANDOMIZED-QUICKSORT持续时间\(T_{\infty}(n)\)的分析过程和RANDOMIZED-SELECT的分析过程完全一致。按照定理9.2的结论，我们得到\(T_{\infty}(n)=\Theta(n)\)。

因此，P-RANDOMIZED-QUICKSORT的并行度为\(T_1(n)/T_{\infty}(n)=O(\lg n)\)。

f

RANDOMIZED-SELECT的并行化版本由P-RANDOMIZED-SELECT给出。其中，给定的P-RANDOMIZED-PARTITION由题目26.3-3的P-PARTITION实现而来。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
P-RANDOMIZED-PARTITION(A, p, r)
  i = RANDOM(p, r)
  exchange A[r] with A[i]
  return P-PARTITION(A, p, r)

P-RANDOMIZED-SELECT(A, p, r, i)
  if p == r
    return A[p]
  q = P-RANDOMIZED-PARTITION(A, p, r)
  k = q – p + 1
  if i == k
    return A[q]
  else if i < k
    return P-RANDOMIZED-SELECT(A, p, q – 1, i)
  else 
    return P-RANDOMIZED-SELECT(A, q + 1, r, i – k)

可见这个算法的串行投影为RANDOMIZED-SELECT，因此它的工作量为\(T_1(n)=\Theta(n)\)。

令示性遍历\(X_k\)表示P-RANDOMIZED-PARTITION划分出来后的元素在于位置\(k\)。那么由于程序的其余部分都需要\(\Theta(\lg n)\)完成（根据题目26.3-3的结论，这里的主要开销就在于P-RANDOMIZED-PARTITION需要花费\(\Theta(\lg n)\)的时间），因此可以对随机变量\(T_{\infty}(n)\)可以写出如下递推式：

\[T_{\infty}(n)=\sum_{i=1}^n X_k\cdot T_{\infty}(\max\{k-1,n-k\})+\Theta(\lg n)\]

可见\(E[X_k]=\dfrac{1}{n}\)，因为每个位置都有可能被选到。那么对式子左右两侧取期望值，那么有

\(\begin{aligned} E[T_{\infty}(n)]&=\sum_{k=1}^n E[X_k]\cdot E[T_{\infty}(\max\{k-1,n-k\})]+\Theta(\lg n)\\ &\le\dfrac{2}{n}\sum_{k=\lfloor n/2\rfloor}^{n-1} E[T_{\infty}(k)]+\Theta(\lg n)\\ \end{aligned}\)

接下来使用代入法证明\(E[T_{\infty}(n)]=O(n^d)\)，即\(\exists c,n_0>0,d\in(0,1)\)，使得\(\forall n\ge n_0\)，都有\(E[T_{\infty}(n)]\le c\cdot n^{d}\)。那么就可以得到

\(\begin{aligned} E[T_{\infty}(n)]&\le\dfrac{2}{n}\sum_{k=\lfloor n/2\rfloor}^{n-1} E[T_{\infty}(k)]+\Theta(\lg n)\\ &\le\dfrac{2}{n}\sum_{k=\lfloor n/2\rfloor}^{n-1} c\cdot k^d+\Theta(\lg n)\\ &=\dfrac{2c}{n}\sum_{k=\lfloor n/2\rfloor}^{n-1} k^d+\Theta(\lg n)\\ &\le \dfrac{2c}{n}\int_{\lfloor n/2\rfloor}^n x^d dx+\Theta(\lg n)\\ &=\dfrac{2c}{n}\cdot\left.\dfrac{x^{d+1}}{d+1}\right|_{x=\lfloor n/2\rfloor}^n+\Theta(\lg n)\\ &=c\cdot\dfrac{2-2^{-d}}{d+1}\cdot n^d + \Theta(\lg n) \end{aligned}\)

考虑关于\(d\)的一元一次不等式\(\dfrac{2-2^{-d}}{d+1}<1\)，可以得到\(d>0\)。

也就是说，无论\(d\)取\((0,1)\)中的什么值，只要第一个项中的\(c\)足够大，它就可以覆盖到\(\Theta(\lg n)\)中的常数，从而最终得到\(T_{\infty}(n)\)

因此，\(T_{\infty}(n)=o(n^d)\)，其中\(d\)是任意正数。

最终可以得到并行度为\(T_1(n)/T_{\infty}(n)=\Theta(n^{1-d})\)。