3.6 Примеры

Умножение матрицы на вектор

Алгоритм параллельного умножения матрицы на вектор состоит в разбиении исходной матрицы на полосы и нахождении для каждого процесса части результирующего вектора.

Разбиение матрицы происходит с помощью функции scatter, для удобства матрица преобразуется в вектор и его части рассылаются по процессам. Результат собирается с помощью функии gather.

#include <stdio.h> 
#include <stdlib.h> 
#include "mpi.h" 
 
int main (int argc, char* argv[]) 
{ 
    // Умножение матрицы на вектор 
    int procs_rank, procs_count; 
    int i, j, n = 1000, local_n; 
    double *local_a, *b, *local_c, *a, *c; 
    MPI_Init(&argc, &argv); 
    MPI_Comm_size(MPI_COMM_WORLD, &procs_count); 
    MPI_Comm_rank(MPI_COMM_WORLD, &procs_rank); 
    local_n = n / procs_count; 
    local_a = (double *) malloc((local_n * n) * sizeof(double)); 
    b       = (double *) malloc(n * sizeof(double)); 
    local_c = (double *) malloc(local_n * sizeof(double)); 
    c       = (double *) malloc(n * sizeof(double)); 
    a       = (double *) malloc((n * n) * sizeof(double)); 
    //инициализируем матрицу в нулевом процессе 
    // и рассылаем его части по процессам 
    if (procs_rank == 0) 
    { 
    //матрицу удобнее рассматривать как одномерный массив но с индексами i*n+j 
        for (i=0; i<n; i++) 
            for (j=0; j<n; j++) 
                a[i*n+j]=rand(); 
        for (i=0; i<n; i++) 
            b[i] = rand(); 
    } 
    //рассылаем вектор b 
    MPI_Bcast(b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
    //разделяем матрицу на горизонтальные полосы 
    //шириной local_a и отпраляем процессам 
    MPI_Scatter(a, n * local_n, MPI_DOUBLE, local_a, n * local_n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
 
    for (i = 0; i < local_n; i++) 
        for (j = 0; j < n; j++) 
        local_c[i] += local_a[i*n+j] * b[j]; 
    //собираем результат в нулевом процессе 
    MPI_Gather(local_c, local_n, MPI_DOUBLE, c, local_n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 
    if (procs_rank == 0) 
    { 
        for (i=0; i<n; i++) 
            printf("%f3.3 \n",c[i]); 
    } 
    MPI_Finalize(); 
    return 0; 
}

Вычисление определенного интеграла

Широко известно [13], [15], что π
4 = arctan(1) - arctan(0) = 01-dx2
1+x. Заменяя вычисление интеграла конечным суммированием, имеем 01-dx-
1+x2 1
n j = 1j = n--1-
1+x2j, где xj = (j - 0,5)∕n, n – число участков суммирования при численном интегрировании.

Площадь каждого участка (вертикальной “полоски” – stripe) вычисляется функцией COMPUTE_INTERVAL как произведение ширины “полоски” (width) на значение функции в центре “полоски”, далее площади суммируются главным процессом (используется равномерная сетка).

С целью уяснения принципов распределения вычислений рекомендуется проанализировать текст COMPUTE_INTERVAL (здесь j – номер участка интегрирования, myrank – номер данного вычислительного узла, intervals – общее число интервалов численного интегрирования, ntasks – общее число вычислительных узлов, значение локальных сумм накапливается в localsum):

localsum = 0.0; 
for (j = myrank; j < intervals; j += ntasks) { 
    x = (j + 0.5) * width; 
    localsum += 4 / (1 + x * x); 
}

Заметим, что нагрузка каждого процесса вычислением конкретного значения столь простой подынтегральной функции было бы крайне нерациональным решением, т.к. время обмена сообщениями (посылка главным процессом значения параметра функции и прием вычисленного значения) сравнимо со временем вычисления функции и существенно повысило бы время выполнения программы (чрезмерно много коротких пересылок); подобный подход является примером излишне тонкозернистого (fine-grained) параллелизма. В целом следует стремиться к использованию возможно меньшего количества максимально длинномерных обменов.

В качестве функции времени совместно с MPI_Wtime применяется стандартная С-функция локального для вычислительного узла времени ftime (функция F_TIME). В программе использованы исключительно (блокирующие) функции обмена “точка-точка” MPI_Send/MPI_Recv, активно используется барьерная функция MPI_Barrier (при вызове MPI_Barrier(MCW) управление в вызывающую программу не возвращается до тех пор, пока все процессы группы не вызовут MPI_Barrier(MCW), где MCW – коммуникатор). Синхронизация с помощью барьеров используется, например, для завершения всеми процессами некоторого этапа решения задачи, результаты которого будут использоваться на следующем этапе. Использование барьера гарантирует, что ни один из процессов не приступит раньше времени к выполнению следующего этапа, пока результат работы предыдущего не будет окончательно сформирован.

#include "mpi.h" 
#include <stdio.h> 
#include <math.h> 
#include <sys/timeb.h> 
double f_time(void); 
#include "f_time.c" 
double compute_interval (int myrank, int ntasks, long intervals); 
int main(int argc, char ** argv) 
{ 
  long intervals=100;

Подключаем модули, определяем функцию времени и задаем количество разбиений для интегрирования;

  int i, myrank,ranksize; 
  double pi, di, send[2], recv[2], 
  t1,t2, t3,t4, t5,t6,t7,

Заводим служебные переменные и указываем интервалы времени;

  pi_prec=4.0*(atanl(1.0)-atan(0.0));

точное решение числа π находим вышеуказанным методом;

  MPI_Status status; 
 
  t1=f_time(); 
  MPI_Init (&argc, &argv); 
  t2=f_time(); 
  MPI_Comm_rank (MPI_COMM_WORLD, &myrank); 
  MPI_Comm_size (MPI_COMM_WORLD, &ranksize);

Определяем ранг и место под MPI систему;

  if (myrank == 0) 
    { 
    intervals = atoi( argv[1] ); 
    printf ("Calculation of PI by numerical integration with %ld intervals\n", intervals); 
    } 
  MPI_Barrier(MPI_COMM_WORLD);

Используем барьер для удостоверения того, что MPI процессы запущены;

  if (myrank == 0) 
    { 
    printf ("Master: Sending # of intervals to MPI-Processes \n"); 
    t3 = MPI_Wtime(); 
    for (i=1; i<ranksize; i++) 
    MPI_Send (&intervals, 1, MPI_LONG, i, 98, MPI_COMM_WORLD); 
    } 
  else 
    { 
    MPI_Recv (&intervals, 1, MPI_LONG, 0, 98, MPI_COMM_WORLD, &status); 
    } 
 
  t4 = MPI_Wtime(); 
  pi = compute_interval (myrank, ranksize, intervals);

Вычисляем число π;

  t5 = MPI_Wtime(); 
  MPI_Barrier (MPI_COMM_WORLD); 
  t6 = MPI_Wtime(); 
  if (myrank == 0) 
  {

После выполнения операций, собираем на ведущем процессоре результаты;

  for (i=1; i<ranksize; i++) 
    { 
    MPI_Recv (&di, 1, MPI_DOUBLE, i, 99, MPI_COMM_WORLD, &status); 
    pi += di; 
    } 
  t7 = MPI_Wtime(); 
  printf ("Master: Has collected sum from MPI-Processes \n"); 
  printf ("\nPi estimation: %.12lf (rel.error= %.5f %%)\n", pi, 1.0e2*(pi_prec-pi)/pi_prec);

оценка числа π;

  printf("%ld tasks used, execution time: %.3lf sec \n",ranksize, t7-t3); 
  printf("\nStatistics:\n"); 
  printf("Master: startup: %.0lf msec\n",t2-t1); 
  printf("Master: time to send # of intervals:%.3lf sec\n",t4-t3); 
  printf("Master: waiting time for sincro after calculation:%.2lf sec\n",t6-t5); 
  printf("Master: time to collect: %.3lf sec\n",t7-t6); 
  printf("Master: calculation time:%.3lf sec\n",t5-t4); 
 
  MPI_Barrier (MPI_COMM_WORLD); 
  for (i=1; i<ranksize; i++) 
    { 
    MPI_Recv(recv, 2, MPI_DOUBLE, i, 100, MPI_COMM_WORLD, &status); 
    printf("process %d: calculation time: %.3lf sec,\twaiting time for sincro.: %.3lf sec\n", i,recv[0],recv[1]); 
    } 
  } 
  else 
  { 
  MPI_Send (&pi, 1, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD); 
  MPI_Barrier (MPI_COMM_WORLD); 
  send[0]=t5-t4; 
  send[1]=t6-t5; 
  MPI_Send (send, 2, MPI_DOUBLE, 0, 100, MPI_COMM_WORLD); 
  } // end work SLAVE 
  MPI_Finalize (); 
  } // end MAIN function

Завершаем основную программу и финализируем MPI;

double compute_interval (int myrank, int ntasks, long intervals) 
{ 
  double x, width, localsum=0.0; 
  long j; 
  width = 1.0 / intervals; 
  for (j=myrank; j<intervals; j+= ntasks) 
    { 
    x = (j + 0.5) * width; 
    localsum += 4 / (1 + x*x); 
    } 
  return (localsum * width); 
  } // end of COMPUTE_INTERVAL function 
  // end of PI_01.C program

Функция вычисления частичной суммы по интервалу;

Функция f_time содержится в файле F_TIME.C (подключается посредством #include "f_time.c"):

double f_time() 
{ 
struct timeb tp; 
  ftime(&tp); 
  return ((double)(tp.time)+1.0e-3*(double)(tp.millitm)); 
}

Возможно точное (double при данной разрядной сетке ЭВМ) значение π определяется как 4.0(atanl(1.0) - atanl(0.0)) и сравнивается с численно вычисленным (относительная ошибка rel.error выводится в процентах).