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

Введение

Дифференциальные уравнения в частных производных представляют собой широко применяемый математический аппарат при разработке моделей в самых разных областях науки и техники. К сожалению, явное решение этих уравнений в аналитическом виде оказывается возможным только в частных простых случаях, и, как результат, возможность анализа математических моделей, построенных на основе дифференциальных уравнений, обеспечивается при помощи приближенных численных методов решения. Объем выполняемых при этом вычислений обычно является значительным и использование высокопроизводительных вычислительных систем является традиционным для данной области вычислительной математики.

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

Рассмотрим в качестве примера проблему численного решения задачи Дирихле для уравнения Пуассона, определяемую как задачу нахождения функции 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 курса Малышев Александр.

Новости

22.10.2012
04.09.2012
05.04.2012
06.03.2012
02.03.2012