Введение
Решение дифференциальных уравнений, является неотъемлемой частью анализа многих математических моделей.
Не всегда есть возможность решить их аналитически, а если есть, то это может быть весьма трудоёмкой задачей.
Использование численных методов помогает избавиться от лишней сложности и получить решение для целого класса задач, не решаемых аналитически.
За счёт использования высокопроизводительных вычислений и большего числа процессоров можно достигнуть большей точности/скорости вычислений за одно и то же время, что позволяет наращивать производительность в зависимости от приоритета задачи.
Постановка задачи
Рассмотрим в качестве примера проблему численного решения задачи Дирихле для уравнения Пуассона, определяемую как задачу нахождения функции 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 год
|
