Ниже приводится исходный код программы, реализующей параллельный алгоритм перемножения матриц на основе алгоритма Фокса:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <mpi.h>
int ProcNum = 0; // Number of available processes
int ProcRank = 0; // Rank of current process
int GridSize; // Size of virtual processor grid
int GridCoords[2]; // Coordinates of current processor in grid
MPI_Comm GridComm; // Grid communicator
MPI_Comm ColComm; // Column communicator
MPI_Comm RowComm; // Row communicator
// Function for simple initialization of matrix elements
void DummyDataInitialization (double* pAMatrix, double* pBMatrix,int Size){
int i, j; // Loop variables
for (i=0; i<Size; i++)
for (j=0; j<Size; j++) {
pAMatrix[i*Size+j] = 1;
pBMatrix[i*Size+j] = 1;
}
}
// Function for random initialization of matrix elements
void RandomDataInitialization (double* pAMatrix, double* pBMatrix,
int Size) {
int i, j; // Loop variables
srand(unsigned(clock()));
for (i=0; i<Size; i++)
for (j=0; j<Size; j++) {
pAMatrix[i*Size+j] = rand()/double(1000);
pBMatrix[i*Size+j] = rand()/double(1000);
}
}
// Function for formatted matrix output
void PrintMatrix (double* pMatrix, int RowCount, int ColCount) {
int i, j; // Loop variables
for (i=0; i<RowCount; i++) {
for (j=0; j<ColCount; j++)
printf("%7.4f ", pMatrix[i*ColCount+j]);
printf("\n");
}
}
// Function for matrix multiplication
void SerialResultCalculation (double* pAMatrix, double* pBMatrix,
double* pCMatrix, int Size) {
int i, j, k; // Loop variables
for (i=0; i<Size; i++) {
for (j=0; j<Size; j++)
for (k=0; k<Size; k++)
pCMatrix[i*Size+j] += pAMatrix[i*Size+k]*pBMatrix[k*Size+j];
}
}
// Function for block multiplication
void BlockMultiplication (double* pAblock, double* pBblock,
double* pCblock, int Size) {
SerialResultCalculation(pAblock, pBblock, pCblock, Size);
}
// Creation of two-dimensional grid communicator
// and communicators for each row and each column of the grid
void CreateGridCommunicators () {
int DimSize[2]; // Number of processes in each dimension of the grid
int Periodic[2]; // =1, if the grid dimension should be periodic
int Subdims[2]; // =1, if the grid dimension should be fixed
DimSize[0] = GridSize;
DimSize[1] = GridSize;
Periodic[0] = 0;
Periodic[1] = 0;
// Creation of the Cartesian communicator
MPI_Cart_create(MPI_COMM_WORLD, 2, DimSize, Periodic, 1, &GridComm);
// Determination of the cartesian coordinates for every process
MPI_Cart_coords(GridComm, ProcRank, 2, GridCoords);
// Creating communicators for rows
Subdims[0] = 0; // Dimensionality fixing
Subdims[1] = 1; // The presence of the given dimension in the subgrid
MPI_Cart_sub(GridComm, Subdims, &RowComm);
// Creating communicators for columns
Subdims[0] = 1;
Subdims[1] = 0;
MPI_Cart_sub(GridComm, Subdims, &ColComm);
}
// Function for memory allocation and data initialization
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("\nEnter size of the initial objects: ");
scanf("%d", &Size);
if (Size%GridSize!= 0) {
printf ("Size of matricies must be divisible by the grid size!\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];
DummyDataInitialization(pAMatrix, pBMatrix, Size);
//RandomDataInitialization(pAMatrix, pBMatrix, Size);
}
}
// Function for checkerboard matrix decomposition
void CheckerboardMatrixScatter (double* pMatrix, double* pMatrixBlock,
int Size, int BlockSize) {
double * MatrixRow = new double [BlockSize*Size];
if (GridCoords[1] == 0) {
MPI_Scatter(pMatrix, BlockSize*Size, MPI_DOUBLE, MatrixRow,
BlockSize*Size, MPI_DOUBLE, 0, ColComm);
}
for (int i=0; i<BlockSize; i++) {
MPI_Scatter(&MatrixRow[i*Size], BlockSize, MPI_DOUBLE,
&(pMatrixBlock[i*BlockSize]), BlockSize, MPI_DOUBLE, 0, RowComm);
}
delete [] MatrixRow;
}
// Data distribution among the processes
void DataDistribution (double* pAMatrix, double* pBMatrix, double*
pMatrixAblock, double* pBblock, int Size, int BlockSize) {
// Scatter the matrix among the processes of the first grid column
CheckerboardMatrixScatter(pAMatrix, pMatrixAblock, Size, BlockSize);
CheckerboardMatrixScatter(pBMatrix, pBblock, Size, BlockSize);
}
// Function for gathering the result matrix
void ResultCollection (double* pCMatrix, double* pCblock, int Size,
int BlockSize) {
double * pResultRow = new double [Size*BlockSize];
for (int i=0; i<BlockSize; i++) {
MPI_Gather(&pCblock[i*BlockSize], BlockSize, MPI_DOUBLE,
&pResultRow[i*Size], BlockSize, MPI_DOUBLE, 0, RowComm);
}
if (GridCoords[1] == 0) {
MPI_Gather(pResultRow, BlockSize*Size, MPI_DOUBLE, pCMatrix,
BlockSize*Size, MPI_DOUBLE, 0, ColComm);
}
delete [] pResultRow;
}
// Broadcasting matrix A blocks to process grid rows
void ABlockCommunication (int iter, double *pAblock, double* pMatrixAblock,
int BlockSize) {
// Defining the leading process of the process grid row
int Pivot = (GridCoords[0] + iter) % GridSize;
// Copying the transmitted block in a separate memory buffer
if (GridCoords[1] == Pivot) {
for (int i=0; i<BlockSize*BlockSize; i++)
pAblock[i] = pMatrixAblock[i];
}
// Block broadcasting
MPI_Bcast(pAblock, BlockSize*BlockSize, MPI_DOUBLE, Pivot, RowComm);
}
// Cyclic shift of matrix B blocks in the process grid columns
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 ParallelResultCalculation (double* pAblock, double* pMatrixAblock,
double* pBblock, double* pCblock, int BlockSize) {
for (int iter = 0; iter < GridSize; iter ++) {
// Sending blocks of matrix A to the process grid rows
ABlockCommunication (iter, pAblock, pMatrixAblock, BlockSize);
// Block multiplication
BlockMultiplication(pAblock, pBblock, pCblock, BlockSize);
// Cyclic shift of blocks of matrix B in process grid columns
BblockCommunication(pBblock, BlockSize);
}
}
// Test printing of the matrix block
void TestBlocks (double* pBlock, int BlockSize, char str[]) {
MPI_Barrier(MPI_COMM_WORLD);
if (ProcRank == 0) {
printf("%s \n", str);
}
for (int i=0; i<ProcNum; i++) {
if (ProcRank == i) {
printf ("ProcRank = %d \n", ProcRank);
PrintMatrix(pBlock, BlockSize, BlockSize);
}
MPI_Barrier(MPI_COMM_WORLD);
}
}
void TestResult (double* pAMatrix, double* pBMatrix, double* pCMatrix,
int Size) {
double* pSerialResult; // Result matrix of serial multiplication
double Accuracy = 1.e-6; // Comparison accuracy
int equal = 0; // =1, if the matrices are not equal
int i; // Loop variable
if (ProcRank == 0) {
pSerialResult = new double [Size*Size];
for (i=0; i<Size*Size; i++) {
pSerialResult[i] = 0;
}
BlockMultiplication(pAMatrix, pBMatrix, pSerialResult, Size);
for (i=0; i<Size*Size; i++) {
if (fabs(pSerialResult[i]-pCMatrix[i]) >= Accuracy)
equal = 1;
}
if (equal == 1)
printf("The results of serial and parallel algorithms are NOT "
" identical. Check your code.");
else
printf("The results of serial and parallel algorithms are "
" identical.");
}
}
// Function for computational process termination
void ProcessTermination (double* pAMatrix, double* pBMatrix,
double* pCMatrix, double* pAblock, double* pBblock, double* pCblock,
double* pMatrixAblock) {
if (ProcRank == 0) {
delete [] pAMatrix;
delete [] pBMatrix;
delete [] pCMatrix;
}
delete [] pAblock;
delete [] pBblock;
delete [] pCblock;
delete [] pMatrixAblock;
}
void main (int argc, char* argv[]) {
double* pAMatrix; // The first argument of matrix multiplication
double* pBMatrix; // The second argument of matrix multiplication
double* pCMatrix; // The result matrix
int Size; // Size of matricies
int BlockSize; // Sizes of matrix blocks on current process
double *pAblock; // Initial block of matrix A on current process
double *pBblock; // Initial block of matrix B on current process
double *pCblock; // Block of result matrix C on current process
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");
// Creating the cartesian grid, row and column communcators
CreateGridCommunicators();
// Memory allocation and initialization of matrix elements
ProcessInitialization (pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock,
pCblock, pMatrixAblock, Size, BlockSize);
DataDistribution(pAMatrix, pBMatrix, pMatrixAblock, pBblock, Size,
BlockSize);
// Execution of Fox method
ParallelResultCalculation(pAblock, pMatrixAblock, pBblock,
pCblock, BlockSize);
ResultCollection(pCMatrix, pCblock, Size, BlockSize);
TestResult(pAMatrix, pBMatrix, pCMatrix, Size);
// Process Termination
ProcessTermination (pAMatrix, pBMatrix, pCMatrix, pAblock, pBblock,
pCblock, pMatrixAblock);
}
MPI_Finalize();
}
В программе можно выделить несколько функциональных модулей:
