Новости
О Центре
Кластер
Обучение
Основной курс по параллельному программированию
Учебные курсы
Магистратура
Дополнительное образование
Работы студентов
Библиотека
Исследования
Конференции
Полезные ссылки
NVIDIA
Контакты
О сайте
Имя:
Пароль:
запомнить:
Забыли пароль? Регистрация

Введение

Решение дифференциальных уравнений, является неотъемлемой частью анализа многих математических моделей.
Не всегда есть возможность решить их аналитически, а если есть, то это может быть весьма трудоёмкой задачей.
Использование численных методов помогает избавиться от лишней сложности и получить решение для целого класса задач, не решаемых аналитически.
За счёт использования высокопроизводительных вычислений и большего числа процессоров можно достигнуть большей точности/скорости вычислений за одно и то же время, что позволяет наращивать производительность в зависимости от приоритета задачи.

Постановка задачи

Рассмотрим в качестве примера проблему численного решения задачи Дирихле для уравнения Пуассона, определяемую как задачу нахождения функции u=u(x,y), удовлетворяющей в области определения уравнению: 
 

На границе области функция u принимает значения:



функции f(x,y) и g(x,y) считаются данными. 

Для решения уравнений на сетке используют метод конечных разностей.
Явная разностная схема решения этого уравнения имеет вид:



Для реализации её могут быть использованы разные вычислительные методы. 
Предлагается использовать алгоритм Гаусса-Зейделя.

Последовательная версия алгоритма

Алгоритм Гаусса-Зейделя подразумевает построчный (либо постолбцовый) обход сетки и вычисление следующего значения функции в точке по схеме:
double A[n][m]; // сетка текущих значений функции, она же результирующая
double F[n][m]; // сетка из производных

...

for(j=1;j<=m-1;j++)
	for(i=1;i<=n-1;i++)
		A[i][j] = 0.25*(A[i-1][j] + A[i+1][j] + A[i][j-1]+ A[i][j+1] - h*h*F[i][j]) 

Параллельная версия алгоритма

Исходная сетка распределяется между процессами по строкам или столбцам полосами, в зависимости от способа хранения и обхода. 
Для определённости разделим по строкам. 


По итогам каждой из итераций процессы обмениваются граничными строками между собой по следующей схеме:

Реализация

// Решать мы можем с двойной или одинарной точностью 
// тип CellPoint хранит значение в точке разностной сетки
// За проверку соответствия типа CellPoint и MPI_CELLPOINT отвечают assert'ы в полном коде программы

typedef float CellPoint;
#define MPI_CELLPOINT MPI_FLOAT

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

typedef CellPoint (*TwoAgrumentsFunction) (CellPoint x, CellPoint y);

// координаты квадратной области, на которой будет решаться диффур
struct Field
{
	CellPoint x;
	CellPoint y;
	CellPoint length;
};

// в качестве тестовой функции взят эллиптический параболоид с максимумом в точке (0,0)

CellPoint g(CellPoint x, CellPoint y)
{
	return 5.0 + (-x*x - y*y)/5.0;
}

CellPoint f(CellPoint x, CellPoint y)
{
	return (CellPoint)-4.0/5.0;}

//Параметры алгоритма задаются извне
struct AlgorithmParameters
{
	UINT gridSize;       // размер сетки
	Field field;         // область на которой решаем
	CellPoint precision; // требуемая точность
	CellPoint stepSize;  // размер шага сетки
};

// для параллельной части программы важно знать ранк текущего процесса и общее число процессов
struct ProcessParameters
{
	int num;
	int rank;
};

// У каждого процесса своя часть данных
struct OwnProcessData
{
	CellPoint * matrixPiece; // рабочая матрица
	CellPoint * derivative;  // матрица производных
	UINT rowNum;             // число строк в полосе
};

//Все оперции синхронизации и обмена данными вынесены в отдельные процедуры

int main(int argc, char ** argv )
{
	MPI_Init(&argc, &argv);

	AlgorithmParameters params;
	ProcessParameters process = GetProcessParameters();

	if(process.rank == MAIN_THREAD_RANK)
	{
		params = GetAlgorithmParametersFromConsole(process);
	}
	BroadcastAlgorithmParameters(¶ms);

	OwnProcessData ownProcData = InitializateWithPredefinedData(params, process);

	ParallelTimer timer;
	timer.Start();
	{
		SolvePDEProblem(params, process, ownProcData );
	}
	timer.Stop();

	CollectResult(  params, process, ownProcData );

	DWORD workTime = timer.GetResult();
	if (process.rank == MAIN_THREAD_RANK)
	{
		std::cout<<"work time: "<a;
		MPI_Barrier(MPI_COMM_WORLD);
#endif
		MPI_Allreduce(&procDelta, &maxDelta, 1, MPI_CELLPOINT, MPI_MAX, MPI_COMM_WORLD);
	} while (maxDelta > params.precision);
}

void CollectResult(AlgorithmParameters params, ProcessParameters process, OwnProcessData ownData)
{
	MPI_Gather ( /*sendbuf*/ ownData.matrixPiece+params.gridSize,
				 /*sendcount*/ params.gridSize * (ownData.rowNum - 2),
				 /*sendtype*/ MPI_CELLPOINT,
				 /*recvbuf*/ matrix+params.gridSize,
				 /*recvcount*/  params.gridSize * (ownData.rowNum - 2),
				 /*recvtype*/ MPI_CELLPOINT,
				 /*root*/     MAIN_THREAD_RANK,
				 /*comm*/     MPI_COMM_WORLD);
}

void BroadcastAlgorithmParameters( AlgorithmParameters * input)
{
	MPI_Bcast( &input->gridSize, 1, MPI_UNSIGNED, MAIN_THREAD_RANK, MPI_COMM_WORLD );
	MPI_Bcast( &input->precision, 1, MPI_CELLPOINT, MAIN_THREAD_RANK, MPI_COMM_WORLD );
	MPI_Bcast( &input->field.length,1, MPI_UNSIGNED, MAIN_THREAD_RANK, MPI_COMM_WORLD );
	MPI_Bcast( &input->field.x, 1, MPI_CELLPOINT, MAIN_THREAD_RANK, MPI_COMM_WORLD );
	MPI_Bcast( &input->field.y, 1, MPI_CELLPOINT, MAIN_THREAD_RANK, MPI_COMM_WORLD );
	MPI_Bcast( &input->stepSize, 1, MPI_CELLPOINT, MAIN_THREAD_RANK, MPI_COMM_WORLD );
}

void SetupMatrix(AlgorithmParameters params, TwoAgrumentsFunction border, TwoAgrumentsFunction derivative )
{
	CellPoint h = params.field.length/(params.gridSize-1);

	matrix[0] = border(params.field.x, params.field.y);
	matrix[params.gridSize-1] = border(params.field.x + params.field.length, params.field.y);
	matrix[params.gridSize*(params.gridSize-1)] =  border(params.field.x, params.field.y + params.field.length);
	matrix[params.gridSize*params.gridSize - 1] = border(params.field.x + params.field.length, params.field.y + params.field.length);

	for ( UINT i = 1; i < params.gridSize -  1; i++ )
	{
		matrix [0 * params.gridSize + i] = border( params.field.x + i * h, params.field.y );
		matrix [i * params.gridSize + 0] = border( params.field.x, params.field.y + i*h);
		matrix [ (params.gridSize - 1 ) * params.gridSize + i ] = border ( params.field.x + i*h, params.field.y + params.field.length ); 
		matrix [i*params.gridSize + (params.gridSize - 1) ] = border ( params.field.x + params.field.length, params.field.y + i*h);
	}

	for ( UINT i = 1; i < params.gridSize - 1; i++)
	{
		for (UINT j = 1; j < params.gridSize - 1; j++ )
		{
			matrix[i*params.gridSize + j] = derivative( params.field.x + j*h, params.field.y + i*h );
		}
	}
}

OwnProcessData InitializateWithPredefinedData(AlgorithmParameters params, 
												 ProcessParameters process)
{
	//  init proc memory
	OwnProcessData ownData = {0};
	ownData.rowNum = ( params.gridSize - 2 )/process.num + 2;
	ownData.matrixPiece = new CellPoint[params.gridSize * ownData.rowNum];
	ownData.derivative =  new CellPoint[params.gridSize * ownData.rowNum];

	if ( process.rank == MAIN_THREAD_RANK )
	{
		//    InitializeMatrix
		matrix = new CellPoint[params.gridSize * params.gridSize];
		SetupMatrix(params, g, f);
		//    get first row
		for ( UINT i = 0; i < params.gridSize; i++ )
		{
			ownData.matrixPiece[i] = matrix [i];
		}
	}
	//Scatter matrix
	MPI_Scatter(/*sendbuf*/  matrix+params.gridSize, 
				/*sendcnt*/ (ownData.rowNum-2)*params.gridSize,
				/*sendtype*/ MPI_CELLPOINT,
				/*rectbuf*/  ownData.matrixPiece+params.gridSize,
				/*recvcnt*/  (ownData.rowNum-2)*params.gridSize,
				/*recvtype*/ MPI_CELLPOINT,
				/*rootrank*/ MAIN_THREAD_RANK,
				/*comm*/     MPI_COMM_WORLD );

	//    send last row
	if ( process.num > 1)
	{
		if (process.rank == MAIN_THREAD_RANK)
		{
			MPI_Send( /*buf*/     matrix + (params.gridSize-1)*params.gridSize,
					 /*count*/    params.gridSize,
					 /*datatype*/ MPI_CELLPOINT,
					 /*dst rank*/ process.num - 1,
					 /*tag*/      5,
					 /*comm*/     MPI_COMM_WORLD);
		}

		if (process.rank == process.num - 1)
		{
			MPI_Status status;
			MPI_Recv( /*buf*/ ownData.matrixPiece + (ownData.rowNum-1)*params.gridSize,
					  /*count*/ params.gridSize,
					  /*datatype*/ MPI_CELLPOINT,
					  /*src rank*/ MAIN_THREAD_RANK,
					  /*tag*/    5,
					  /*comm*/   MPI_COMM_WORLD,
					  /*status*/ &status );

		}
	}
	else
	{
		for (int i = 0; i < params.gridSize; i++)
		{
			ownData.matrixPiece[ (ownData.rowNum-1)*params.gridSize + i] = matrix[ (params.gridSize-1)*params.gridSize + i];
		}
	}
	for ( UINT i = 0; i < (ownData.rowNum-2)*params.gridSize; i++)
	{
		ownData.derivative[i] = ownData.matrixPiece[i];
	}
	return ownData;
}

void ExchangeBoundRows(ProcessParameters process, 
					   CellPoint * rows, 
					   UINT height, 
					   UINT gridSize)
{
	MPI_Status status;
	int nextProcNum = (process.rank == process.num - 1)? MPI_PROC_NULL:process.rank + 1;
	int prevProcNum = (process.rank == MAIN_THREAD_RANK)?MPI_PROC_NULL:process.rank - 1;

	MPI_Sendrecv( rows + gridSize*(height - 2), 
				  gridSize,
				  MPI_CELLPOINT, 
				  nextProcNum, 
				  4, 
				  rows,
				  gridSize,
				  MPI_CELLPOINT,
				  prevProcNum,
				  4,
				  MPI_COMM_WORLD,
				  &status
				  );
	MPI_Sendrecv( rows + gridSize, 
				  gridSize,
				  MPI_CELLPOINT, 
				  prevProcNum, 
				  5, 
				  rows + gridSize*(height - 1),
				  gridSize,
				  MPI_CELLPOINT,
				  nextProcNum,
				  5,
				  MPI_COMM_WORLD,
				  &status
				  );
}

CellPoint CalculateIteration(OwnProcessData ownData, AlgorithmParameters params)
{
	CellPoint dm, dmax, temp;
	dmax = (CellPoint)0;
	for ( UINT i = 1; i < ownData.rowNum - 1; i++ )
	{
		for ( UINT j = 1; j < params.gridSize - 1; j++)
		{
			temp = ownData.matrixPiece[params.gridSize*i + j];
			ownData.matrixPiece[params.gridSize*i + j] = (CellPoint)0.25 * (
														   ownData.matrixPiece[ params.gridSize*i + j + 1] +
														   ownData.matrixPiece[ params.gridSize*i + j - 1] +
														   ownData.matrixPiece[ params.gridSize*(i+1) + j] +
														   ownData.matrixPiece[ params.gridSize*(i-1) + j] -
														   //params.stepSize*params.stepSize*ownData.derivative[params.gridSize*i + j]
														   params.stepSize*params.stepSize*f(i,j)
			);
			dm = fabs(ownData.matrixPiece[params.gridSize*i + j] - temp);
			if (dmax < dm) dmax = dm;
		}
	}
	return dmax;
}


void SolvePDEProblem(AlgorithmParameters params, 
					 ProcessParameters process,
					 OwnProcessData ownData)
{
	CellPoint procDelta = 0;
	CellPoint maxDelta = 0;
	int Iterations = 0;
#ifdef _DEBUG
	if (process.rank == MAIN_THREAD_RANK) std::cout<<"SolvePDEProblem:"<>a;
		MPI_Barrier(MPI_COMM_WORLD);
#endif
		MPI_Allreduce(&procDelta, &maxDelta, 1, MPI_CELLPOINT, MPI_MAX, MPI_COMM_WORLD);
	} while (maxDelta > params.precision);
}

Оценка времени выполнения

Для анализа будем использовать следующие показатели:
n - количество узлов по каждой из координат области D.

m - число операций, выполняемых методом для одного узла сетки.

k - количество итераций метода до выполнения условия остановки.

T1 - время решения задачи в 1 процесс.

Tp - время решения задачи запущенной в p процессов.

S - ускорение (speedup). Ускорение определяется из отношения: S=T1/Tp. Числитель в этом выражении: параллельная программа запущеная в 1 процесс.

Вычислительная трудоёмкость последовательного алгоритма T1 = k*m*n*n.

Для p потоков:Tp=k*m*n*n/p.

Из элементарных соображений получаем, что: S -> p

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



где S - время работы последовательного алгоритма, p - число процессов, K(i) - число итераций алгоритма ( зависит от задаваемой точности), альфа - латентность, бетта - пропускная способность, m - порядок матрицы

Результаты экспериментов

Характеристика оборудования:
OS: Windows 7 x64

CPU: AMD Athlon 64 X2 3800+
RAM: 4GB DDR2

Точночть: 0.001
Сторона сетки   P1     P2     P1/P2
100             391    250    1,564
500             4953   2812   1,761
1000            26766  16000  1,672
1500            65438  35906  1,822
2000            123797 67078  1,845

Автор

Блохин Олег, группа 83-03, 2010 год

Новости

22.10.2012
04.09.2012
05.04.2012
06.03.2012
02.03.2012