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

Распараллеливание в 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, второй подход показывает даже лучшее время, чем первый.

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

Результаты

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

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

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

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