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


и принимающей значения g(x,y) на границе области (f и g являются функциями, задаваемыми при постановке задачи). В качестве области задания D функции u=u(x,y), далее будет использоваться единичный квадрат.

Необходимо создать программный комплекс, реализующий метод Гаусса-Зейделя для решения задачи Дирихле для уравнения Пуассона и выполнить сравнение времени работы алгоритмов.
Описание Алгоритма
Последовательная версия

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

где h - величина обратная к размерности сетки.
Как правило, сетки берутся квадратными, то есть, число узлов по каждой из размерностей одинаковою.
Последовательность сеток получаемая таким алгоритмом равномерно сходится к решению задачи Дирихле.
Параллельная версия
Параллельная схема основана на делении исходной сетки на меньшие, ленточным способом (полосами).

А схема передачи сообщений между процессами может быть представлена на схеме:

Реализация алгоритма
Последовательная версия
void SolveSystem(double* Matrix,double eps,int m,int n)
{
double dmax= 0; //Наибольшее изменение в матирце Matrix
do
{
dmax=0;
for (int i=m+1; i < (n-1)*m - 1; i++)
{
if ( (i%m == 0)||(i%m == m-1) ) continue;
double temp = Matrix[i];
Matrix[i] = 0.25*(Matrix[i-1]+Matrix[i+1]
+Matrix[i-m]+Matrix[i+m]);
double dm = fabs(temp-Matrix[i]);
if ( dmax< dm ) dmax= dm;
}
} while( dmax> eps);
}
Параллельная версия
//Инициализируем MPI
MPI_Init(&argc,&argv);
//Ранг процесса
int ProcRank;
//Число процессов
int ProcNum;
//Определим общее число процессов в рамках нашей задачи
MPI_Comm_size(MPI_COMM_WORLD,&ProcNum);
//Процесс узнает свой ранг в коммуникаторе MPI_COMM_WORLD
MPI_Comm_rank(MPI_COMM_WORLD,&ProcRank);
//Размеры сетки должны храниться на всех процессах
// Причём число строк (n) обязательно кратно числу процессов!
int m=0,n=0;
MPI_Status Status;
//Просто матрица, используется только 0 процессом как начальная
double * Matrix=new double;
//Матрица с значениями функции в узлах сетки - для каждого процесса своя
double * FunctionValues;
//Матрица в которую считаются значения функции
double* Fvalues = new double;
//Узнать точность
double eps = atof(argv[3]);
//Прочитать и распарсить файл c граничными условиями 0 процессом
if (ProcRank == 0)
{
double *left,*top,*right,*bottom;
ParseFileDim(argv[1],&m,&n);
delete Matrix;
//Память под границы
left = new double[n];
right = new double[n];
top = new double[m];
bottom = new double[m];
//И под саму матрицу
Matrix = new double [m * n];
ParseFile(argv[1],left,top,right,bottom);
CreateMatrix(Matrix,left,top,right,bottom,m,n);
}
//Рассылка размеров матрицы
MPI_Bcast(&m,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
//Выделение памяти на значения функции
FunctionValues = new double [m * n / ProcNum];
if ( ProcRank == 0 )
{
delete Fvalues;
Fvalues = new double[m * n];
ParseForValues(argv[2],Fvalues,m,n);
}
MPI_Scatter(Fvalues,m * n / ProcNum,MPI_DOUBLE,FunctionValues,m * n / ProcNum,MPI_DOUBLE,0,MPI_COMM_WORLD);
//Локальная матрица в которой и будут происходить вычисления для всех процессов
double* LocalMatrix = new double[m * n / ProcNum + 2 * m];
//Разослать в LocaMatrix с нужной позиции
MPI_Scatter(Matrix, m * n / ProcNum, MPI_DOUBLE,&LocalMatrix[m],m * n / ProcNum,MPI_DOUBLE,0,MPI_COMM_WORLD);
//dmax - максимальное изменение значения в матрице
//Условие останова цикла основано на том что dmax < eps
double dm,dmax;
//Переменные для передачи сообщения
int NextProc = (ProcRank + 1)%ProcNum;
int PrevProc;
if (ProcRank == 0)PrevProc = ProcNum - 1;
else PrevProc = ProcRank - 1;
//Вычислительный цикл
do
{
dm = 0.0;
MPI_Sendrecv(&LocalMatrix[m * n /ProcNum],m,MPI_DOUBLE,NextProc,0,LocalMatrix,m,MPI_DOUBLE,PrevProc,0,MPI_COMM_WORLD,&Status);
MPI_Sendrecv(&LocalMatrix[m],m,MPI_DOUBLE,PrevProc,0,&LocalMatrix[m * n /ProcNum + m],m,MPI_DOUBLE,NextProc,0,MPI_COMM_WORLD,&Status);
if ( ProcRank == 0 )
{
DoGausseSeidelStep_first(LocalMatrix,FunctionValues,&dm,m,n / ProcNum + 2);
}
if ( ProcRank == ProcNum - 1)
{
DoGausseSeidelStep_last(LocalMatrix,FunctionValues,&dm,m,n / ProcNum + 2);
}
if (( ProcRank != 0) && ( ProcRank != ProcNum - 1))
{
DoGausseSeidelStep(LocalMatrix,FunctionValues,&dm,m, n / ProcNum + 2);
}
MPI_Allreduce(&dm,&dmax,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
}while (dmax > eps);
Оценки сложности
Для анализа будем использовать следующие показатели:
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: Pentium Dual-Core T4400 2.20 GHz
RAM: 4GB DDR3
Латентность: 0.0019358 с
Пропускная способность: 35.08 Мб/с
Разработка велась в MSVS2008 с встроенным компилятором, в режиме релиза была включена полная оптимизация(ключ /Ox)
S1 - последовательная версия, P1 - параллельная в 1 процесс, P2 - параллельная в 2 процесса, P2* - теоретическая оценка времени работы
Debug
|
|
S1
|
P2
|
P2*
|
S1/P2
|
100
|
0,21
|
1,559
|
1,029
|
0,135
|
200
|
0,73
|
1,876
|
1,517
|
0,389
|
400
|
2,915
|
3,324
|
2,768
|
0,877
|
800
|
13,025
|
9,368
|
7,956
|
1,390
|
1600
|
65,073
|
44,151
|
34,140
|
1,474
|
3200
|
277,798
|
160,937
|
138,899
|
1,726
|
Release
|
|
S1
|
P1
|
P2
|
P2*
|
S1/P2
|
100
|
0,033
|
0,042
|
1,135
|
0,940
|
0,029
|
200
|
0,135
|
0,169
|
1,285
|
1,219
|
0,105
|
400
|
0,567
|
0,711
|
1,641
|
1,594
|
0,346
|
800
|
2,393
|
2,97
|
3,332
|
2,640
|
0,718
|
1600
|
18,361
|
20,745
|
15,844
|
10,784
|
1,159
|
3200
|
81,323
|
93,699
|
54,293
|
42,535
|
1,498
|
Автор
Работу выполнил студент 3 курса Малышев Александр.
|