Параллельная версия

// GaussParallel.cpp : Defines the entry point for the console application.

#include "stdafx.h"

//-----------------------Параллельная версия метода Гаусса----------------------

//------------------------------------------------------------------------------

#include "stdafx.h"

#include < stdio.h>

#include < stdlib.h>

#include < conio.h>

#include < time.h>

#include < math.h>

#include < mpi.h>

enum solution {NO_SOLUTION,ONE_SOLUTION,MANY_SOLUTIONS};

int g_size; //размерность исходной системы

double *g_matrix; //матрица коэффициентов системы в виде 1 вектора

double *g_vector; //вектор свободных членов

double *g_result; //вектор значений неизвестных

solution g_isConsistent; //признак совместности системы

double *g_processMatrix; //подматрица коэффициентов системы на конкретном процессе

double *g_processVector; //элементы вектора свободных членов на конкретном процессе

double *g_processResult; //компоненты вектора искомых значений неизвестных на конкретном процессе

int g_rowsCount; /* число строк исходной системы, хранящихся на конкретном процессе

(пересчитывается во время распределения исходных данных) */

double g_startTime, //время начала работы алгоритма

g_finishTime, //время окончания работы алгоритма

g_solvingTime; //время завершения работы алгоритма

int *g_pivotRowIndex; //индекс ведущей строки системы на каждой итерации

int *g_processIterationNumber; /* номер итерации, на которой

соответствующая строка системы, расположенная на процессе, выбрана в качестве ведущей */

int g_processCount; //число используемых процессов

int g_processRank; //ранг текущего процесса

int *g_processFirstRowGlobalIndex; //глобальный индекс первой строки системы на конкретном процессе

int *g_processRowsCount; //число строк системы на конкретном процессе

double default_value = 0.0;

const int BUFFER_SIZE = 65536;

/// Загрузка данных из файла

void LoadTask()

{

FILE* file;

fopen_s(&file, "..\\..\\input.txt", "r");

int count = 0;

/// Считаем количество строк

char *s = new char [BUFFER_SIZE];

while (!feof(file))

{

fgets(s, BUFFER_SIZE, file);

count++;

}

delete[] s;

g_size = count-1;

/// Выделяем память под исходные данные

g_matrix = new double [g_size*g_size];

g_vector = new double [g_size];

g_result = new double [g_size];

fseek(file, 0, SEEK_SET);

char *buf = new char [BUFFER_SIZE]; ///< Буфер строки

/// Формируем сисему линейных уравнений

for (int i=0; i< g_size; i++)

{

int k = 1; ///< Индекс пробела после числа в строке

int pos = 1; ///< Позиция пробела перед числом в строке

/// Читаем строку из файла в буфер

fgets(buf, BUFFER_SIZE, file);

/// Запоминаем строку

for (int j = 0; j< =g_size; j++)

{

/// Читаем до пробела

while (buf[k]&&(buf[k]!=' ')) k++;

char *number = new char [k-pos+1]; ///< Буфер для хранения числа

for (int m = pos; m< k; m++)

number[m-pos] = buf[m]; ///< Записываем число из общей строки

number[k-pos] = '\0'; ///< Последний символ - конец строки

char *p = NULL;

if (j!=g_size) /// Если число - элемент метрицы

g_matrix[g_size*i+j] = strtod(number, &p);

else /// Если нет - записываем в вектор свободных коэффициентов

g_vector[i] = strtod(number, &p);

k++;

pos = k;

delete[] number;

}

delete[] buf;

//Закрываем файл исходных данных

fclose(file);

}

void ProcessInitialization()

{

int RestRows; ///< Число строк, ещё не распределённых по процессам

if (g_processRank == 0)

{

printf("Loading task...");

LoadTask();

}

/// Рассылка размера системы всем процессам

MPI_Bcast(&g_size, 1, MPI_INT, 0, MPI_COMM_WORLD);

/// Определение размера части данных, расположенных на конкретном процессе

RestRows = g_size;

g_rowsCount = floor((double)RestRows/(double)(g_processCount)); ///< Каждый процесс должен посчитать минимум такое количество срок.

if((RestRows - g_rowsCount*g_processCount)< =g_processRank)

g_rowsCount++;

/// Выделение памяти под часть данных, расположенных на конкретном процессе

g_processMatrix = new double [g_rowsCount*g_size];

g_processVector = new double [g_rowsCount];

g_processResult = new double [g_rowsCount];

g_pivotRowIndex = new int [g_size];

g_processIterationNumber = new int [g_rowsCount];

g_processFirstRowGlobalIndex = new int [g_processCount];

g_processRowsCount = new int [g_processCount];

for (int i=0; i< g_rowsCount; i++)

g_processIterationNumber[i] = -1;

}

void InitialSystemDistribution()

{

int *FirstIndexSending; ///< Индекс первого элемента матрицы, передаваемого процессу

int *SendingSize; ///< Число элементов матрицы, передаваемых процессу

int CountProcessRows; ///< Число строк, отдаваемых текущему процессу

int RestRows; ///< Количество ещё не разосланных строк

FirstIndexSending = new int [g_processCount];

SendingSize = new int [g_processCount];

RestRows = g_size;

/// Устанавливаем параметры для передачи строк

for (int i=0; i< g_processCount; i++)

{

int PreviousSize = (i==0) ? 0 : SendingSize[i-1];

int PreviousIndex = (i==0) ? 0 : FirstIndexSending[i-1];

int Portion = (i==0) ? 0 : CountProcessRows; /// Количество строк для процесса

RestRows -= Portion;

CountProcessRows = RestRows/(g_processCount-i);

SendingSize[i] = CountProcessRows*g_size;

FirstIndexSending[i] = PreviousIndex+PreviousSize;

}

/// Рассылаем матрицу процессам

MPI_Scatterv(g_matrix, SendingSize, FirstIndexSending, MPI_DOUBLE, g_processMatrix,

SendingSize[g_processRank], MPI_DOUBLE, 0, MPI_COMM_WORLD);

RestRows = g_size;

/// Считаем параметры для пересылки вектора свободных коэффициентов

for (int i=0; i< g_processCount; i++)

{

int PreviousSize = i==0 ? 0 : g_processRowsCount[i-1];

int PreviousIndex = i==0 ? 0 : g_processFirstRowGlobalIndex[i-1];

int Portion = i==0 ? 0 : g_processRowsCount[i-1];

RestRows -= Portion;

g_processRowsCount[i] = RestRows/(g_processCount-i);

g_processFirstRowGlobalIndex[i] = PreviousIndex+PreviousSize;

}

/// Рассылаем вектор cвободных коэффициентов

MPI_Scatterv(g_vector, g_processRowsCount, g_processFirstRowGlobalIndex, MPI_DOUBLE, g_processVector,

g_processRowsCount[g_processRank], MPI_DOUBLE, 0, MPI_COMM_WORLD);

g_isConsistent = ONE_SOLUTION;

/// Рассылаем признак наличия решения

MPI_Bcast(&g_isConsistent, 1, MPI_INT, 0, MPI_COMM_WORLD);

delete [] SendingSize;

delete [] FirstIndexSending;

}

//------------------------------------------------------------------------------

/// Исключение неизвестных в столбце матрицы с заданным номером

void ColumnElimination(int IterationNumber, double *PivotSystemRow)

{

double multiplier;

for (int i=0; i< g_rowsCount; i++)

{

if (g_processIterationNumber[i] == -1)

{

multiplier = g_processMatrix[i*g_size+IterationNumber] / PivotSystemRow[IterationNumber];

for (int j=IterationNumber; j< g_size; j++)

{

g_processMatrix[i*g_size + j] -= PivotSystemRow[j]*multiplier;

}

g_processVector[i] -= PivotSystemRow[g_size]*multiplier;

}

//------------------------------------------------------------------------------

//Прямой ход метода Гаусса

void ForwardGauss()

{

int PivotRowPosition; // индекс ведущей строки на конкретном процессе

/* Структура, хранящая

локально максимальный по модулю элемент и

ранг соответствующего процесса */

struct { double MaxValue; int ProcRank; } ProcessPivot, //локальные значения

Pivot; //глобальные значения

/* pPivotRow хранит строку матрицы коэффициентов,

и соответствующий элемент вектора свободных членов */

double* pPivotRow = new double [g_size];

// Цикл гауссовых преобразований

for (int i=0; i< g_size; i++)

{

// Определение локальной ведущей строки

double MaxValue = 0;

int tmp_index = -1;

for (int j=0; j< g_rowsCount; j++)

{

tmp_index = j;

if ((g_processIterationNumber[j] == -1)&&

(MaxValue < fabs(g_processMatrix[i+g_size*j])))

{

MaxValue = fabs(g_processMatrix[i+g_size*j]);

PivotRowPosition = j;

}

//Запоминание локальных максимальных значений и рангов соответствующих процессов в единую структуру

ProcessPivot.MaxValue = MaxValue;

ProcessPivot.ProcRank = g_processRank;

/// Определение максимального по модулю элемента

MPI_Allreduce(&ProcessPivot, &Pivot, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD);

///Вычисление номера ведущей строки всей системы

if ( g_processRank == Pivot.ProcRank )

{

/// Если самое большое значение по модулю == 0 -> нулевой столбец -> много решений

if (Pivot.MaxValue == 0)

{

g_isConsistent = MANY_SOLUTIONS;

MPI_Barrier(MPI_COMM_WORLD);

MPI_Send(&g_isConsistent, MANY_SOLUTIONS, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD);

g_processIterationNumber[tmp_index] = i;

g_pivotRowIndex[i] = g_processFirstRowGlobalIndex[g_processRank] + PivotRowPosition;

continue;

}

else

{

g_processIterationNumber[PivotRowPosition] = i; /* Запоминаем номер итерации, на которой

строка с локальным номером

является ведущей для всей системы */

g_pivotRowIndex[i] = g_processFirstRowGlobalIndex[g_processRank] + PivotRowPosition;

}

//Рассылка всем процессам вычисленного номера ведущей строки системы

MPI_Bcast(&g_pivotRowIndex[i], 1, MPI_INT, Pivot.ProcRank, MPI_COMM_WORLD);

if (g_processRank == Pivot.ProcRank)

{

//Заполнение ведущей строки системы

for (int j=0; j< g_size; j++)

pPivotRow[j] = g_processMatrix[PivotRowPosition*g_size + j];

pPivotRow[g_size] = g_processVector[PivotRowPosition];

}

//Рассылка ведущей строки всем процессам

MPI_Bcast(pPivotRow, g_size+1, MPI_DOUBLE, Pivot.ProcRank, MPI_COMM_WORLD);

//Исключение неизвестных в столбце с номером i

ColumnElimination(i, pPivotRow);

}

void Find_Rank_and_Position(int RowIndex, int &IterationProcessRank, int &IterationPivotRowIndex)

{

/// Определяем ранг процесса, содержащего заданную строку

for (int i=0; i< g_processCount-1; i++)

{

if ((g_processFirstRowGlobalIndex[i]< =RowIndex) && (RowIndex< g_processFirstRowGlobalIndex[i+1]))

IterationProcessRank = i;

}

if (RowIndex >= g_processFirstRowGlobalIndex[g_processCount-1])

IterationProcessRank = g_processCount-1;

/// Определяем номер строки на процессе

IterationPivotRowIndex = RowIndex - g_processFirstRowGlobalIndex[IterationProcessRank];

}

void BackGauss()

{

int IterationProcessRank; // Ранг процесса, хранящего текущую ведущую строку

int IterationPivotRowIndex; // Номер текущей ведущей строки на содержащем её процессе

double IterResult; // Вычисленное значение текущей неизвестной

double val;

// Основной цикл обратного вычислительного процесса

for (int i=g_size-1; i>=0; i--)

{

// Определяем ранг процесса, содержащего текущую ведущую строку, и номер этой строки на процессе

Find_Rank_and_Position(g_pivotRowIndex[i], IterationProcessRank, IterationPivotRowIndex);

// Вычисляем значение неизвестной

if (g_processRank == IterationProcessRank)

{

if (g_processMatrix[IterationPivotRowIndex*g_size+i] == 0)

{

if (g_processVector[IterationPivotRowIndex] == 0)

IterResult = 0.0;

else

{

g_isConsistent = NO_SOLUTION;

MPI_Barrier(MPI_COMM_WORLD);

MPI_Send(&g_isConsistent, 1, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD);

break;

}

else

IterResult = g_processVector[IterationPivotRowIndex]/g_processMatrix[IterationPivotRowIndex*g_size+i];

g_processResult[IterationPivotRowIndex] = IterResult;

}

//Рассылка всем остальным процессам найденного значения переменной

MPI_Bcast(&IterResult, 1, MPI_DOUBLE, IterationProcessRank, MPI_COMM_WORLD);

//Корректировка элементов вектора свободных членов

for (int j=0; j< g_rowsCount; j++)

if (g_processIterationNumber[j] < i)

{

val = g_processMatrix[g_size*j + i] * IterResult;

g_processVector[j] -= val;

}

/// Решение исходной СЛУ

void SolvingLinearSystem()

{

ForwardGauss();

BackGauss();

MPI_Gatherv(g_processResult, g_processRowsCount[g_processRank], MPI_DOUBLE, g_result,

g_processRowsCount, g_processFirstRowGlobalIndex, MPI_DOUBLE, 0, MPI_COMM_WORLD);

}

/// Удаление всех массивов, под которые была выделена память в процессе работы программы

void DeletingProcessesData()

{

if (g_processRank == 0)

{

delete [] g_matrix;

delete [] g_vector;

delete [] g_result;

}

delete [] g_processMatrix;

delete [] g_processVector;

delete [] g_processResult;

delete [] g_pivotRowIndex;

delete [] g_processIterationNumber;

delete [] g_processFirstRowGlobalIndex;

delete [] g_processRowsCount;

}

void _tmain(int argc, char* argv[])

{

MPI_Init(&argc, &argv);

MPI_Comm_rank (MPI_COMM_WORLD, &g_processRank);

MPI_Comm_size (MPI_COMM_WORLD, &g_processCount);

printf("process #%d started \n", g_processRank);

MPI_Barrier(MPI_COMM_WORLD);

if (g_processRank == 0)

printf("Gauss method of solving linear systems. Parallel version...\n\n");

MPI_Barrier(MPI_COMM_WORLD);

ProcessInitialization();

InitialSystemDistribution();

MPI_Barrier(MPI_COMM_WORLD);

g_startTime = clock();

SolvingLinearSystem();

MPI_Barrier(MPI_COMM_WORLD);

g_finishTime = clock();

g_solvingTime = (g_finishTime-g_startTime)/CLOCKS_PER_SEC;

if (g_processRank == 0)

{

if (g_isConsistent == ONE_SOLUTION)

{

printf("\nSystem has one solution. \n");

}

if (g_isConsistent == NO_SOLUTION)

printf("\nSystem has no solutions.\n");

else if (g_isConsistent == MANY_SOLUTIONS)

{

printf("\nSystem has infinitely many solutions. \n");

}

printf("\nTime of solving the system: %f\n", g_solvingTime);

}

//Удаление всех массивов

DeletingProcessesData();

if (g_processRank == 0)

{

printf("Press any key to exit...");

_getch();

}

MPI_Finalize();

}

Параллельная версия

Новости