Распараллеливать этот алгоритм, чтобы ускорить его работу

Question

Как вызов, меня попросили сделать параллельный алгоритм для инвертирования матрицы. В основном я смотрел на this paper и this SO question, проводя исследования для этого.

Прежде чем я попытался написать свой собственный код, я наткнулся на someone else's implementation.

Я исхожу из объектива-c фона, поэтому я сразу же подумал об использовании GCD для этой задачи. Я также наткнулся на то, что называется POSIX, которое выглядит более низкоуровневым и может быть уместным для этой задачи, если GCD не будет работать - я не знаю.

Моя наивная попытка распараллеливать это была только для замены каждого цикла на dispatch_apply, который работал (произведение исходного и обратного выражает идентичную матрицу). Тем не менее, это только замедлило ситуацию значительно (около 20 раз, как медленно с первого взгляда). Я вижу, что есть SO questions on GCD and for-loops, но меня в основном интересует, какой может быть лучший подход к этому, а не ссылки на те ответы, которые я уже прочитал. Возможно ли, что проблема заключается в том, как я создаю очередь отправки, или, может быть, в том, что я использую только одну очередь отправки?

#include <stdio.h> 
#include <dispatch/dispatch.h> 
#include <stdlib.h> 
#include <math.h> 
#include <time.h> 

#define PARALLEL true 

void invertMatrixNonParallel(double **matrix, long n); 
void invertMatrixParallel(double **matrix, long n, dispatch_queue_t q); 

void invertMatrixParallel(double **matrix, long n, dispatch_queue_t q) 
{ 
    __block double r; 
    __block long temp; 

    dispatch_apply(n, q, ^(size_t i) { 
     dispatch_apply(n, q, ^(size_t j) { 
      matrix[i][j + n] = (j == i) ? 1 : 0; 
     }); 
    }); 
    /* using gauss-jordan elimination */ 

    dispatch_apply(n, q, ^(size_t j) { 
     temp=j; 

     /* finding maximum jth column element in last (n-j) rows */ 

     dispatch_apply(n - j - 1, q, ^(size_t i) { 
      if (matrix[i + j + 1][j] > matrix[temp][j]) 
      { 
       temp = i + j + 1; 
      } 
     }); 

     /* swapping row which has maximum jth column element */ 

     if(temp!=j) 
     { 
      double *row = matrix[j]; 
      matrix[j] = matrix[temp]; 
      matrix[temp] = row; 
     } 

     /* performing row operations to form required identity matrix out of the input matrix */ 
     dispatch_apply(n, q, ^(size_t i) { 
      r = matrix[i][j]; 

      if (i == j) 
      { 
       dispatch_apply(2 * n, q, ^(size_t k) { 
        matrix[i][k]/=r ; 
       }); 
      } 
      else 
      { 
       dispatch_apply(2 * n, q, ^(size_t k) { 
        matrix[i][k]-=(matrix[j][k]/matrix[j][j])*r ; 
       }); 
      } 
     }); 
    }); 
} 

void invertMatrixNonParallel(double **matrix, long n) 
{ 
    double temporary, r; 
    long i, j, k, temp; 

    for (i = 0; i < n; ++i) 
    { 
     for (j = n; j < n * 2; ++j) 
     { 
      matrix[i][j] = (j == i + n) ? 1 : 0; 
     } 
    } 
    /* using gauss-jordan elimination */ 

    for(j=0; j<n; j++) 
    { 
     temp=j; 

     /* finding maximum jth column element in last (n-j) rows */ 

     for(i=j+1; i<n; i++) 
      if(matrix[i][j]>matrix[temp][j]) 
       temp=i; 

     /* swapping row which has maximum jth column element */ 

     if(temp!=j) 
     { 
      for(k=0; k<2*n; k++) 
      { 
       temporary=matrix[j][k] ; 
       matrix[j][k]=matrix[temp][k] ; 
       matrix[temp][k]=temporary ; 
      } 
     } 

     /* performing row operations to form required identity matrix out of the input matrix */ 

     for(i=0; i<n; i++) 
     { 
      if(i!=j) 
      { 
       r=matrix[i][j]; 
       for(k=0; k<2*n; k++) 
        matrix[i][k]-=(matrix[j][k]/matrix[j][j])*r ; 
      } 
      else 
      { 
       r=matrix[i][j]; 
       for(k=0; k<2*n; k++) 
        matrix[i][k]/=r ; 
      } 
     } 
    } 
} 

#pragma mark - Main 

int main(int argc, const char * argv[]) 
{ 
    long i, j, k; 
    const long n = 5; 
    const double range = 10.0; 
    __block double **matrix; 
    __block double **invertedMatrix = malloc(sizeof(double *) * n); 

    matrix = malloc(sizeof(double *) * n); 
    invertedMatrix = malloc(sizeof(double *) * n); 
    for (i = 0; i < n; ++i) 
    { 
     matrix[i] = malloc(sizeof(double) * n); 
     invertedMatrix[i] = malloc(sizeof(double) * n * 2); 
     for (j = 0; j < n; ++j) 
     { 
      matrix[i][j] = drand48() * range; 
      invertedMatrix[i][j] = matrix[i][j]; 
     } 
    } 

    clock_t t; 

#if PARALLEL 
    dispatch_queue_t q1 = dispatch_queue_create("com.example.queue1", DISPATCH_QUEUE_CONCURRENT); 
    t = clock(); 
    invertMatrixParallel(invertedMatrix, n, q1); 
#else 
    t = clock(); 
    invertMatrixNonParallel(invertedMatrix, n); 
#endif 

    t = clock() - t; 
    double time_taken = ((double)t * 1000)/CLOCKS_PER_SEC; // in seconds 

    printf("\n%s took %f milliseconds to execute \n\n", (PARALLEL == true) ? "Parallel" : "Non-Parallel", time_taken); 

    printf("Here's the product of the inverted matrix and the original matrix\n"); 
    double product[n][n]; 
    for (i = 0; i < n; ++i) 
    { 
     for (j = 0; j < n; ++j) 
     { 
      double sum = 0; 
      for (k = 0; k < n; ++k) 
      { 
       sum += matrix[i][k] * invertedMatrix[k][j + n]; 
      } 
      product[i][j] = sum; 
     } 
    } 

    // should print the identity matrix 
    for (i = 0; i < n; ++i) 
    { 
     for (j = 0; j < n; ++j) 
     { 
      printf("%5.2f%s", product[i][j], (j < n - 1) ? ", " : "\n"); 
     } 
    } 

    return 0; 
} 

Выход для параллельного:

Parallel took 0.098000 milliseconds to execute 

Для непараллельных:

Non-Parallel took 0.004000 milliseconds to execute 

Для обоих:

Here's the product of the inverted matrix and the original matrix 
1.00, -0.00, -0.00, 0.00, -0.00 
0.00, 1.00, 0.00, 0.00, 0.00 
0.00, -0.00, 1.00, -0.00, 0.00 
-0.00, -0.00, -0.00, 1.00, 0.00 
0.00, 0.00, 0.00, 0.00, 1.00 

Пожалуйста, нет ответов, которые только ссылки, я не буду используя только SO как последнее средство.

Answer 1

0) Как уже упоминалось в комментариях, нужна большая матрица. Создание параллельного потока занимает некоторое накладное время, поэтому вы не можете сделать более быструю параллельную версию, если она занимает слишком мало времени. Даже если вам удастся добиться лучшей производительности для малой матрицы, трудно точно измерить ее.

1)

dispatch_apply(n, q, ^(size_t i) { 
    dispatch_apply(n, q, ^(size_t j) { 
     matrix[i][j + n] = (j == i) ? 1 : 0; 
    }); 
}); 

Существует не так много смысла распараллеливания каждого вложенного цикла. Нет смысла добавлять каждую отдельную операцию по очереди в очереди отправки, поскольку она по-прежнему требует некоторых накладных расходов, поэтому лучше добавить некоторые нетривиальные блоки.

dispatch_apply(n, q, ^(size_t i) { 
    for (j = n; j < n * 2; ++j) { 
     matrix[i][j + n] = (j == i) ? 1 : 0; 
    } 
}); 

Достаточно.

2) Вам нужно узнать о безопасности потоков и хорошо понимать ваш алгоритм, или вы можете столкнуться с непредсказуемым и невоспроизводимым неправильным поведением вашего приложения. Я не уверен, что существует много циклов, которые можно сопоставить эффективно и безопасно, за исключением упомянутой выше инициализации, и одну, отмеченную/* выполнением операций строки, чтобы сформировать требуемую матрицу идентификации из входной матрицы. */

Итак, вы вероятно, нужно найти какой-то конкретный алгоритм инверсии параллельной матрицы.