Реализация алгоритма Фокса

Ниже приведена полная программ, реализующая алгоритм Фокса для перемножения матриц. Сначала приведена функция, которая создает различные коммуникаторы и ассоциированную с ними информацию. Так как ее реализация требует большого количества переменных, которые можно использовать в других функциях, их целесообразно поместить в структуру, чтобы облегчить передачу среди различных функций:

typedef struct {

  int p; /* Общее число процессов */

  MPI_Comm comm; /* Коммуникатор для сетки */

  MPI_Comm row_comm; /* Коммуникатор строки */

  MPI_Comm col_comm; /* Коммуникатор столбца */

  int q; /* Порядок сетки */

  int my_row; /* Номер строки */

  int my_col; /* Номер столбца */

  int my_rank; /* Ранг процесса в коммуникаторе сетки */

} GRID-INFO-TYPE;

 

/* Пространство для сетки выделено

   в вызывающей процедуре.*/

void Setup_grid(GRID_INFO_TYPE* grid) {

  int old_rank;

  int dimensions[2];

  int periods[2];

  int coordinates[2];

  int varying-coords[2];

 

  /* Настройка глобальной информации о сетке */

  MPI_Comm_size(MPI_COMM_WORLD, &(grid->p));

  MPI_Comm_rank(MPI_COMM_WORLD, &old_rank);

  grid->q = (int) sqrt((double) grid->p);

  dimensions[0] = dimensions[1] = grid->q;

  periods[0] = periods[1] = 1;

  MPI_Cart_create(MPI_COMM_WORLD, 2, dimensions,

      periods, 1, &(grid->comm));

  MPI_Comm_rank(grid->comm, &(grid->my_rank));

  MPI_Cart_coords(grid->comm, grid->my_rank, 2,

     coordinates);

  grid->my_row = coordinates[0];

  grid->my_col = coordinates[1];

 

  /* Настройка коммуникаторов для строк и столбцов */

  varying_coords[0] = 0; varying_coords[1] = 1;

  MPI_Cart_sub(grid->comm, varying_coords,

     &(grid->row_comm));

  varying_coords[0] = 1; varying_coords[1] = 0;

  MPI_Cart_sub(grid->comm, varying_coords,

     &(grid->col_comm));

} /* Setup_grid */

Следующая функция выполняет фактическое умножение. Пусть пользователь сам создает определения типа и функции для локальных матриц. При этом определение типа находится в
LOCAL_MATRIX_TYPE, а соответствующий производный тип в
DERIVED_LOCAL_MATRIX. Существуют также три функции:
Local_matrix_multiply, Local_matrix_allocate, Set_to_zero.
Можно также предположить, что память для параметров была выделена в вызывающей функции, и все параметры, кроме локальной матрицы произведения local_C, уже инициализированы:

void Fox(int n, GRID_INFO_TYPE* grid,

     LOCAL-MATRIX-TYPE* local-A,

     LOCAL_MATRIX_TYPE* local_B,

     LOCAL_MATRIX_TYPE* local_C)

{

  LOCAL_MATRIX_TYPE* temp_A;

  int step;

  int bcast_root;

  int n_bar; /* порядок подматрицы = n/q */

  int source;

  int dest;

  int tag = 43;

  MPI_Status status;

 

  n_bar = n/grid->q;

  Set_to_zero(local_C);

 

  /* Вычисление адресов для циклического сдвига B */

  source = (grid->my_row + 1) % grid->q;

  dest = (grid->my_row + grid->q-1) % grid->q;

 

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

  temp_A = Local_matrix_allocate(n_bar);

  for (step = 0; step < grid->q; step++) {

    bcast_root = (grid->my_row + step) % grid->q;

    if (bcast_root == grid->my_col) {

      MPI_Bcast(local_A, 1, DERIVED_LOCAL_MATRIX,

         bcast_root, grid->row_comm);

      Local_matrix_multiply(local_A, local_B, local_C);

    } else {

      MPI_Bcast(temp_A, 1, DERIVED_LOCAL_MATRIX,

         bcast_root, grid->row_comm);

      Local_matrix_multiply(temp_A, local_B, local_C);

  }

  MPI_Send(local_B, 1, DERIVED_LOCAL_MATRIX, dest, tag,

     grid->col_comm);

  MPI_Recv(local_B, 1, DERIVED_LOCAL_MATRIX, source,

     tag, grid->col_comm, &status);

  } /*for*/

}/*Fox*/