next up previous contents
Next: Одномерное температурное уравнение Up: PVM - параллельная виртуальная Previous: Пример с ошибкой   Contents

Матричное умножение

В следующем примере программируется алгоритм матричного умножения, предложенный Фоксом (Fox)[8]. Программа
mmult будет вычислять $C=AB$, где $C$, $A$ и $B$ - все квадратные матрицы. С целью упрощения предполагается, что для вычисления результата будут использоваться $m\times m$ задач. Каждая задача будет вычислять подблок результирующей матрицы $C$. Размер блока и значение $m$ передаются программе как аргументы командной строки. Матрицы $A$ и $B$ также устанавливаются в виде блоков, распределенных в среде из $m^{2}$ задач.

Предположим, что существует ``сеть'' из $m\times m$ задач. Каждая задача ($t_{i,j}$, где $0\leq i,j<m$) первоначально содержит блоки $C_{i,j}$, $A_{i,j}$ и $B_{i,j}$. Первым шагом алгоритма задачи, расположенные по диагонали ($t_{i,j}$, где $i=j$), передают свои блоки $A_{i,i}$ всем остальным задачам в строке $i$. После передачи $A_{i,i}$ каждая задача вычисляет $A_{i,i}\times B_{i,j}$ и добавляет результат к $C_{i,j}$. Следующим шагом блоки $B$ поворачиваются по столбцам. Так, задача $t_{i,j}$ передает свой блок $B$ задаче $t_{i-1,j}$ (задача $t_{0,j}$ передает свой блок $B$ задаче $t_{m-1,j}$). Затем задачи возвращаются к первому шагу: $A_{i,i+1}$ широковещательно передается всем другим задачам в строке $i$ - алгоритм продолжается. После $m$ итераций матрица $C$ будет содержать $A\times B$, а матрица $B$ вернется к исходному состоянию.

В PVM не накладывается ограничений на то, как задача может связываться с любой другой задачей. Однако для данной программы хотелось бы представлять задачи в виде двухмерной модели тора. Чтобы учесть задачи, каждая из них включается в группу mmult. Групповые идентификаторы используются для отображения задач на наш тор. При включении в группу первой задачи она получает групповой идентификатор 0. В программе mmult, задача с нулевым групповым идентификатором порождает остальные задачи и передает параметры матричного умножения этим задачам. Такими параметрами являются m и bklsize - квадратный корень числа блоков и соответственно размер блока. После того, как все задачи порождены, а параметры переданы, вызывается pvm_barrier() - чтобы удостовериться, что все задачи присоединились к группе. Если барьер не организован, то последующий вызов pvm_gettid() может завершиться ненормально, так как некая задача к тому моменту может быть еще не в составе группы.

После того как барьер организован, сохраняются идентификаторы задач из строки myrow в массиве. Это достигается путем определения групповых идентификаторов всех задач в строке и запрашивания PVM о соответствующих им идентификаторах задач. Далее, с помощью malloc() выделяется место для блоков матриц. В действительном программном приложении, может случиться так, что матрицы уже распределены. Далее программа определяет строку и столбец блока $C$, который будет вычисляться. Это основывается на значениях групповых идентификаторов - в диапазоне от $0$ до $m-1$ включительно. Так, если допустить построчное отображение групповых идентификаторов задач, целочисленное деление mygid/m даст строку задач, а mygid mod m - даст столбец. С использованием подобного отображения определяется групповой идентификатор задачи, расположенной в торе ``выше'' или ``ниже'', и сохраняется соответственно в up и down.

Далее блоки инициализируются вызовом InitBlock(): $A$ - случайными числами, $B$ - идентичной матрицей, а $C$ - нулями. Это позволит нам верифицировать вычисление в конце программы простой проверкой равенства $A=C$.

Наконец, выполняется главный вычислительный цикл матричного умножения. Сначала диагональные задачи широковещательно передают имеющиеся блоки $A$ другим задачам в своих строках. Отметим, что массив myrow в действительности содержит идентификаторы задач, выполняющих широковещательные передачи. При повторных вызовах pvm_mcast() сообщения будут передаваться всем задачам из массива, исключая делающую вызов задачу. Для задач mmult эта процедура великолепно срабатывает - чтобы напрасно не обрабатывать ``лишнее'' сообщение, поступающее самой широковещательно передающей задаче при ``лишнем'' pvm_recv(). Как широковещательно передающая задача, так и задачи, принимающие блок, вычисляют $A\times B$ с использованием и диагонального блока $B$ и блоков $B$, принадлежащих задачам.

После того как подблоки умножены и результат добавлен к блоку $C$, вертикально сдвигаются блоки $B$. Особенным образом блок $B$ упаковывается в сообщение, передается задаче с идентификатором задачи up и затем принимается новый блок $B$ от задачи с идентификатором задачи down.

Обратите внимание на то, что используются различные теги сообщений при передачах блоков $A$ и блоков $B$ при различных циклических итерациях. Полностью указываются и идентификаторы задач при выполнении pvm_recv(). Возникает соблазн применять специальные символы в полях pvm_recv(), однако, такая практика может быть опасной. К примеру, если некорректно определить значение для up и указать специальный символ при pvm_recv() вместо down, то это может привести к ``неосознанной'' передаче сообщений не тем задачам. В примере сообщения однозначно направлены, тем самым уменьшая вероятность возникновения ошибок при приеме сообщений от не тех задач, т.е. ошибочных фаз алгоритма.

Как только вычисление завершено, проверяется $A=C$ - только для того, чтобы убедиться в правильности матричного умножения, т.е. корректности значений $C$. Такую проверку не возможно осуществить, например, с помощью подпрограмм из библиотеки матричного умножения.

В вызове pvm_lvgroup() нет необходимости, поскольку PVM выгрузит отработавшие задачи и удалит их из группы. Однако хорошим примером является ``явный'' выход из группы перед вызовом pvm_exit(). Команда консоли PVM reset ``сбросит'' все группы PVM. Команда pvm_gstat выведет на экран статус любой из существующих групп.

/*

Матричное умножение

*/

/* определения и прототипы библиотеки PVM */

#include <pvm3.h>

#include <stdio.h>

/* максимальное число потомков, которые

   будут порождаться этой программой */

#define MAXNTIDS 100

#define MAXNROW 10

/* теги сообщений */

#define ATAG 2

#define BTAG 3

#define DIMTAG 5

void

InitBlock(float *a, float *b, float *c, int blk,

  int row, int col) {

 int len, ind;

 int i, j;

 srand_pvm_mytid__

 len = blk*blk;

 for (ind = 0; ind < len; ind++)

 { a[ind] = (float)(rand()%1000)/100.0; c[ind] = 0.0; }

  for (i = 0; i < blk; i++) {

   for (j = 0; i < blk; j++) {

    if (row == col)

      b[j*blk+i] = (i==j)? 1.0 : 0.0;

    else

      b[j*blk+i] = 0.0;

   }

  }

}

void

BlockMult(float* c, float* a, float* b, int blk) {

 int i,j,k;

 for (i = 0; i < blk; i++)

  for (j = 0; i < blk; j++)

   for (k = 0; k < blk; k++)

    for (j = 0; i < blk; j++)

     c[i*blk+j] += (a[i*blk+k] * b[k*blk+j]);

}

 

int

main(int argc, char* argv[])

{

 /* количество задач для порождения, 3 используются

   по умолчанию */

 int ntask = 2;

 /* код возврата для вызовов PVM */

 int info;

 /* свой идентификатор задачи и групповой

   идентификатор */

 int mytid, mygid;

 /* массив идентификаторов задач-потомков*/

 int child[MAXNTIDS-1];

 int i, m, blksize;

 /* массив идентификаторов задач в своей строке*/

 int myrow[MAXROW];

 float *a, *b, *c, *atmp;

 int row, col, up, down;

 /* поиск своего идентификатора задачи */

 mytid = pvm_mytid();

 pvm_setopt(PvmRoute, PvmRouteDirect);

 /* проверка на ошибки*/

 if (mytid < 0) {

  /* вывод на экран сообщения об ошибке*/

  pvm_perror(argv[0]);

  /* выход из программы */

  return -1;

 }

 /* присоединение к группе mmult*/

 mygid = pvm_joingroup("mmult");

 if (mygid < 0) {

  pvm_perror(argv[0]); pvm_exit(); return -1;

 }

 /* если свой групповой идентификатор не равен 0,

   то нужно породить другие задачи */

 if (mygid == 0) {

   /* определение числа задач для порождения */

   if (argc == 3) {

    m = atoi(argv[1]);

    blksize = atoi(argv[2]);

   }

   if (argc < 3) {

    fprintf(stderr, "usage: mmult m blk\n");

    pvm_lvgroup("mmult"); pvm_exit(); return -1;

   }

   /* удостоверение, что ntask - допустимо */

   ntask = m*m;

   if ((ntask < 1) || (ntask > MAXNTIDS)) {

   fprintf(stderr, "ntask = %d not valid.\n",

           ntask);

   pvm_lvgroup("mmult"); pvm_exit(); return -1;

   }

   /* порождать не нужно, если имеется

    только одна задача */

   if (ntask == 1) goto barrier;

   /* порождение задач-потомков*/

   info = pvm_spawn("mmult", (char**)0,

        PvmTaskDefault, (char*)0, ntask-1, child);

   /* удостоверение, что порождение

     произошло успешно */

   if (info != ntask-1) {

   pvm_lvgroup("mmult"); pvm_exit(); return -1;

   }

   /* передача размерности матрицы */

   pvm_initsend(PvmDataDefault);

   pvm_pkint(&m, 1, 1);

   pvm_pkint(&blksize, 1, 1);

   pvm_mcast(child, ntask-1, DIMTAG);

 }

 else {

   /* прием размерности матрицы */

   pvm_recv(pvm_gettid("mmult", 0), DIMTAG);

   pvm_pkint(&m, 1, 1);

   pvm_pkint(&blksize, 1, 1);

   ntask = m*m;

 }

 /* удостоверение, что все задачи

   присоединились к группе */

 barrier:

 info = pvm_barrier("mmult",ntask);

 if (info < 0) pvm_perror(argv[0]);

 /* поиск идентификаторов задач в своей строке */

 for (i = 0; i < m; i++)

 myrow[i] = pvm_gettid("mmult", (mygid/m)*m + i);

 /* распределение памяти для локальных блоков */

 a = (float*)malloc(sizeof(float)*blksize+blksize);

 b = (float*)malloc(sizeof(float)*blksize+blksize);

 c = (float*)malloc(sizeof(float)*blksize+blksize);

 atmp = (float*)malloc(sizeof(float)*blksize+blksize);

 /* проверка достоверности указателей */

 if (!(a && b && c && atmp)) {

   fprintf(stderr, "%s: out of memory!\n", argv[0]);

   free(a); free(b); free(c); free(atmp);

   pvm_lvgroup("mmult"); pvm_exit(); return -1;

 }

 /* поиск строки и столбца своего блока */

 row = mygid/m; col = mygid % m;

 /* определение соседей сверху и снизу */

 up = pvm_gettid("mmult",

    ((row)?(row-1):(m-1))*m+col);

 down = pvm_gettid("mmult",

    ((row) == (m-1))?col:(row+1)*m+col));

 /* инициализация блоков */

 InitBlock(a, b, c, blksize, row, col);

 /* выполнение матричного умножения */

 for (i = 0; i < m; i++) {

 /* широковещательная передача блока матрицы A */

 if (col == (row + i)%m) {

   pvm_initsend(PvmDataDefault);

   pvm_pkfloat(a, blksize*blksize, 1);

   pvm_mcast(myrow, m, (i+1)*ATAG);

   BlockMult(c,a,b,blksize);

 }

 else {

   pvm_recv(pvm_gettid("mmult", row*m + (row +i)%m),

      (i+1)*ATAG);

   pvm_upkfloat(atmp, blksize*blksize);

   BlockMult(c,atmp,b,blksize);

 }

 /* поворот столбца B */

 pvm_initsend(PvmDataDefault);

 pvm_pkfloat(b, blksize*blksize, 1);

 pvm_send(up, (i+1)*BTAG);

 pvm_recv(down, (i+1)*BTAG);

 pvm_upkfloat(b, blksize*blksize, 1);

 }

 /* проверка */

 for (i = 0; i < blksize*blksize; i++)

 if (a[i] != c[i])

 printf("Error a[%d] (%g) != c[%d] (%g) \n", i,

    a[i], i, c[i]);

 printf("Done.\n");

 free(a); free(b); free(c); free(atmp);

 pvm_lvgroup("mmult");

 pvm_exit();

 return 0;

}

next up previous contents
Next: Одномерное температурное уравнение Up: PVM - параллельная виртуальная Previous: Пример с ошибкой   Contents
2004-06-22