Основы параллельного программирования

Основы параллельного программирования

Лабораторная работа № 4

Решение систем линейных алгебраических уравнений итерационными методами

Цель - практическое освоение методов решения систем линейных алгебраических уравнений итерационными методами. В лабораторной работе приведены два примера. Один – решение СЛАУ методом простой итерации; второй – решение СЛАУ методом сопряженных градиентов. В этом разделе даны формулы, по которым осуществляются вычисления в обоих методах и последовательные программы решения СЛАУ этими методами.


ПРИМЕР 4.1. Решение СЛАУ методом простой итерации

Здесь приведены формулы и последовательные программы данного метода.

Дана система линейных алгебраических уравнений:

a₀₀x₀+a₀₁x₁+ ... +a_0nx_n = f₀
a₁₀x₀+a₁₁x₁+ ... +a_1nx_n = f_{1 (1)}
. . . . . . . . . . . . . . . .
a_n0x₀+a_n1x₁+ ... +a_nnx_n = f_n

Краткая запись:

Ax = f (2)
Корни этой системы методом простой итерации вычисляются по следующей формуле:
(3)
Здесь верхний индекс обозначает номер итерационного шага, i - номер корня или правой части. τ – итерационный шаг.
Завершение вычислений по условию:


(4)

Замечание.

(5)
Числитель выражения (4) вычисляется по аналогичной формуле.
Ниже приведена программа для последовательного алгоритма:

#include

#include

#define N 100

double A[N][N], F[N], mf, X[N], X1[N], S[N], msf;

main()
{ int i, j;

long int dt;
double t, e;
struct timeval tv1, tv2;
/* Генерация данных. Здесь задается матрица с элементами равными 1,
* по диагонали равными 2. Матрица берется хорошо обусловленной.
* При правой части равной N+1 все корни будут равными 1. */
mf = 0;
for(i = 0; i < N; i++)
{ for(j = 0; j < N; j++)
(i==j)?(A[i][j]=2):(A[i][j]=1);
F[i] = N+1;
/* Сразу вычисляем сумму квадратов элементов вектора F,
* т.е. подкоренное выражение формулы (4). */
mf += F[i]*F[i];
}

/* Задаем шаг t и e. */

t = 0.01;
e = 0.00001;

/* Задаем начальное приближение корней. В Х1 хранятся значения корней

* к+1-й итерации. */
for(i = 0; i < N; i++)
X1[i] = 0.6;

/* Засекаем время начала вычислений. */

gettimeofday(&tv1,NULL);

do
{ for(i = 0; i < N; i++)

/* В Х хранятся значения корней к-й итерации. */
X[i] = X1[i];

for(msf=0,i = 0; i < N; i++)

{ for(S[i] = 0,j = 0; j < N; j++)
S[i] += A[i][j] * X[j];
X1[i] = X[i] - t*(S[i] - F[i]);
/* Вычисляем сумму квадратов элементов невязки, т.е. подкоренное
* выражение числителя формулы (4). */
msf += (S[i] – F[i])*(S[i] – F[i]);
}
}
while(msf/mf > e*e); /* Проверка условия по формуле (3). */

/* Засекаем время конца вычислений. */

gettimeofday(&tv2,NULL);
dt = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;

/* Выводим на экран время вычислений в секундах */

printf(" Time= %d\n",dt);
/* Для контроля выводим 4-е первых корня */
for(i = 0; i < 4; i++)
printf(" %f\n",X[i]);
return(0);
}

Обратите внимание на функцию замера времени gettimeofday(&tv,NULL). Эту функцию можно использовать и для параллельных программ.

ПРИМЕР 2.2. Решение СЛАУ методом сопряженных градиентов

Здесь приведены формулы и последовательные программы данного метода.

Дана система линейных алгебраических уравнений:

Ax = f (6)
Этот метод предназначен для симметричных матриц: a_ij = a_ji
Выбирается начальное приближение x₀. r₀ = f – Ax₀; z₀ = r₀;
α_k = (r_k-1,r_k-1)/(Az_k-1,z_k-1) - коэффициент (здесь к – это номер
итерационного шага);
x_k = x_k-1+ α_kz_k-1 - вектор решений на к-й текущей итерации;
r_k = r_k-1- α_kAz_k-1 - вектор невязки на к-й текущей итерации;
β_k = (r_k,r_k)/(r_k-1,r_k-1) - коэффициент;
z_k = r_k+ β_kz_k-1 - вектор спуска на к-й текущкй итерации (сопряженное
направление).
Выход из итерационного процесса:
(7)

Ниже приведена программа для последовательного алгоритма.

Обратите на то как вычисляются скалярные произведения векторов.

#include

#include
#include

#define M 100

#define E 0.00001

#define T 0.01

double A[M][M], F[M], Xk[M], Zk[M];

double Rk[M], Sz[M], alf, bet, mf;
double Spr, Spr1, Spz;

int main()

{ int i, j, v;
struct timeval tv1, tv2;
long int dt1;

/* Генерация данных. Здесь задается матрица с элементами равными 1,

* по диагонали равными 2. Матрица берется хорошо обусловленной и
* симметричной. При правой части равной M+1 все корни будут равными 1 */
for(mf=0,i = 0; i < M; i++)
{ for(j = 0; j < M; j++)
{ if(i == j)
A[i][j] = 2.0;
else
A[i][j] = 1.0;
}
F[i] = M + 1;
/* Сразу вычисляем сумму квадратов элементов вектора F,
* т.е. подкоренное выражение формулы (7). */
mf += F[i] * F[i];
}

/* Задаем начальное приближение корней. В Х_k хранятся значения корней

* к-й итерации. */
for(i = 0; i < M; i++)
Xk[i] = 0.2;

/* Задаем начальное значение r₀ и z₀. */

for(i = 0; i < M; i++)
{ for(Sz[i]=0,j = 0; j < M; j++)
Sz[i] += A[i][j] * Xk[j];
Rk[i] = F[i] - Sz[i];
Zk[i] = Rk[i];
}

/* Засекаем время начала вычислений. */

gettimeofday(&tv1, NULL);

do
{

/* Вычисляем числитель и знаменатель для коэффициента
* α_k = (r_k-1,r_k-1)/(Az_k-1,z_k-1) */
Spz = 0;
Spr = 0;
for(i = 0; i < M; i++)
{ for(Sz[i]=0, j = 0; j < M; j++)
Sz[i] += A[i][j] * Zk[j];
Spz += Sz[i] * Zk[i];
Spr += Rk[i] * Rk[i];
}
alf = Spr/Spz; /* α_k */
/* Вычисляем вектор решения: x_k = x_k-1+ α_kz_{k-1,
вектор невязки:} r_k = r_k-1- α_kAz_k-1 и числитель для β_k равный (r_k,r_k) */
Spr1 = 0;
for(i = 0; i < M; i++)
{ Xk[i] += alf*Zk[i];
Rk[i] -= alf*Sz[i];
Spr1 += Rk[i]*Rk[i];
}
/* Вычисляем β_k */
bet = Spr1/Spr;

/* Вычисляем вектор спуска: z_k = r_k+ β_kz_k-1 */

for(i = 0; i < M; i++)
Zk[i] = Rk[i] + bet*Zk[i];
}
/* Проверяем условие выхода из итерационного цикла */
while(Spr1/mf > E*E);

/* Засекаем время конца вычислений. Время решения задачи выводим на экран */

gettimeofday(&tv2, (struct timezone*)0);
dt1 = (tv2.tv_sec - tv1.tv_sec) * 1000000 + tv2.tv_usec - tv1.tv_usec;
printf("Time = %d\n",dt1);

/* Для контроля выводим 8 первых корней */

printf(" %f %f %f %f %f %f %f %f\n",Xk[0],Xk[1],Xk[2],Xk[3],Xk[4],Xk[5],
Xk[6], Xk[7]);

return(0);

}
ЗАДАНИЕ
Тщательно изучить программы примеров. Скомпилировать и запустить все программы примеров на одном компьютере.