Размещено на http://www.allbest.ru/

Распараллеливание в OpenMP среднеквадратических приближений неполиномиальными сплайнами минимального дефекта

Винник М.П. mvinnik92@gmail.com 541 группа Математико-механический факультет СПБГУ

профессор Бурова Ирина Герасимовна burovaig@mail.ru каф. Вычислительной математики Математико-механический факультет СПБГУ



Аннотация

В данной работе рассматривается применение неполиномиальных сплайнов минимального дефекта к задаче построения среднеквадратического приближения. Исследуются различные варианты оптимизации решения методом релаксации системы линейных алгебраических уравнений, возникающей в процессе построения среднеквадратического приближения. Проведен сравнительный анализ различных вариантов распараллеливания программ.

среднеквадратический приближение линейный уравнение



Введение

Практически всегда экспериментально полученные данные содержат различные погрешности и помехи. При решении такого рода задач появляется необходимость корректировки данных. Базисные функции второй степени минимального дефекта приводят к решению системы линейных алгебраических уравненний с ленточной матрицей с шириной ленты равной пяти и дают непрерывно дифференцируемое решение задачи среднеквадратического приближения.



1.Построение базисных сплайнов

1.1 О построении полиномиальных В-сплайнов второй степени

Пусть - натуральное число, на промежутке задана равномерная сетка узлов

с шагом

так что

Базисные сплайны, , удовлетворяющие условиям:

,

будем строить решая систему уравнений с параметрами относительно . На промежутке система имеет вид:

Аналогичные системы уравнений имеем на соседних промежутках.

Параметры подбираем так, чтобы .

При значениях параметров

находим:

Нетрудно получить формулу базисного сплайна:

1.2 О построении тригонометрических сплайнов минимального дефекта

Для построения тригонометрических базисных сплайнов, таких, что рассмотрим три системы уравнений с параметрами .

Таким образом, при

на промежутке



и если ---

.

Значения параметров находим из условия: . Таким образом, получаем , , и далее на промежутке находим

Объединяя формулы , найденные на трех соседних промежутках с одинаковым номером , получаем:

1.3 Вычисление элементов матрицы

Как известно, задача среднеквадратического приближения приводит к задаче решения системы линейных уравнений с матрицей Грама.
В частном случае, при , структура матрицы Грама будет иметь следующий вид:

Вычисляя скалярные произведения, в случае полиномиальных B-сплайнов, получаем следующие значения:

Вычислив интегралы, получим:

Элементы вектора правой части таковы:
Вычисляя скалярные произведения в случае тригонометрических сплайнов минимального дефекта, получаем следующие значения:

Элементы вектора в правой части F таковы:



2. Распараллеливание вычислений и оптимизация

2.1 Хранение данных

При решении задач, требующих хранения информации большого объема, встает вопрос о способе хранения данных. В нашем случае матрица системы уравнений имеет ленточный пятидиагональный вид, и ранее было получено, что различных элементов в ней всего 6.

В данной работе было рассмотрено два подхода к хранению данных.

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

Второй подход состоит в том, что реализована специально для этой задачи функция, которая по номеру элемента матрицы определяет его значение. Использование такого подхода гарантирует нам очень большой выигрыш в памяти (вместо элементов мы храним только 6 для любых размеров матрицы), но дает также сильный проигрыш в скорости.

Был проведен сравнительный анализ времени выполнения программ, иллюстрирующих оба подхода (таблица 2).

Теперь рассмотрим функцию вычисления элемента:

get_elem (int i, int j,int n,double h,bool zero_on_diag)

{

int low=min(i,j);

if ((i==j)&(i>=3)&(i<=n)){

if (zero_on_diag) return 0;

return 41/20.0*h;

}

if ((abs(i-j)==1)&(low>=2)&(low<=n)){

//if (zero_on_diag) return 47/60.0*h/(41/20.0*h);

if (zero_on_diag) return 47/60.0*h/get_elem(i,i,n,h,false);

return -47/60.0*h;

}

if ((abs(i-j)==2)&(low>=1)&(low<=n)){

if (zero_on_diag) return -31/120.0*h/get_elem(i,i,n,h,false);

return 31/120.0*h;

}

if ((i==1)&(j==1)){

if (zero_on_diag) return 0.0;

return 31/20.0*h;

}

if ((i==2)&(j==2)){

if (zero_on_diag) return 0.0;

return 2*h;

}

if ((i==n+1)&(j==n+1)){

if (zero_on_diag) return 0;

return h/2;

}

if ((i==n+2)&(j==n+2)){

if (zero_on_diag) return 0;

return 1/20.0*h;

}

if ((i==2)&(j==1)){

if (zero_on_diag)return 77/120.0*h/(2.0*h);;

return -77/120.0*h;

}

if((i==1)&(j==2)){

if (zero_on_diag)return 77/120.0*h/(31/20.0*h);

return -77/120.0*h;

}

if ((i==n+2)&(j==n+1)){

if (zero_on_diag) return 17/120.0*h/(1/20.0*h);

return -17/120.0*h;

}

if((i==n+1)&(j==n+2)){

if (zero_on_diag) return 17/120.0*h/(1/2.0*h);

return -17/120.0*h;

}

if (zero_on_diag) return -0.0;

return 0.0;

}

Это первый вариант функции вычисления элемента матрицы.

Функцию получения элемента матрицы также можно реализовать с помощью конструкции switch, что даст нам небольшое преимущество в скорости. К сожалению, в C++ предусмотрена подстановка в аргументы только константных выражений, поэтому функция выглядит довольно громоздко и в ней есть повторяющиеся участки кода:

double get_elem_case (int i, int j,int n,double h,bool zero_on_diag)

{

int low=min(i,j);

switch(i)

{

case 1:

if ((i==1)&(j==1)){

if (zero_on_diag) return 0.0;

return 31/20.0*h;

}

if ((abs(i-j)==2)&(low>=1)&(low<=n)){

if (zero_on_diag) return -31/120.0*h/get_elem_case(i,i,n,h,false);

return 31/120.0*h;

}

if((i==1)&(j==2)){

if (zero_on_diag)return 77/120.0*h/(31/20.0*h);

return -77/120.0*h;

}

break;

case 2:

if ((i==2)&(j==1)){

if (zero_on_diag)return 77/120.0*h/(2.0*h);;

return -77/120.0*h;

}

if ((i==2)&(j==2)){

if (zero_on_diag) return 0.0;

return 2*h;

}

if ((abs(i-j)==1)&(low>=2)&(low<=n)){

if (zero_on_diag) return 47/60.0*h/get_elem_case(i,i,n,h,false);

return -47/60.0*h;

}

if ((abs(i-j)==2)&(low>=1)&(low<=n)){

if (zero_on_diag) return -31/120.0*h/get_elem_case(i,i,n,h,false);

return 31/120.0*h;

}

break;

default:

if ((i==j)&(i>=3)&(i<=n)){

if (zero_on_diag) return 0;

return 41/20.0*h;

}

if ((abs(i-j)==1)&(low>=2)&(low<=n)){

if (zero_on_diag) return 47/60.0*h/get_elem_case(i,i,n,h,false);

return -47/60.0*h;

}

if ((abs(i-j)==2)&(low>=1)&(low<=n)){

if (zero_on_diag) return -31/120.0*h/get_elem_case(i,i,n,h,false);

return 31/120.0*h;

}

if ((i==n+1)&(j==n+1)){

if (zero_on_diag) return 0;

return h/2;

}

if ((i==n+2)&(j==n+2)){

if (zero_on_diag) return 0;

return 1/20.0*h;

}

if ((i==n+2)&(j==n+1)){

if (zero_on_diag) return 17/120.0*h/(1/20.0*h);

return -17/120.0*h;

}

if((i==n+1)&(j==n+2)){

if (zero_on_diag) return 17/120.0*h/(1/2.0*h);

return -17/120.0*h;

}

break;

}

if (zero_on_diag) return -0.0;

return 0.0;

}

Сравнительный анализ времени выполнения приведен в таблице 3.

2.2 Поиск оптимального параметра

Пусть система приведена к виду . Оптимальное значение параметра находится по формуле:

Для поиска спектрального радиуса, в программе реализован степенной метод:

lambda2=lambda+1;

double temp2=0;

while(abs(lambda-lambda2)>eps)

{

k++;

lambda2=lambda;

for (i=1;i<=n+2;i++)

{

npH2[i]=npH[i];

npH[i]=0;

}

for (i=1;i<=n+2;i++)

{

for (j=1;j<=n+2;j++)

{

npH[i]=npH[i]+

npH2[j]*get_elem(i,j,n,h,false);

}

}

lambda=npH[n/2]/npH2[n/2];

}

omega=2/(1+sqrt(1-lambda*lambda));

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

lambda2=lambda+1;

double temp2=0;

while(abs(lambda-lambda2)>eps)

{

#pragma omp parallel for

for (i2=0;i2<iter_size;i2++){

k++;

lambda2=lambda;

for (i=1;i<=n+2;i++)

{

npH2[i]=npH[i];

npH[i]=0;

}

for (i=1;i<=n+2;i++)

{

for (j=1;j<=n+2;j++)

{

npH[i]=npH[i]+
npH2[j]*get_elem(i,j,n,h,false);

}

}

lambda=npH[n/2]/npH2[n/2];

}

}

omega=2/(1+sqrt(1-lambda*lambda));

Второй подход состоит в распараллеливании только внутреннего цикла, в котором матрица умножается на вектор.

Сравнительный анализ времени выполнения в таблице (5).

2.3 Решение системы методом верхней релаксации

Рассмотрим решение системы уранвнеий , где x - вектор-решение, а

.

Расчетная формула имеет вид:



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

for (i=1;i<=n+2;i++)

{

dtmp2=0;

for (j=1;j<=n+2;j++)

{

dtmp2=dtmp2+get_elem(i,j,n,h,true)*x[j];//x_k+1[i]

}

x[i]=omega*F[i]+omega*dtmp2-x1[i]*(omega-1);

}

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

do{

k2++;

#pragma omp parallel for shared(x,omega,F,n) private(i,j,dtmp2,x1)

for (int i5=0;i5<iter_size;i5++){

//вычисления

}

}

while(dtemp1>pogr);

Второй подход предполагает распараллеливание только внутреннего цикла с тем же набором общих и разделяемых переменных.

do{

k2++;

//x1=x;//запоминаем предыдущую итерацию

for (int i4=1;i4<=n+2;i4++)

{

x1[i4]=x[i4];

}

#pragma omp parallel for shared(x,omega,F,n)

private(i,j,dtmp2,x1)

for (int i=1;i<=n+2;i++)

{

//вычисления

}

}

while(dtemp1>pogr);

На матрицах размера nxn, n<5000, второй подход показывает даже лучшее время, чем первый.

Сравнительный анализ времени выполнения:

	Обычный цикл	Параллельный цикл (I)	Параллельный цикл (II)
100	172	187	141
500	3338	2512	2324
1000	9937	9375	9220
2000	38517	36535	38652
5000	243298	231962	249818



Результаты

Целью работы было произвести оптимизацию времени расчетов. Были произведены замеры времени на различных матрицах. Результаты приведены в сводной таблице. Здесь и далее все измерения проводились при погрешности метода релаксации , равной и при априорном фиксированном количестве итераций (для построения эффективной параллельной программы) равном . Время в таблицах указано в миллисекундах.

1)Сравнительный анализ производительности (римская цифра I указывает на обычный способ хранения матрицы, а II - на оптимизированный)

N(размер матрицы)	Последовательная программа(I)	Последовательная программа(II)	Оптим. программа(I)	Оптим. программа(II)
100	31	156	47	97
250	109	827	78	301
500	436	4415	218	1193
750	568	7878	364	2640
1000	1466	12386	827	4704
2000	2824	43539	2871	18160
5000	18439	264077	18081	133401

2)Сравнительный анализ способов хранения данных: (римская цифра I указывает на обычный способ хранения матрицы, а II - на оптимизированный)

N(размер матрицы)	Последовательная программа(I)	Последовательная программа(II)	Оптим. программа(I)	Оптим. программа(II)
100	31	156	47	172
250	109	827	78	796
500	436	4415	218	2886
750	568	7878	364	6147
1000	1466	12386	827	10967
2000	2824	43539	2871	40544
5000	18439	264077	18081	252008

3)Сравнительный анализ времени доступа к элементу: (общее время доступа ко всем элементам матрицы)

N(размер матрицы)	без switch	switch	встроенное
250	6	3	0
500	17	12	0
750	51	48	0
1000	61	49	0
2000	220	196	15
5000	1404	1355	62
10000	5725	5727	234

4)Сравнительный анализ степенного метода:

	обычный	оптимизированный (1)	оптимизированный (2)
100	94	141	47
500	951	1001	562
1000	2450	2867	1388
2000	5709	8689	4711
5000	19172	49301	18027
10000	52931	174428	51261

Целью приближения функции сплайнами была очистка функции от “помех”, искажающих истинное значение.

Результат среднеквадратического приближения функции

полиномиальными В-сплайнами второй степени при и можно

видеть на рис.1

Рис. 1



Заключение

Благодаря особенностям распараллеливания, в оптимизированных программах реализуется выигрыш в точности вычислений (а благодаря многопоточности, одновременно и выигрыш во времени). Оптимальная функция доступа к элементу матрицы обеспечит нам еще больший выигрыш в скорости. Точная оценка числа итераций в общем случае является достаточно серьезной проблемой, а выигрыш в точности обеспечен на любой матрице. Программа демонстрирует мощь Open MP как средства оптимизации и языка C++, как средства реализации объемных вычислений.



Литература

Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. 1984 .

Рябенький В.С. Введение в вычислительную математику. 2000.

Фаддеев Д.К. Фаддеева В.Н., Вычислительные методы линейной алгебры. 2009.

Бурова И.Г. Евдокимова Т.О., Приближения неполиномиальными сплайнами минимального дефекта. 2007.

Размещено на Allbest.ru