Параллельная реализация исключения Гаусса с помощью pthreads

Question

Мне нужно реализовать параллельную версию следующего алгоритма исключения Гаусса с использованием pthreads.

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 

Я понимаю, последовательное внедрение и Pthreads, но не получает ни единого намека на то, как построить логику для параллельной версии этого (критерии распределения работы Thread, который петлями распараллеливание и т.д.). Любая подсказка или отправная точка помогут мне продолжить.

Answer 1

Работа для каждой строки в матрице требует, чтобы все предыдущие строки были закончены, поэтому вы не можете разделить эту работу таким образом.

Однако в пределах одной строки каждый столбец может обрабатываться параллельно (с оговоркой, что исходное значение k-й колонки должно быть сохранено и использовано при расчете других столбцов). Это соответствует вашим значениям 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 для подмножества значений j от 0 до n - 1. Этим потокам будет нужен барьер синхронизации после обработки каждой строки, потому что результаты всех столбцов для предыдущей строки необходимы для обработки следующей строки. Вы можете использовать для этого pthread_barrier_wait().

Отмечено, что не каждый столбец (значение j) имеет равную работу. Столбец j обрабатывается этим циклом j раз (поэтому столбец 0 выполняет 0 раз, а столбец n - 1 выполняет n - 1 раз). Вы хотите присвоить номера столбцов потокам, чтобы каждый поток максимально приближался к одинаковой работе для каждой строки - это можно сделать, назначив столбцы циклически. Например. если у вас есть три темы, A, B и C и 10 столбцов от 0 до 9, вы бы назначить их:

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