Распараллеливание алгоритма Фокса

#include <Mpi.h>
#include <stdio.h>
 
int ProcNum = 0;    // Количество доступных процессов 
int ProcRank = 0;   // Ранг текущего процесса
int GridSize;       // Размер виртуальной решетки процессов
int GridCoords[2];      // Координаты текущего процесса в процессной
// решетке
MPI_Comm GridComm;      // Коммуникатор в виде квадратной решетки
MPI_Comm ColComm;       // коммуникатор – столбец решетки
MPI_Comm RowComm;   
 
// Функция для параллельного умножения матриц
void ParallelResultCalculation(double* pAblock, double* pMatrixAblock, double* pBblock, double* pCblock, int BlockSize) {
    for (int iter = 0; iter < GridSize; iter ++) 
    {
        // Рассылка блоков матрицы A по строкам процессной решетки
        ABlockCommunication (iter, pAblock, pMatrixAblock, BlockSize);
        // Умножение блоков
        BlockMultiplication(pAblock, pBblock, pCblock, BlockSize);
        // Циклический сдвиг блоков матрицы B в столбцах процессной 
        // решетки
        BblockCommunication(pBblock, BlockSize);
    }
}
 
 
// Циклический сдвиг блоков матрицы В вдоль столбца процессной 
// решетки 
void BblockCommunication (double *pBblock, int BlockSize) 
{
    MPI_Status Status;
    int NextProc = GridCoords[0] + 1;
    if ( GridCoords[0] == GridSize-1 ) 
    NextProc = 0;
    
    int PrevProc = GridCoords[0] - 1;
    if ( GridCoords[0] == 0 ) 
        PrevProc = GridSize-1;
        
    MPI_Sendrecv_replace( pBblock, BlockSize*BlockSize, MPI_DOUBLE,NextProc, 0, PrevProc, 0, ColComm, &Status);
}
 
// Умножение матричных блоков
void BlockMultiplication (double *pAblock, double *pBblock, double *pCblock, int BlockSize) 
{
  // Вычисление произведения матричных блоков
  for (int i=0; i<BlockSize; i++) 
    {
        for (int j=0; j<BlockSize; j++) 
        {
            double temp = 0;
            for (int k=0; k<BlockSize; k++ )
                temp += pAblock [i*BlockSize + k] * pBblock [k*BlockSize + j];
            pCblock [i*BlockSize + j] += temp;
        }
    }
}
 
 
// Рассылка блоков матрицы А по строкам решетки процессов 
void ABlockCommunication (int iter, double *pAblock, double* pMatrixAblock, int BlockSize) {
 
  // Определение ведущего процесса в строке процессной решетки 
  int Pivot = (GridCoords[0] + iter) % GridSize;
  
  // Копирование передаваемого блока в отдельный буфер памяти
  if (GridCoords[1] == Pivot) 
  {
    for (int i=0; i<BlockSize*BlockSize; i++)
        pAblock[i] = pMatrixAblock[i];
  }
  
  // Рассылка блока
  MPI_Bcast(pAblock, BlockSize*BlockSize, MPI_DOUBLE, Pivot,RowComm);
}
 
 
// Функция для выделения памяти и инициализации исходных данных
void ProcessInitialization (double* &pAMatrix, double* &pBMatrix, double* &pCMatrix, double* &pAblock, double* &pBblock, double* &pCblock, double* &pTemporaryAblock, int &Size, int &BlockSize ) {
    if (ProcRank == 0) 
    {
        do
        {
            printf("\nВведите размер матриц: ");
            scanf("%d", &Size);
      
            if (Size%GridSize != 0) 
            {
                printf ("Размер матриц должен быть кратен размеру сетки! \n");
            }
        }
        while (Size%GridSize != 0);
      }
    MPI_Bcast(&Size, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
    BlockSize = Size/GridSize;
 
    pAblock = new double [BlockSize*BlockSize];
    pBblock = new double [BlockSize*BlockSize];
    pCblock = new double [BlockSize*BlockSize];
    pTemporaryAblock = new double [BlockSize*BlockSize];
 
    for (int i=0; i<BlockSize*BlockSize; i++) 
    {
        pCblock[i] = 0;
    }
    if (ProcRank == 0) 
    {
        pAMatrix = new double [Size*Size];
        pBMatrix = new double [Size*Size];
        pCMatrix = new double [Size*Size];
        RandomDataInitialization(pAMatrix, pBMatrix, Size);
     } 
}
 
// Создание коммуникатора в виде двумерной квадратной решетки 
// и коммуникаторов для каждой строки и каждого столбца решетки
void CreateGridCommunicators() 
{
      int DimSize[2];   // Количество процессов в каждом измерении 
                        // решетки
      int Periodic[2];  // =1 для каждого измерения, являющегося 
                        // периодическим 
      int Subdims[2];   // =1 для каждого измерения, оставляемого 
                        // в подрешетке
      DimSize[0] = GridSize; 
      DimSize[1] = GridSize;
      Periodic[0] = 0;
      Periodic[1] = 0;
 
      // Создание коммуникатора в виде квадратной решетки 
      MPI_Cart_create(MPI_COMM_WORLD, 2, DimSize, Periodic, 1, &GridComm);
 
      // Определение координат процесса в решетке 
      MPI_Cart_coords(GridComm, ProcRank, 2, GridCoords);
      
      // Создание коммуникаторов для строк процессной решетки
      Subdims[0] = 0;  // Фиксация измерения
      Subdims[1] = 1;  // Наличие данного измерения в подрешетке
      MPI_Cart_sub(GridComm, Subdims, &RowComm);
      
      // Создание коммуникаторов для столбцов процессной решетки
      Subdims[0] = 1;
      Subdims[1] = 0;
      MPI_Cart_sub(GridComm, Subdims, &ColComm);
}
 
 
int main ( int argc, char * argv[] ) 
{
  double* pAMatrix;     // Первый аргумент матричного умножения
  double* pBMatrix;     // Второй аргумент матричного умножения
  double* pCMatrix;     // Результирующая матрица
  int Size;         // Размер матриц
  int BlockSize;    // Размер матричных блоков, расположенных 
                    // на процессах
  double *pAblock;      // Блок матрицы А на процессе
  double *pBblock;      // Блок матрицы В на процессе
  double *pCblock;      // Блок результирующей матрицы С на процессе
  double *pMatrixAblock;
  double Start, Finish, Duration;
 
  setvbuf(stdout, 0, _IONBF, 0);
 
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);
  MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank);
 
  GridSize = sqrt((double)ProcNum);
  if (ProcNum != GridSize*GridSize) {
    if (ProcRank == 0) {
      printf ("Number of processes must be a perfect square \n");
    }
  }
  else {
    if (ProcRank == 0)
      printf("Parallel matrix multiplication program\n");
 
    // Создание виртуальной решетки процессов и коммуникаторов 
    // строк и столбцов
    CreateGridCommunicators();
  
    // Выделение памяти и инициализация элементов матриц
    ProcessInitialization ( pAMatrix, pBMatrix, pCMatrix, pAblock,
      pBblock, pCblock, pMatrixAblock, Size, BlockSize );
    // Блочное распределение матриц между процессами
    DataDistribution(pAMatrix, pBMatrix, pMatrixAblock, pBblock, Size, 
      BlockSize);
 
    // Выполнение параллельного метода Фокса
    ParallelResultCalculation(pAblock, pMatrixAblock, pBblock, 
      pCblock, BlockSize);
 
    // Сбор результирующей матрицы на ведущем процессе
    ResultCollection(pCMatrix, pCblock, Size, BlockSize);
 
    // Завершение процесса вычислений
    ProcessTermination (pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock, 
      pCblock, pMatrixAblock);
  }
 
  MPI_Finalize();
  return 0;
}