В данном разделе приведена лишь часть реализации, которая является наиболее значимой с точки зрения разработки.
Основная функция
int ProcNum = 0; // Количество используемых процессов
int ProcRank = 0; // Ранг текущего процесса
int GridSize; // Размер виртуальной решетки
int GridCoords[2]; // координаты текущего процесса в решетке
MPI_Comm GridComm; // коммуникатор-решетка
MPI_Comm ColComm; // Коммуникатор столбца
MPI_Comm RowComm;
void main(int argc, char* argv[]) {
double* pAMatrix; // Матрица А
double* pBMatrix; // Матрица В
double* pCMatrix; // Результат
int Size; // Размер матриц
int BlockSize; // Размер блоков
double* pAblock; // блок матрицы А
double* pBblock; // блок матрицы В
double* pCblock; // блок результирующей матрицы
double* TmpBlock;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);
MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank);
GridSize = (int)sqrt((double)ProcNum);
// Создаем коммуникаторы
Create_Communicators();
// выделение памяти и инициализация элементов массивов
ProcessInitialization( pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock, pCblock, TmpBlock, Size, BlockSize );
//разделение данных
DataDistribution(pAMatrix, pBMatrix, TmpBlock, pBblock, Size, BlockSize);
// Выполнение параллельного метода Фокса
ParallelResultCalculation(pAblock, TmpBlock, pBblock, pCblock, BlockSize);
//сбор результата
ResultCollection(pCMatrix, pCblock, Size, BlockSize);
// Завершение процесса вычислений
ProcessTermination(pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock, pCblock, TmpBlock);
}
MPI_Finalize();
}


CreateGridCommunicators. Данная функция создает коммуникатор в виде двумерной квадратной решетки, определяет координаты каждого процесса в этой решетке, а также создает коммуникаторы отдельно для каждой строки и каждого столбца.

Функция ProcessInitialization. Данная функция определяет параметры решаемой задачи (размеры матриц и их блоков), выделяет память для хранения данных и осуществляет ввод исходных матриц (или формирует их при помощи какого-либо датчика случайных чисел).

Функции DataDistribution и ResultCollection. После задания исходных матриц на нулевом процессе необходимо осуществить распределение исходных данных. Для этого предназначена функция DataDistribution.

void Distribute(double* Matrix, double* Block, int Size, int BlockSize)
{
double* Buff = new double [BlockSize*Size];// буфер для рассылки блоков
if (GridCoords[1] == 0)
MPI_Scatter(Matrix, BlockSize*Size, MPI_DOUBLE, Buff, BlockSize*Size, MPI_DOUBLE, 0, ColComm);//распределение строк и столбцов между первыми проц-ми
for (int i=0; i {
MPI_Scatter(&Buff[i*Size], BlockSize, MPI_DOUBLE, &Block[i*BlockSize], BlockSize, MPI_DOUBLE, 0, RowComm);//дальнейшее распределение строк между процессами
}
delete [] Buff;
}

void DataDistribution(double* pAMatrix, double* pBMatrix, double* TmpBlock, double* pBblock, int Size, int BlockSize)
{
Distribute(pAMatrix, TmpBlock, Size, BlockSize);
Distribute(pBMatrix, pBblock, Size, BlockSize);
}

Для выполнения сбора результирующей матрицы из блоков предназначена функция ResultCollection.

void ResultCollection(double* pCMatrix, double* pCblock, int Size, int BlockSize)
{
double * Buff = new double [Size*BlockSize];
for (int i=0; i { // сборка матрицы С из блоков
MPI_Gather( &pCblock[i*BlockSize], BlockSize, MPI_DOUBLE, &Buff[i*Size], BlockSize, MPI_DOUBLE, 0, RowComm);
}
if (GridCoords[1] == 0)
{ // сборка из строк
MPI_Gather(Buff, BlockSize*Size, MPI_DOUBLE, pCMatrix, BlockSize*Size, MPI_DOUBLE, 0, ColComm);
} delete [] Buff;
}

Для непосредственного выполнения параллельного алгоритма Фокса умножения матриц предназначена функция ParallelResultCalculation, которая реализует логику работы алгоритма.

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);
}
}

Более подробное описание и реализацию некоторых методов можно найти на сайте intuin.ru