Программирование для многопроцессорных систем в стандарте MPI - Шпаковский Г.И., Серикова Н.В.. Программирование для многопроцессорных систем в стандарте MPI -. Организация вычислений в многопроцессорных системах
 Скачать 1.61 Mb.
|
|
Глава 10. СИСТЕМЫ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ 10.1. МЕТОДЫ РЕШЕНИЯ СЛАУ Основная задача вычислительной алгебры – решение систем ли- нейных алгебраических уравнений (СЛАУ): a 11 x 1 + a 12 x 2 + ... +a 1n x n = b 1 , a 21 x 1 + a 22 x 2 + ... +a 2n x n = b 2 , a n1 x 1 + a n2 x 2 + ... +a nn x n = b n . Предполагается, что матрица А неособенная, т. е. 0 detA ≠ , и реше- ние единственно. 214 Прямые и итерационные методы. Численные методы решения СЛАУ делятся на две большие группы: прямые и итерационные [24,25]. Прямые методы при отсутствии ошибок округления за конеч- ное число арифметических операций позволяют получить точное ре- шение x * . В итерационных методах задается начальное приближение x 0 и строится последовательность {x k } ∗ ⎯ ⎯ ⎯ → ⎯ ∞ → x k , где k – номер ите- рации. В действительности итерационный процесс прекращается, как только x k становится достаточно близким к x * Итерационные методы привлекательнее с точки зрения объема вы- числений и требуемой памяти, когда решаются системы с матрицами высокой размерности. При небольших порядках системы используют прямые методы либо прямые методы в сочетании с итерационными методами. В этой главе рассмотрены два метода решения СЛАУ – Якоби и Гаусса–Зейделя. Метод Якоби рассмотрен в силу того, что вычисле- ния компонент вектора внутри одной итерации независимы друг от друга и параллелизм очевиден, поэтому легко написать MPI про- грамму. Метод Гаусса–Зейделя более эффективен, поскольку требует за- метно меньше итераций. Высокая сходимость к решению в этом ме- тоде достигается за счет выбора матрицы расщепления, лучше ап- проксимирующей матрицу А. Однако в методе Гаусса–Зейделя вы- числение каждой компоненты вектора зависит от компонент вектора, вычисленных на этой же итерации. Поэтому на первый взгляд метод носит последовательный характер. В этой главе предложен парал- лельный алгоритм для модифицированного метода Гаусса–Зейделя. Метод простой итерации (метод Якоби). Пусть требуется ре- шить систему n R b x, ; b Ax ∈ = . Представим + − + + = A D A A , где D – диагональная матрица с диагональными членами матрицы A; − A – часть матрицы A, лежащая ниже центральной диагонали; + A – часть матрицы A, лежащая выше центральной диагонали. Тогда ( ) b x A D A = + + + − , или ( ) b x A A Dx + + + − − = 215 Запишем итерационный метод в виде ( ) K , , , k ; k k 2 1 0 b x A A Dx 1 = + + − = + − + Разрешим его относительно 1 k x + : K 0,1, k ; b D )x A (A D x 1 k 1 1 k = + + − = − + − − + Нетрудно убедиться, что метод Якоби в координатной форме есть не что иное, как разрешение каждого из уравнений системы относи- тельно одной из компонент вектора. Из первого уравнения системы выражается 1 x и его значение принимается за значение x k 1 1 + . Из вто- рого уравнения определяется 2 x и его значение принимается за x k 1 2 + и т. д. Переменные в правой части этих соотношений при этом полага- ются равными их значениям на предыдущей итерации. Метод Гаусса–Зейделя . Пусть решаемая система представлена в виде b x ) A D A ( = + + + − Итерационная схема Гаусса–Зейделя следует из этого представления системы: ,... , , k ; k k 2 1 0 x A b x ) D A ( 1 = − = + + + − , или ,... , , k ; k k k 2 1 0 b x A x A Dx 1 1 = + − − = + + − + Приведем метод Гаусса–Зейделя к стандартному виду: ,... , k ; k k 1 0 b ) D A ( x A ) D A ( x 1 1 1 = + + + − = − − + − − + Стандартная форма метода позволяет выписать его итерационную матрицу и провести над ней очевидные преобразования: A ) D A ( I ) A D A ( ) D A ( A ) D A ( B 1 1 1 − + − − = − − − − + − − = + − + − − = Представим метод Гаусса–Зейделя в координатной форме для системы общего вида: ,... , , k ; ,n i ; a x a x a b x ii n i j k j ij i j k j ij i k i 2 1 0 1 1 1 1 1 1 = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = ∑ ∑ + = − = + + . (10.2) ) 1 10 ( 1 1 1 1 1 1 ,n. i ; n i j xkj aij i j xkj aij bi aii xki = ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∑ + = − ∑ − = − = + 216 Координатная форма метода Гаусса–Зейделя отличается от координатной формы метода Якоби лишь тем, что первая сумма в правой части итерационной формулы содержит компоненты вектора решения не на k-й, а на (k+1)-й итерации. 10.2. ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ РЕШЕНИЯ СЛАУ 10.2.1. Последовательный алгоритм метода простой итерации Для реализации метода простой итерации должны быть заданы матрица СЛАУ А, вектор свободных членов В, начальное приближе- ние вектора Х, точность вычислений ε . Тогда новые значения вектора Х вычисляются по формуле (10.1). Критерий завершения процесса вычислений n i x x x x k i k i k k ≤ ≤ < − = − + + 1 , max 1 1 ε , где x k – приближенное значение решения на k-м шаге численного ме- тода. Ниже приведен текст процедуры Iter_Jacoby для вычисления но- вых значений компонент вектора Xпо формуле 10.1. int ind(int i,int j,int SIZE) /* процедура обращения к элементу матрицы */ { return (i*(SIZE+1)+j); } void Iter_Jacoby(double *A, double *X, double *X_old, int size) /* задана матрица А, начальное приближение вектора Х_old, размерность матрицы size, вычисляем новое значение вектора Х */ { unsigned int i, j; double Sum; for (i = 0; i < size; ++i) { Sum = 0; for (j = 0; j < i; ++j) Sum += A[ind(i,j,size)] * X_old[j]; for (j = i+1; j < size; ++j) Sum += A[ind(i,j,size)] * X_old[j]; X[i]=(A[ind(i,size,size)] – Sum) / A[ind(i,i,size)]; } } Рис. 10.1. Процедура вычисления значений вектора по методу простой итерации 217 Тогда процедура решения СЛАУ методом простой итерации бу- дет следующей. unsigned long int SolveSLAE(double *A, double *X, double Error, int size) /* задана матрица А размерности size+1 (матрица+столбец свободных членов), начальное приближение вектора Х, погрешность вычислений Error */ { double X_old; /* предыдущее значение вектора Х */ int i,Iter = 0; /* число итераций */ double dNorm, dVal; X_old = malloc(sizeof(double) * size); /* выделяем память для X_old */ do{ ++Iter; /*сохраняем предыдущее значение вектора Х */ memcpy(X_old, X, size); /* вычисляем новое значение вектора Х */ Iter_Jacoby(A, X, X_old, size); dNorm = 0; /* считаем норму погрешности */ for (i = 0; i < size; ++i) { dVal = fabs(X[i] – X_old[i]); if (dNorm < dVal) dNorm = dVal; } }while(Error < dNorm); /* цикл до достижения необходимой точности */ free(X_old ); return Iter; } Рис. 10.2. Последовательная программа реализации метода простой итерации В основном цикле процедуры сохраняется значение вектора, вы- численное на предыдущей итерации цикла, новое значение вектора вычисляется в процедуре Iter_Jacoby, затем вычисляется норма по- грешности. Если достигнута необходимая точность вычислений, осу- ществляется выход из цикла,. 10.2.2. Параллельный алгоритм метода простой итерации Следующая система уравнений описывает метод простой итера- ции: X 1 k+1 =f1(x 1 k ,x 2 k ,. . .x n k ) X 2 k+1 =f2(x 1 k ,x 2 k ,. . .x n k ) … X n k+1 =fn(x 1 k ,x 2 k ,. . .x n k ) 218 Вычисления каждой координаты вектора зависят от значений век- тора, вычисленных на предыдущей итерации, и не зависят между со- бой. Поэтому, естественно, для параллельной реализации можно ред- ложить следующий алгоритм. Количество координат вектора, которые вычисляются в каждом процессе, определим следующим образом. size = (MATR_SIZE/numprocs)+((MATR_SIZE % numprocs) > myid ? 1 : 0 ); где size – количество координат вектора, вычисляемых в данном про- цессе, myid – собственный номер процесса, numprocs – количество процессов приложения, MATR_SIZE – размерность вектора. Тогда каждый процесс может вычислять свое количество size но- вых значений координат вектора. После этого процессы могут обме- няться вновь вычисленными значениями, что позволит главному про- цессу произвести оценку точности вычислений. Процедура Iter_Jacoby для вычисления значений вектора методом простой итерации в параллельной версии изменится следующим об- разом: в ней вычисляется size значений вектора Х с номера first. void Iter_Jacoby(double *X_old, int size, int MATR_SIZE, int first) /* задана матрица А, начальное приближение вектора Х_old, размерность матрицы MATR_SIZE, количество вычисляемых элементов вектора в данном процессе size, вычисляем новые значение вектора Х с номера first */ { int i, j; double Sum; for (i = 0; i < size; ++i) { Sum = 0; for (j = 0; j < i+first; ++j) Sum += A[ind(i,j,MATR_SIZE)] * X_old[j]; for (j = i+1+first; j < MATR_SIZE; ++j) Sum += A[ind(i,j,MATR_SIZE)] * X_old[j]; X[i+first]=(A[ind(i,MATR_SIZE,MATR_SIZE)] – Sum) / A[ind(i,i+first, MATR_SIZE)]; } } Рис. 10.3. Процедура вычисления значений вектора по методу простой итерации (параллельная версия) Ниже представлена процедура SolveSLAE, которая реализует ме- тод простой итерации. 219 void SolveSLAE(int MATR_SIZE, int size, double Error) /* задана матрица А, размерность матрицы MATR_SIZE, количество элементов, вычисляемых в данном процессе size, погрешность вычислений Error*/ { double *X_old ; int Iter = 0, i, Result, first; double dNorm = 0, dVal; /* определение номера элемента first вектора Х, с которого будут вычисляться новые значения в данном процессе*/ MPI_Scan(&size, &first,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); first-= size; /* заполнение массива sendcounts значениями size от каждого процесса*/ MPI_Allgather(&size,1,MPI_INT,sendcounts, 1, MPI_INT, MPI_COMM_WORLD); /* заполнение массива displs – расстояний между элементами вектора, распределенных по процессам */ displs[0] = 0; for (i = 0; i X_old = (double *)malloc(sizeof(double) * MATR_SIZE); do { ++Iter; /*сохраняем предыдущее значение вектора Х */ memcpy(X_old, X, MATR_SIZE); /* вычисляем новое значение вектора Х */ Iter_Jacoby(X_old, size, MATR_SIZE, first); /* Рассылаем вектор Х всем процессам */ MPI_Allgatherv(&X[first], size, MPI_DOUBLE, X, sendcounts, displs, MPI_DOUBLE, MPI_COMM_WORLD ); /* Считаем норму */ if (myid == Root) { for (i = 1; i <= MATR_SIZE; ++i) { dVal = fabs(X[i] – X_old[i]); if (dNorm < dVal) dNorm = dVal; } Result = Error < dNorm; } /* рассылаем результат всем процессам */ MPI_Bcast(&Result, 1, MPI_INT, Root, MPI_COMM_WORLD); }while(Result); free(X_old); return ; } Рис. 10.4. Процедура реализации метода простой итерации (параллельная версия) 220 Прежде всего определяется номер элемента вектора first, с которого вычисляются новые значения в каждом процессе. Переменная first будет в общем случае иметь разные значения, так как количество эле- ментов вектора size различно в процессах. Вычислить значение пере- менной удобно вызовом коллективной функции MPI_Scan: MPI_Scan(&size, &first,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); first-= size; Затем в цикле (аналогично последовательному варианту) сохраня- ется значение вектора, вычисленного на предыдущей итерации; вызо- вом процедуры Iter_Jacoby вычисляются новые значения вектора, осуществляется рассылка вычисленных значений вектора от всех – всем: MPI_Allgatherv(&X[first], size, MPI_DOUBLE, X, sendcounts, displs, MPI_DOUBLE, MPI_COMM_WORLD ); Для выполнения коллективной функции MPI_ Allgatherv необхо- димо заполнить массивы sendcounts и displs (заполнение массивов лучше выполнить до цикла). После обмена каждый процесс имеет но- вый вектор Х. В корневом процессе вычисляется норма погрешности и сравнивается с заданной. Условие выхода из цикла Result ≠ 0, по- этому корневой процесс рассылает значение Result всем процессам: MPI_Bcast(&Result, 1, MPI_INT, Root, MPI_COMM_WORLD). Ниже представлена главная программа метода простой итерации. #include "mpi.h" #include #include #include #include /* myid – номер процесса, numprocs – количество процессов, Root – номер корневого процесса, sendcounts, displs – массивы для организации обменов */ int myid,numprocs,Root=0, *sendcounts, *displs; /* AB – исходная матрица метода, может быть задана любым способом X – вектор, содержит начальное приближение */ double *AB, *A,* X; int main(int argc,char *argv[]) /* MATR_SIZE – размерность системы, size – количество строк матрицы в каждом процессе, SIZE – количество элементов матрицы в каждом процессе */ { int i, size, MATR_SIZE, SIZE; double Error; /* точность вычислений */ MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&numprocs); 221 MPI_Comm_rank(MPI_COMM_WORLD,&myid); if (myid=Root) { /* получаем размерность системы MATR_SIZE */ /* выделяем память для матрицы AB */ AB = (double *)malloc(sizeof(double) * MATR_SIZE * (MATR_SIZE + 1)); /* получение матрицы АB размерности MATR_SIZE+1(матрица+столбец свободных членов), задание точности вычислений Error */ } /* рассылка значений всем процессам */ MPI_Bcast(&MATR_SIZE, 1, MPI_INT, Root, MPI_COMM_WORLD); MPI_Bcast(&Error, 1, MPI_DOUBLE, Root, MPI_COMM_WORLD); /* выделяем память для вектора X */ X = (double *)malloc(sizeof(double) * MATR_SIZE); if (myid=Root) { /*получаем начальное значение для вектора Х */ } /* рассылка вектора Х всем процессам */ MPI_Bcast(X, MATR_SIZE, MPI_DOUBLE, Root, MPI_COMM_WORLD); /* определение количества элементов вектора, которые будут вычисляться в каждом процессе (равно количеству строк матрицы, находящихся в данном процессе) */ size = (MATR_SIZE/numprocs)+((MATR_SIZE % numprocs) > myid ? 1 : 0 ); /* выделяем память для матрицы А в каждом процессе*/ A = (double *)malloc(sizeof(double) * (MATR_SIZE+1)*size); displs= (int *)malloc(numprocs*sizeof(int)); sendcounts = (int *)malloc(numprocs*sizeof(int)); /* рассылка частей матрицы по процессам */ SIZE = (MATR_SIZE+1) * size; MPI_Gather(&SIZE,1,MPI_INT,еndcounts,1,MPI_INT,Root, MPI_COMM_WORLD); displs[0] = 0; for (i = 1;i (MATR_SIZE+1) * size, MPI_DOUBLE,Root, MPI_COMM_WORLD ); /* решение СЛАУ методом простой итерации */ SolveSLAE(MATR_SIZE, size, Error); /* освобождение памяти */ free(sendcounts); free(displs); free(AB); free(A); free(X); MPI_Finalize(); return 0; } Рис. 10.5. Главная программа параллельного алгоритма метода простой итерации 222 В корневом процессе происходит задание размерности исходной системы MATR_SIZE, заполняется матрица значений системы AB (матрица + столбец свободных членов), начальное приближение век- тора решения Х, точность вычислений Error. После инициализации MPI, определения количества процессов в приложении nimprocs, соб- ственного номера процесса myid каждый процесс получает от корне- вого процесса заданную размерность исходной матрицы, точность вычислений, начальное приближение вектора Х: MPI_Bcast(&MATR_SIZE, 1, MPI_INT, Root, MPI_COMM_WORLD); MPI_Bcast(&Error, 1, MPI_DOUBLE, Root, MPI_COMM_WORLD); MPI_Bcast(X, MATR_SIZE, MPI_DOUBLE, Root, MPI_COMM_WORLD); После этого каждый процесс определяет количество координат вектора size, которые будут вычисляться в данном процессе (разные процессы, возможно, будут иметь разные значения size ): size = (MATR_SIZE/numprocs)+((MATR_SIZE % numprocs) > myid ? 1 : 0 ); Далее необходимо разделить исходную матрицу AB по процес- сам: по size строк матрицы в каждый процесс: MPI_Scatterv(AB, sendcounts, displs, MPI_DOUBLE, A, (MATR_SIZE+1) * size, MPI_DOUBLE,Root, MPI_COMM_WORLD ); Для выполнения функции MPI_Scatterv необходимо заполнить массивы sendcounts, displs. В первом из них находится количество элементов матрицы SIZE каждого процесса, во втором – расстояния между элементами вектора, распределенного по процессам, причем SIZE=(MATR_SIZE+1)*size, где MATR_SIZE+1 – количество элементов в строке матрицы, size – количество строк матрицы в процессе. Тогда заполнение первого мас- сива sendcounts выполняется вызовом функции: MPI_Gather(&SIZE, 1, MPI_INT, sendcounts, 1, MPI_INT, Root, MPI_COMM_WORLD); Массив displs заполняется следующим образом: displs[0] = 0; for (i = 1;i 223 10.2.3. Параллельный алгоритм метода Гаусса–Зейделя. Отличие метода Гаусса–Зейделя от метода простой итерации за- ключается в том, что новые значения вектора вычисляются не только на основании значений предыдущей итерации, но и с использованием значений уже вычисленных на данной итерации (формула (10.2)). Текст последовательной программы для вычисления новых значений компонент вектора представлен ниже. void GaussZeidel (double *A, double *X, int size) /* задана матрица А, начальное приближение вектора Х, размерность матрицы size, вычисляем новое значение вектора Х */ { unsigned int i, j; double Sum; for (i = 0; i < size; ++i) { Sum = 0; for (j = 0; j < i; ++j) Sum += A[ind(i,j,size)] * X[j]; for (j = i+1; j < size; ++j) Sum += A[ind(i,j,size)] * X[j]; X[i]=(A[ind(i,size,size)] – Sum) / A[ind(i,i,size)]; } } Рис. 10.6. Процедура вычисления значений вектора по методу Гаусса-Зейделя Следующая система уравнений описывает метод Гаусса-Зейделя. X 1 k+1 =f1(x 1 k ,x 2 k ,. . .x n k ) X 2 k+1 =f2(x 1 k+1 ,x 2 k ,. . .x n k ) … X n k+1 =fn(x 1 k+1 ,x 2 k+1 ,. . . x n-1 k+1 ,x n k ) Вычисления каждой координаты вектора зависят от значений, вы- численных на предыдущей итерации, и значений координат вектора вычисленных на данной итерации. Поэтому нельзя реализовывать параллельный алгоритм, аналогичный методу простой итерации: каждый процесс не может начинать вычисления пока, не закончит вы- числения предыдущий процесс. Можно предложить следующий модифицированный метод Гаус- са–Зейделя для параллельной реализации. Разделим вычисления коор- динат вектора по процессам аналогично методу простой итерации. Бу- дем в каждом процессе вычислять свое количество координат вектора по методу Гаусса-Зейделя, используя только вычисленные значения 224 вектора данного процесса. Различие в параллельной реализации по сравнению с методом простой итерации заключается только в проце- дуре вычисления значений вектора (вместо процедуры Iter_Jacoby используем процедуру GaussZeidel). void GaussZeidel(int size, int MATR_SIZE, int first) /* задана матрица А, размерность матрицы MATR_SIZE, количество вычисляемых элементов вектора в данном процессе size, вычисляем новые значения вектора Х с номера first, используя значения вектора Х */ { int i, j; double Sum; for (i = 0; i < size; ++i) { Sum = 0; for (j = 0; j < i+first; ++j) Sum += A[ind(i,j,MATR_SIZE)] * X[j]; for (j = i+1+first; j < MATR_SIZE; ++j) Sum += A[ind(i,j,MATR_SIZE)] * X[j]; X[i+first]=(A[ind(i,MATR_SIZE,MATR_SIZE)] – Sum) / A[ind(i,i+first, MATR_SIZE)]; } } Рис. 10.7. Процедура вычисления значений вектора по методу Гаусса–Зейделя (параллельная версия) КОНТРОЛЬНЫЕ ВОПРОСЫ И ЗАДАНИЯ К ГЛАВЕ 10 Контрольные вопросы к 10.1 1. В чем различие между прямыми и итерационными методами решения СЛАУ? 2. В чем различие между методоми простой итерации и Гаусса–Зейделя для ре- шения СЛАУ? 3. Почему метод Гаусса–Зейделя эффективнее метода простой итерации? Контрольные вопросы к 10.2 1. Сравните распределение работ между процессами в методе простой итерации и в задаче криптографии в предыдущей главе. 2. В чем заключается распараллеливание метода простой итерации? 3. Почему для распределения матрицы системы по процессам нужно использо- вать коллективную функцию MPI_Scatterv, а не MPI_Scatter? 4. Для чего необходимы массивы sendcounts, displs в процедуре SolveSLAE? 5. Объясните суть переменной first в процедуре Iter_Jacoby параллельной реали- зации. Как иначе можно вычислить ее значение? 6. Как изменится параллельная программа метода простой итерации в случае, если не использовать коллективные функции? 7. В чем различие между методом Гаусса–Зейделя и его модификацией, предло- женной в данной главе? |