Лаб.раб. ГАУСС. ГАУСС. Лабораторная работа 3 Решение системы линейных алгебраических уравнений методом Гаусса

• Название: Лабораторная работа 3 Решение системы линейных алгебраических уравнений методом Гаусса
• Анкор: Лаб.раб. ГАУСС
• Дата: 06.09.2022
• Размер: 225.38 Kb.
• Имя файла: ГАУСС.pdf
• Тип: Лабораторная работа; #664828

1
Лабораторная работа 3
Решение системы линейных алгебраических
уравнений методом Гаусса
Суть классического метода Гаусса заключается в следующем. Пусть в системе уравнений
)
0
(
1
n
,
1
n
)
0
(
n
1 2
)
0
(
12 1
)
0
(
11
a x
a
,
x a
x a






)
0
(
1
n
,
2
n
)
0
(
n
2 2
)
0
(
22 1
)
0
(
21
a x
a
,
x a
x a






(1)





)
0
(
1
n
,
n n
)
0
(
nn
2
)
0
(
2
n
1
)
0
(
1
n a
x a
,
x a
x a






первый элемент
0
a
)
0
(
11

и назовем его ведущим элементом первой строки. Поделим все элементы этой строки на
)
0
(
11
a и исключим
1
x из всех последующих строк, начиная со второй, путем вычитания первой (преобразованной), умноженной на коэффициент при
1
x в соответствующей строке. Получим
)
1
(
1
n
,
1
n
)
1
(
n
1 2
)
1
(
12 1
a x
a
,
x a
x






)
1
(
1
n
,
2
n
)
1
(
n
2 2
)
1
(
22
a x
a
,
x a
0











)
1
(
1
n
,
n n
)
1
(
nn
2
)
1
(
2
n a
x a
,
x a
0






Если
0
a
)
1
(
22

, то, продолжая аналогичное исключение, приходим к системе уравнений с верхней треугольной матрицей
)
1
(
1
n
,
1
n
)
1
(
n
1 2
)
1
(
12 1
a x
a
,
x a
x






)
2
(
1
n
,
2
n
)
2
(
n
2 2
a x
a
,
x
0






)
3
(
1
n
,
3
n
)
3
(
n
3 3
a x
a
,
x
0 0












)
n
(
1
n
,
n n
a x
,
0 0







Далее из нее в обратном порядке аналогичным образом исключим с помощью последней строки все n
x во всех строках, расположенных выше последней. Продолжая этот процесс, получим на месте исходной матрицы единичную матрицу, а в столбце правых частей будут находиться значения искомых компонент решения i
x :
)
n
(
1
n
,
n n
a x


n
)
1
n
(
n
,
2
n
)
1
n
(
n
,
1
n
1
n x
a a
x












n
)
1
(
n
,
1 2
)
1
(
2
,
1
)
1
(
1
n
,
1 1
x a
x a
a x







2
Процесс приведения исходной системы к системе с треугольной матрицей называется прямым ходом метода Гаусса, а нахождение неизвестных - обратным.
Рассмотрим один из вариантов организация параллельных вычислений при решении системы уравнений y
x
A



методом Гаусса, где A – квадратная матрица размерности М, x

– вектор-столбец неизвестных, y

– вектор-столбец правых частей системы. К матрице А добавляется вектор правых частей и эта расширенная матрица A

















1
n n
nn
2
n
1
n
1
n
,
2 21 22 21 1
n
,
1
n
1 12 11
a a
a a
a a
a a
a a
a a
A







разрезается на N полос равной ширины (для удобства выберем значение N кратным М), каждая из которых загружается в свой процессор. Далее в нулевом (активном) процессоре первая (ведущая) строка матрицы A делится на первый (ведущий) элемент и она передается всем остальным процессорам. Во всех процессорах с помощью этой строки элементарными преобразованиями все элементы первого столбца матрицы А (за исключением первой строки активного процессора) преобразовываются в нулевые. Затем в активном процессоре ведущей становится следующая строка, ведущим – ее первый ненулевой элемент и процесс продолжается до исчерпания строк в активном процессоре. После этого активным становится следующий (первый) процессор и все повторяется снова до тех пор, пока не исчерпаются все ведущие строки во всех процессорах. В результате в расширенной матрица A матрица А приведется к верхней треугольной матрице, распределенной по всем процессорам. На этом прямой ход метода Гаусса заканчивается. Далее активным становится (N-1) – й процессор и аналогичным образом организуется обратный ход, в результате которого матрица А приводится к единичной. На этом этапе каждая ведущая строка состоит только из одного ненулевого элемента и после приведения матрицы А к единичной на месте последнего столбца расширенной матрицы A оказываются значения вектора-столбца искомого решения x

Порядок выполнения работы
1. Разберите приведенную ниже параллельную программу, написанную на MPI, для решения описанной выше задачи. Число строк в матрице А задайте равным 480 (для удобства счета на 2, 4,
8 и 16 процессорах) а значение N равным целому числу от деления М на число процессоров без остатка. В качестве тестового примера задайте значения элементов матрица А произвольным образом, но так, чтобы ее определитель был отличен от нуля. Задайте произвольные значения вектора-столбца x

, тогда значения вектора-столбца y

найдутся из соотношения x
A
y



. Решите прямую задачу y
x
A



и сравните полученные значения вектора-столбца x

(точнее, его части в каждом процессоре) с исходными.
2. Вставьте в программу замеры времени счетной части программы (t_c) в каждом процессоров и времени (t_mpi), которое занимают обмены с соседними процессорами (“накладные расходы”).
Проведите расчеты на двух, четырех, восьми и шестнадцати процессорах и постройте зависимость времени счета от числа процессоров. Постройте в виде таблицы зависимости t_c и t_mpi от числа процессоров и определите “оптимальное” число процессоров, поддерживающее баланс между временными затратами на счет и “накладными расходами”.
3. Разберите приведенную ниже параллельную программу на OpenMP, в которой полагается, что матрица A целиком располагается в памяти узла. Постройте зависимость коэффициента

3 ускорения от числа ядер (потоков) на узле.
/* Решение СЛАУ методом Гаусса. Программа на MPI, распределение
данных - горизонтальными полосами. (Запуск задачи на 8-ми процессах).
*/
#include
#include
/* Каждая ветвь задает размеры своих полос матрицы MA и вектора правой части.
(Предполагаем, что размеры данных делятся без остатка на количество компьютеров.) */
#define M 400
#define N 50
/* Описываем массивы для полос исходной матрицы - MA и вектор V для приема данных.
Для простоты, вектор правой части уравнений присоединяем дополнительным столбцом к матрице коэффициентов. В этом дополнительном столбце и получим результат. */ double MA[N][M+1], V[M+1], MAD, R; int main(int args, char **argv)
{ int size, MyP, i, j,k, m, p; double t0,dt,t1,t2,t3,t4,dt1,dt2,dt3,dt4;
MPI_Status stat;
/* Инициализация библиотеки */
MPI_Init(&args, &argv);
/* Каждая ветвь узнает размер системы */
MPI_Comm_size(MPI_COMM_WORLD, &size);
/* и свой номер (ранг) */
MPI_Comm_rank(MPI_COMM_WORLD, &MyP);
/* Каждая ветвь генерирует свою полосу матрицы A и свой отрезок вектора правой части, который присоединяется дополнительным столбцом к A. Нулевая ветвь генерирует нулевую полосу, первая ветвь - первую полосу и т.д. (По диагонали исходной матрицы - числа =
2,остальные числа = 1).
*/ for(i = 0; i < N; i++) for(j = 0; j < M; j++) if((N*MyP+i) == j) MA[i][j] = 2.0; else
MA[i][j] = 1.0; for(j=0;j { MA[i][M]=0;for(j = 0; j < M; j++)
MA[i][M]+=MA[i][j]*V[j];}
/* Каждая ветвь засекает начало вычислений и производит вычисления */ t0=MPI_Wtime();
/* Прямой ход */
/* Цикл p - цикл по компьютерам. Все ветви, начиная с нулевой,

4 последовательно приводят к диагональному виду свои строки. Ветвь, приводящая свои строки к диагональному виду, назовем активной, строка, с которой производятся вычисления, так же назовем активной. */ for(p = 0; p < size; p++)
{
/* Цикл k - цикл по строкам. (Все ветви "крутят" этот цикл). */ for(k = 0; k < N; k++)
{ if(MyP == p)
{
/* Активная ветвь с номером MyP == p приводит свои строки к диагональному виду.
Активная строка k передается ветвям, с номером большим чем MyP */
MAD = 1.0/MA[k][N*p+k]; for(j = M; j >= N*p+k; j--)
MA[k][j] = MA[k][j] * MAD; t1=MPI_Wtime(); for(m = p+1; m < size; m++)
MPI_Send(&MA[k][0], M+1, MPI_DOUBLE, m, 1,MPI_COMM_WORLD); dt1=MPI_Wtime()-t1; for(i = k+1; i < N; i++)
{ for(j = M; j >= N*p+k; j--)
MA[i][j] = MA[i][j]-MA[i][N*p+k]*MA[k][j];
}
}
/* Работа принимающих ветвей с номерами MyP > p */ else if(MyP > p)
{t2=MPI_Wtime();
MPI_Recv(V, M+1, MPI_DOUBLE, p, 1,MPI_COMM_WORLD, &stat); dt2=MPI_Wtime()-t2; for(i = 0; i < N; i++)
{ for(j = M; j >= N*p+k; j--)
MA[i][j] = MA[i][j]-MA[i][N*p+k]*V[j];
}
}
} /* for k */
} /* for p */
/* Обратный ход */
/* Циклы по p и k аноалогичны, как и при прямом ходе. */ for(p = size-1; p >= 0; p--)
{ for(k = N-1; k >= 0; k--)
{
/* Работа активной ветви */ if(MyP == p)
{t3=MPI_Wtime(); for(m = p-1; m >= 0; m--)
MPI_Send(&MA[k][M], 1, MPI_DOUBLE, m,1,MPI_COMM_WORLD); dt3=MPI_Wtime()-t3; for(i = k-1; i >= 0; i--)
MA[i][M] -= MA[k][M]*MA[i][N*p+k];
}
/* Работа ветвей с номерами MyP < p */ else
{ if(MyP < p)

5
{t4=MPI_Wtime();
MPI_Recv(&R, 1, MPI_DOUBLE, p, 1,MPI_COMM_WORLD, &stat); dt4=MPI_Wtime()-t4; for(i = N-1; i >= 0; i--)
MA[i][M] -= R*MA[i][N*p+k];
}
}
} /* for k */
} /* for p */
/* Все ветви засекают время и печатают его */ dt = MPI_Wtime()-t0;t0=dt1+dt2+dt3+dt4; printf("MyP = %d Time = %13.4e los time=%13.4e\n", MyP, dt,t0);
/* Все ветви печатают, для контроля, свои первые четыре значения корня */ printf("MyP = %d %f %f %f %f\n",MyP,MA[0][M],MA[1][M],MA[2][M],MA[3][M]); dt=(double)N*MyP; printf("toch = %f %f %f %f\n",-(dt+1)/2,-(dt+2)/2,-(dt+3)/2,-(dt+4)/2);
/* Все ветви завершают выполнение */
MPI_Finalize(); return(0);
}
/*
* Решение СЛАУ методом Гаусса. Параллельная программа на OpenMP.
*/
#include
#include
#include
#include
#define M 1600
/* Описываем массивы для расширенной матрицы - MA и вектор V - решение системы. */ double MA[M][M + 1], V[M]; int main(int argc, char **argv)
{ int i, j, p; double a; struct timeval tv_start, tv_end; gettimeofday(&tv_start, NULL);
/* Заполняется матрица системы */ for (i = 0; i < M; i++) for(j = 0; j < M; j++) if (i == j)
MA[i][j] = 2.0; else
MA[i][j] = 1.0;
/* Задается решение системы. */ for (j = 0; j < M; j++)
V[j] = -(double)(j + 1) / 2.;
/* Вычисляется вектор правой части, который записывается в последний столбец расширенной матрицы */

6 for(i = 0; i < M; i++) {
MA[i][M] = 0.0; for(j = 0; j < M; j++)
MA[i][M] += MA[i][j] * V[j];
}
/* Прямой ход */ for(p = 0; p < M; p++) { a = MA[p][p]; for(i = p; i <= M; i++)
MA[p][i] = MA[p][i] / a;
/* Цикл k - цикл по строкам. (Все ветви "крутят" этот цикл). */
#pragma omp parallel for private (i, j, a) for(j = p + 1; j < M; j++) { a = MA[j][p]; for(i = p; i <= M; i++)
MA[j][i] = MA[j][i] - a * MA[p][i];
}
}
/* Обратный ход */ for(p = M - 1; p >= 0; p--) {
#pragma omp parallel for private (i, j) for(j = p - 1; j >=0; j--) { for(i = M; i > j; i--)
MA[j][i] = MA[j][i]- MA[j][p] * MA[p][i];
}
} gettimeofday(&tv_end, NULL);
/* Вычисленные значения решения расположены в последнем столбце расширенной матрицы MA */
/* Выдаются для контроля первые четыре значения корня */ printf("Вычисленные значения: %f %f %f %f\n", MA[0][M], MA[1][M], MA[2][M], MA[3][M]); printf("Точные значения: %f %f %f %f\n", V[0], V[1], V[2], V[3]); p = 0; for (i = 0; i < M && p == 0; i++) if (fabs(MA[i][M] - V[i]) > 1E-6) p = 1; if (p) { printf("Ошибка! Вычисленные значения не совпадают с точными\n");
} printf("Elapsed time = %.6f sec.\n", (tv_end.tv_sec * 1E6 + tv_end.tv_usec - tv_start.tv_sec * 1E6 + tv_start.tv_usec) / 1E6); return 0;
}