用线程并行实现Gauss消去

Question

我需要实现并行版本的以下高斯消除算法使用线程。

procedure GAUSSIAN ELIMINATION (A, b, y)
begin
    for k := 0 to n − 1 do /* Outer loop */
    begin
        for j := k + 1 to n − 1 do
            A[k, j] := A[k, j]/A[k, k]; /* Division step */
        y[k] := b[k]/A[k, k];
        A[k, k] := 1;
        for i := k + 1 to n − 1 do
        begin
            for j := k + 1 to n − 1 do
                A[i, j] := A[i, j] − A[i, k] × A[k, j]; /* Elimination step */
            b[i] := b[i] − A[i, k] × y[k];
            A[i, k] := 0;
        endfor; /* Line 9 */
    endfor; /* Line 3 */
end GAUSSIAN ELIMINATION

我理解顺序实现和线程，但对于如何构造并行版本的逻辑(线程工作分配准则，循环并行化等)，我没有得到任何提示。任何线索或起点都会帮助我继续下去。

Answer 1

矩阵中每一行的工作都需要完成前面的所有行，因此不能这样划分工作。

但是，在单个行中，可以并行处理每一列(请注意，必须保存k-th列的原始值，并将其用于计算其他列)。这对应于您的j值。

我相信您可以重新安排算法以使这更容易，这样在j上只有一个循环

procedure GAUSSIAN ELIMINATION (A, b, y)
begin
    for k := 0 to n − 1 do /* Outer loop */
    begin
        for j := k + 1 to n − 1 do
        begin
            A[k, j] := A[k, j]/A[k, k]; /* Division step */
            for i := k + 1 to n − 1 do
                A[i, j] := A[i, j] − A[i, k] × A[k, j]; /* Elimination step */
        endfor;
        y[k] := b[k]/A[k, k];
        A[k, k] := 1;
        for i := k + 1 to n − 1 do
        begin
            b[i] := b[i] − A[i, k] × y[k];
            A[i, k] := 0;
        endfor;
    endfor;
end GAUSSIAN ELIMINATION

注意，j上的循环主体只访问j和k列中的值--这是可以并行完成的循环。然后您还可以进一步注意，外部循环的第二部分不依赖于第一部分，因此可以将外部循环分成两个部分：

procedure GAUSSIAN ELIMINATION (A, b, y)
begin
    for k := 0 to n − 1 do /* Outer loop */
    begin
        for j := k + 1 to n − 1 do
        begin
            A[k, j] := A[k, j]/A[k, k]; /* Division step */
            for i := k + 1 to n − 1 do
                A[i, j] := A[i, j] − A[i, k] × A[k, j]; /* Elimination step */
        endfor;
    endfor;

    for k := 0 to n − 1 do /* Outer loop, second pass */
    begin
        y[k] := b[k]/A[k, k];
        A[k, k] := 1;
        for i := k + 1 to n − 1 do
        begin
            b[i] := b[i] − A[i, k] × y[k];
            A[i, k] := 0;
        endfor;
    endfor;
end GAUSSIAN ELIMINATION

您可以预先创建线程，每个线程负责在j上执行从0到n - 1的j值子集的循环。在处理每一行后，这些线程将需要一个同步屏障，因为处理下一行需要前一行的所有列的结果。为此您可以使用pthread_barrier_wait()。

您将注意到，并非每一列(j的值)都具有相同的工作。列j由该循环j次数处理(因此列0执行0次，列n - 1执行n - 1时间)。您将希望为线程分配列号，这样每个线程对每一行都有尽可能接近相同的工作--这可以通过以循环方式分配列来完成。例如：如果您有三个线程A、B和C，以及从0到9的10列，您将分配它们：

Thread A: 0, 3, 6, 9
Thread B: 1, 4, 7,
Thread C: 2, 5, 8,

然后，线程函数看起来如下所示：

for k := 0 to n − 1 do /* Outer loop */
begin
    call pthread_barrier_wait(row_barrier);
    for j := k + 1 + thread_number to n − 1 step n_threads do
    begin
        A[k, j] := A[k, j]/A[k, k]; /* Division step */
        for i := k + 1 to n − 1 do
            A[i, j] := A[i, j] − A[i, k] × A[k, j]; /* Elimination step */
    endfor;
endfor;

call pthread_barrier_wait(phase1_barrier);

其主要功能如下：

procedure GAUSSIAN ELIMINATION (A, b, y)
begin
    call start_threads;
    call pthread_barrier_wait(phase1_barrier);

    for k := 0 to n − 1 do /* Outer loop, second pass */
    begin
        y[k] := b[k]/A[k, k];
        A[k, k] := 1;
        for i := k + 1 to n − 1 do
        begin
            b[i] := b[i] − A[i, k] × y[k];
            A[i, k] := 0;
        endfor;
    endfor;
end GAUSSIAN ELIMINATION