Код:
// Iterations of the back substitution stage
for (int i=Size-1; i>=0; i--) {
// Calculating the rank of the process, which holds the pivot row
FindBackPivotRow(pParallelPivotPos[i], Size, IterProcRank, IterPivotPos);
// Calculating the unknown
if (ProcRank == IterProcRank) {
IterResult = pProcVector[IterPivotPos]/pProcRows[IterPivotPos*Size+i];
pProcResult[IterPivotPos] = IterResult;
}
// Broadcasting the value of the current unknown
MPI_Bcast(&IterResult, 1, MPI_DOUBLE, IterProcRank, MPI_COMM_WORLD);
// Updating the values of the vector b
for (int j=0; j<RowNum; j++)
if ( pProcPivotIter[j] < i ) {
val = pProcRows[j*Size + i] * IterResult;
pProcVector[j]=pProcVector[j] - val;
}
}
}
void TestDistribution(double* pMatrix, double* pVector, double* pProcRows,
double* pProcVector, int Size, int RowNum) {
if (ProcRank == 0) {
printf("Ishodnaya matrica: \n");
PrintMatrix(pMatrix, Size, Size);
printf("Ishodnoy Vector: \n");
PrintVector(pVector, Size);
}
MPI_Barrier(MPI_COMM_WORLD);
for (int i=0; i<ProcNum; i++) {
if (ProcRank == i) {
printf("\nProcRank = %d \n", ProcRank);
printf(" Matrix Stripe:\n");
PrintMatrix(pProcRows, RowNum, Size);
printf(" Vector: \n");
PrintVector(pProcVector, RowNum);
}
MPI_Barrier(MPI_COMM_WORLD);
}
}
// Function for the execution of the parallel Gauss algorithm
void ParallelResultCalculation(double* pProcRows, double* pProcVector,
double* pProcResult, int Size, int RowNum) {
ParallelGaussianElimination (pProcRows, pProcVector, Size, RowNum);
ParallelBackSubstitution (pProcRows, pProcVector, pProcResult, Size, RowNum);
}
}
// Function for computational process termination
void ProcessTermination (double* pMatrix, double* pVector, double* pResult,
double* pProcRows, double* pProcVector, double* pProcResult) {
if (ProcRank == 0) {
delete [] pMatrix;
delete [] pVector;
delete [] pResult;
}
delete [] pProcRows;
delete [] pProcVector;
delete [] pProcResult;
delete [] pParallelPivotPos;
delete [] pProcPivotIter;
delete [] pProcInd;
delete [] pProcNum;
}
// Function for testing the result
void TestResult(double* pMatrix, double* pVector, double* pResult, int Size) {
/* Buffer for storing the vector, that is a result of multiplication
of the linear system matrix by the vector of unknowns */
double* pRightPartVector;
// Flag, that shows wheather the right parts vectors are identical or not
int equal = 0;
double Accuracy = 1.e-6; // Comparison accuracy
if (ProcRank == 0) {
pRightPartVector = new double [Size];
for (int i=0; i<Size; i++) {
pRightPartVector[i] = 0;
for (int j=0; j<Size; j++) {
pRightPartVector[i] += pMatrix[i*Size+j]*pResult[pParallelPivotPos[j]];
}
}
for (int i=0; i<Size; i++) {
if (fabs(pRightPartVector[i]-pVector[i]) > Accuracy)
equal = 1;
}
if (equal == 1)
printf("The result of the parallel Gauss algorithm is NOT correct. Check your code.");
else
printf("The result of the parallel Gauss algorithm is correct.");
delete [] pRightPartVector;
}
}
int main(int argc, char* argv[])
{
double* pMatrix; // Matrix of the linear system
double* pVector; // Right parts of the linear system
double* pResult; // Result vector
double *pProcRows; // Rows of the matrix A
double *pProcVector; // Block of the vector b
double *pProcResult; // Block of the vector x
int Size; // Size of the matrix and vectors
int RowNum; // Number of the matrix rows
double start, finish, duration;
setvbuf(stdout, 0, _IONBF, 0);
MPI_Init ( &argc, &argv );
MPI_Comm_rank ( MPI_COMM_WORLD, &ProcRank );
MPI_Comm_size ( MPI_COMM_WORLD, &ProcNum );
if (ProcRank == 0)
printf("Paralelniy algoritm Gausa\n");
// Memory allocation and data initialization
ProcessInitialization(pMatrix, pVector, pResult,
pProcRows, pProcVector, pProcResult, Size, RowNum);
// The execution of the parallel Gauss algorithm
DataDistribution(pMatrix, pProcRows, pVector, pProcVector, Size, RowNum);
ParallelResultCalculation(pProcRows, pProcVector, pProcResult, Size, RowNum);
ResultCollection(pProcResult, pResult);
TestResult(pMatrix, pVector, pResult, Size);
// Computational process termination
ProcessTermination(pMatrix, pVector, pResult, pProcRows, pProcVector, pProcResult);
MPI_Finalize();