Лаб.раб. ГАУСС. ГАУСС. Лабораторная работа 3 Решение системы линейных алгебраических уравнений методом Гаусса
![]()
|
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]+=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; } |