Методич.указание для СРО. Решение систем линейных алгебраических уравнений итерационными методами
 Скачать 36.08 Kb.
|
Методическое указание для выполнения СРОпо дисциплине: «Архитектура параллельных вычислений» Тема: Решение систем линейных алгебраических уравнений итерационными методамиЦель - практическое освоение методов решения систем линейных алгебраических уравнений итерационными методами. В лабораторной работе приведены два примера. Один – решение СЛАУ методом простой итерации; второй – решение СЛАУ методом сопряженных градиентов. В этом разделе даны формулы, по которым осуществляются вычисления в обоих методах и последовательные программы решения СЛАУ этими методами. ПРИМЕР 4.1. Решение СЛАУ методом простой итерации Здесь приведены формулы и последовательные программы данного метода. Дана система линейных алгебраических уравнений: a00x0+a01x1+ ... +a0nxn = f0 a10x0+a11x1+ ... +a1nxn = f1 (1) . . . . . . . . . . . . . . . . an0x0+an1x1+ ... +annxn = fn Краткая запись: 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) Этот метод предназначен для симметричных матриц: aij = aji Выбирается начальное приближение x0. r0 = f – Ax0; z0 = r0; αk = (rk-1,rk-1)/(Azk-1,zk-1) - коэффициент (здесь к – это номер итерационного шага); xk = xk-1+ αkzk-1 - вектор решений на к-й текущей итерации; rk = rk-1- αkAzk-1 - вектор невязки на к-й текущей итерации; βk = (rk,rk)/(rk-1,rk-1) - коэффициент; zk = rk+ βkzk-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; /* Задаем начальное значение r0 и z0. */ 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 = (rk-1,rk-1)/(Azk-1,zk-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 */ /* Вычисляем вектор решения: xk = xk-1+ αkzk-1, вектор невязки: rk = rk-1- αkAzk-1 и числитель для βk равный (rk,rk) */ 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; /* Вычисляем вектор спуска: zk = rk+ βkzk-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); } ЗАДАНИЕ Тщательно изучить программы примеров. Скомпилировать и запустить все программы примеров на одном компьютере. |