Учебное пособие саранск издательство свмо 2013 2 удк 004. 42 Ббк з97 Авторский знак о753
 Скачать 6.58 Mb.
|
|
Программа. Параллельная программа умножения матрицы на вектор #include #include #include "mpi.h" using namespace std; int main(int argc, char *argv[]) { int numprocs,myid; int *a,*b,*c; int n,p,l; int N=10, M=10; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD,&myid); MPI_Status status; if (myid==0) { a=new int[N*M]; b=new int[M]; c=new int[N]; for (int i=0;i b[i]=1;} for (int j=0;j 33 a[j]=1; } n=int(N/numprocs); for (int k = 1; k MPI_Send(&n, 1, MPI_INT, k,0,MPI_COMM_WORLD); MPI_Send(&a[k*n*M],n*M, MPI_INT, k,1,MPI_COMM_WORLD); MPI_Send(&b[0],M, MPI_INT, k,2,MPI_COMM_WORLD); } } else { MPI_Recv(&n, 1, MPI_INT,0,0,MPI_COMM_WORLD,&status); a=new int[n*M]; b=new int[M]; c=new int[n]; MPI_Recv(a, n*M, MPI_INT,0,1,MPI_COMM_WORLD,&status); MPI_Recv(b, M, MPI_INT,0,2,MPI_COMM_WORLD,&status); } MPI_Barrier(MPI_COMM_WORLD); for (int i=0;i c[i]=0; } MPI_Barrier(MPI_COMM_WORLD); for (int i=0;i for (int j=0;j c[i]+=a[j+i*n]*b[j]; } } MPI_Barrier(MPI_COMM_WORLD); if (myid==0) { for (int i=1; i MPI_Recv(&c[i*n], n, MPI_INT,i,3,MPI_COMM_WORLD,&status); } } else { MPI_Send(c, n, MPI_INT, 0,3,MPI_COMM_WORLD); } if (myid==0) { FILE *f=fopen("result.txt", "w"); for (int i=0; i fprintf(f,"%d\n",c[i]); } fclose(f); 34 } return 0; } 2.7.3 Алгоритм умножения матриц Умножение матрицы размера и матрицы размера приводит к получению матрицы размера , каждый элемент которой определяется в соответствии с выражением: Как следует из выражения, каждый элемент результирующей матрицы есть скалярное произведение соответствующих строки матрицы и столбца матрицы . Этот алгоритм предполагает выполнение операций умножения и столько же операций сложения элементов исходных матриц. Будем предполагать далее, что все матрицы являются квадратными и имеют размер Последовательный алгоритм умножения матриц представляется тремя вложенными циклами: double MatrixA[Size][Size]; double MatrixB[Size][Size]; double MatrixC[Size][Size]; int i,j,k; for (i=0; i for (j=0; j MatrixC[i][j] = 0; for (k=0; k MatrixC[i][j] = MatrixC[i][j] + MatrixA[i][k]*MatrixB[k][j]; } } } Этот алгоритм является итеративным и ориентирован на последовательное вычисление строк матрицы . Действительно, при выполнении одной итерации внешнего цикла (цикла по переменной ) вычисляется одна строка результирующей матрицы (рисунок 2.7). 35 Рис. 2.7. Вычисление строки матрицы На первой итерации цикла по переменной используется первая строка матрицы и все столбцы матрицы для того, чтобы вычислить элементы первой строки результирующей матрицы . Рассмотрим ленточную схему разделения данных. Из определения операции матричного умножения следует, что вычисление всех элементов матрицы может быть выполнено независимо друг от друга. Как результат, возможный подход для организации параллельных вычислений состоит в использовании в качестве базовой подзадачи процедуры определения одного элемента результирующей матрицы . Для проведения всех необходимых вычислений каждая подзадача должна содержать по одной строке матрицы и одному столбцу матрицы . Общее количество получаемых при таком подходе подзадач оказывается равным (по числу элементов матрицы ). Рассмотрев предложенный подход, можно отметить, что достигнутый уровень параллелизма является в большинстве случаев избыточным. Обычно при проведении практических расчетов такое количество сформированных подзадач превышает число имеющихся процессоров и делает неизбежным этап укрупнения базовых задач. В этом плане может оказаться полезной агрегация вычислений уже на шаге выделения базовых подзадач. Возможное решение может состоять в объединении в рамках одной подзадачи всех вычислений, связанных не с одним, а с несколькими элементами результирующей матрицы . Для дальнейшего рассмотрения определим базовую задачу как процедуру вычисления всех элементов одной из строк матрицы . Такой подход приводит к снижению общего количества подзадач до величины . Для выполнения всех необходимых вычислений базовой подзадаче должны быть доступны одна из строк матрицы и все столбцы матрицы . Простое решение этой проблемы – дублирование матрицы во всех подзадачах – является, как правило, неприемлемым в силу больших затрат памяти для хранения данных. Поэтому организация вычислений должна быть построена таким образом, чтобы в каждый текущий момент времени подзадачи содержали лишь часть данных, необходимых для проведения расчетов, а доступ к остальной части данных обеспечивался бы при помощи передачи данных между процессорами. 36 Для вычисления одной строки матрицы необходимо, чтобы в каждой подзадаче содержалась строка матрицы , и был обеспечен доступ ко всем столбцам матрицы . Алгоритм представляет собой итерационную процедуру, количество итераций которой совпадает с числом подзадач. На каждой итерации алгоритма каждая подзадача содержит по одной строке матрицы и одному столбцу матрицы . При выполнении итерации проводится скалярное умножение содержащихся в подзадачах строк и столбцов, что приводит к получению соответствующих элементов результирующей матрицы . По завершении вычислений в конце каждой итерации столбцы матрицы должны быть переданы между подзадачами с тем, чтобы в каждой подзадаче оказались новые столбцы матрицы и могли быть вычислены новые элементы матрицы . При этом данная передача столбцов между подзадачами должна быть организована таким образом, чтобы после завершения итераций алгоритма в каждой подзадаче последовательно оказались все столбцы матрицы . Возможная простая схема организации необходимой последовательности передач столбцов матрицы между подзадачами состоит в представлении топологии информационных связей подзадач в виде кольцевой структуры. В этом случае на каждой итерации подзадача , , будет передавать свой столбец матрицы подзадаче с номером (в соответствии с кольцевой структурой подзадача передает свои данные подзадаче с номером ) – рисунок 2.8. После выполнения всех итераций алгоритма необходимое условие будет обеспечено – в каждой подзадаче поочередно окажутся все столбцы матрицы . На рисунке 2.8 представлены итерации алгоритма матричного умножения для случая, когда матрицы состоят из четырех строк и четырех столбцов . В начале вычислений в каждой подзадаче , , располагаются -я строка матрицы и -й столбец матрицы . В результате их перемножения подзадача получает элемент результирующей матрицы . Далее подзадачи осуществляют обмен столбцами, в ходе которого каждая подзадача передает свой столбец матрицы следующей подзадаче в соответствии с кольцевой структурой информационных взаимодействий. Далее выполнение описанных действий повторяется до завершения всех итераций параллельного алгоритма [4,7]. 37 Рис. 2.8. Общая схема передачи данных для первого параллельного алгоритма матричного умножения при ленточной схеме разделения данных Программа. Параллельная программа умножения матрицы на матрицу #include #include #include #include "mpi.h" #include const int N=4; int main(int argc, char *argv[]) { int r,q,myid,numprocs; int i0; int *b,*c,*loc_a,*loc_c; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&numprocs); MPI_Comm_rank(MPI_COMM_WORLD,&myid); MPI_Status status; q=N/numprocs; b=new int [N*N]; c=new int [N*N]; loc_c=new int[N*N]; loc_a=new int[q]; for(int i=0;i { c[i]=0; loc_c[i]=0; } if(myid==0) { 38 for (int j=0;j { for(r=0;r { loc_a[r]=1; } MPI_Send(&loc_a[0],q*N,MPI_INT,j,0,MPI_COMM_WORLD); } for(int i=0;i { b[i]=1; } } MPI_Recv(&loc_a[0],q*N,MPI_INT,0,0,MPI_COMM_WORLD,&status); for(r=0;r { MPI_Bcast(&b[r*N],N,MPI_INT,0,MPI_COMM_WORLD); i0=myid*q; for(int i=0;i { for(int j=0;j { loc_c[r*N+i0]+=loc_a[i*N+j]*b[r*N+j]; } i0++; } MPI_Reduce(loc_c,c,N*N,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); } if(myid==0) { FILE *f=fopen("result.txt", "w"); for(int i=0;i { for(int j=0;j { fprintf(f,"%d\t",c[j*N+i]); } fprintf(f,"\n"); } fclose(f); } 39 return 0; MPI_Finalize();} 2.7.4 Алгоритм решения СЛАУ методом Гаусса Метод Гаусса – широко известный прямой алгоритм решения систем линейных уравнений, для которых матрицы коэффициентов являются плотными. Если система линейных уравнений невырожденна, то метод Гаусса гарантирует нахождение решения с погрешностью, определяемой точностью машинных вычислений. Основная идея метода состоит в приведении матрицы A посредством эквивалентных преобразований к треугольному виду, после чего значения искомых неизвестных могут быть получены непосредственно в явном виде. Метод Гаусса включает последовательное выполнение двух этапов. На первом этапе – прямой ход метода Гаусса – исходная система линейных уравнений при помощи последовательного исключения неизвестных приводится к верхнему треугольному виду c Ux На обратном ходе метода Гаусса (второй этап алгоритма) осуществляется определение значений неизвестных. Из последнего уравнения преобразованной системы может быть вычислено значение переменной 1 n x , после этого из предпоследнего уравнения становится возможным определение переменной 2 n x и т.д. При внимательном рассмотрении метода Гаусса можно заметить, что все вычисления сводятся к однотипным вычислительным операциям над строками матрицы коэффициентов системы линейных уравнений. Как результат, в основу параллельной реализации алгоритма Гаусса может быть положен принцип распараллеливания по данным. В качестве базовой подзадачи можно принять тогда все вычисления, связанные с обработкой одной строки матрицы и соответствующего элемента вектора . Рассмотрим общую схему параллельных вычислений и возникающие при этом информационные зависимости между базовыми подзадачами. Для выполнения прямого хода метода Гаусса необходимо осуществить итерацию по исключению неизвестных для преобразования матрицы коэффициентов к верхнему треугольному виду. Выполнение итерации , , прямого хода метода Гаусса включает ряд последовательных действий. Прежде всего, в самом начале итерации необходимо выбрать ведущую строку, которая при использовании метода главных элементов определяется поиском строки с 40 наибольшим по абсолютной величине значением среди элементов столбца , соответствующего исключаемой переменной . Поскольку строки матрицы распределены по подзадачам, для поиска максимального значения подзадачи с номерами , , должны обменяться своими элементами при исключаемой переменной . После сбора всех необходимых данных в каждой подзадаче может быть определено, какая из подзадач содержит ведущую строку и какое значение является ведущим элементом. Далее для продолжения вычислений ведущая подзадача должна разослать свою строку матрицы и соответствующий элемент вектора всем остальным подзадачам с номерами , . Получив ведущую строку, подзадачи выполняют вычитание строк, обеспечивая тем самым исключение соответствующей неизвестной . При выполнении обратного хода метода Гаусса подзадачи выполняют необходимые вычисления для нахождения значения неизвестных. Как только какая-либо подзадача , , определяет значение своей переменной , это значение должно быть разослано всем подзадачам с номерами . Далее подзадачи подставляют полученное значение новой неизвестной и выполняют корректировку значений для элементов вектора . Выделенные базовые подзадачи характеризуются одинаковой вычислительной трудоемкостью и сбалансированным объемом передаваемых данных. В случае, когда размер матрицы, описывающей систему линейных уравнений, оказывается большим, чем число доступных процессоров (т.е. ), базовые подзадачи можно укрупнить, объединив в рамках одной подзадачи несколько строк матрицы. Возможно использовать ленточную циклическую схему для распределения данных между укрупненными подзадачами. В этом случае матрица делится на наборы (полосы) строк вида (рисунок 2.9): . Рис. 2.9. Пример использования ленточной циклической схемы разделения строк матрицы между тремя процессорами 41 Сопоставив схему разделения данных и порядок выполнения вычислений в методе Гаусса, можно отметить, что использование циклического способа формирования полос позволяет обеспечить лучшую балансировку вычислительной нагрузки между подзадачами. Распределение подзадач между процессорами должно учитывать характер выполняемых в методе Гаусса коммуникационных операций. Основным видом информационного взаимодействия подзадач является операция передачи данных от одного процессора всем процессорам вычислительной системы. Как результат, для эффективной реализации требуемых информационных взаимодействий между базовыми подзадачами топология сети передачи данных должны иметь структуру гиперкуба или полного графа [8,9]. 2.7.5 Алгоритм параллельной пузырьковой сортировки Последовательный алгоритм пузырьковой сортировки сравнивает и обменивает соседние элементы в последовательности, которую нужно отсортировать. Для последовательности ) ,..., 2 , 1 ( n a a a алгоритм сначала выполняет 1 n базовых операций "сравнения-обмена" для последовательных пар элементов ) , 1 ( )..., 3 , 2 ( ), 2 , 1 ( n a n a a a a a В результате после первой итерации алгоритма самый большой элемент перемещается ("всплывает") в конец последовательности. Далее последний элемент в преобразованной последовательности может быть исключен из рассмотрения, и описанная выше процедура применяется к оставшейся части последовательности ) ' 1 ,..., ' 2 , ' 1 ( n a a a Как можно увидеть, последовательность будет отсортирована после 1 n итерации. // Последовательный алгоритм пузырьковой сортировки void BubbleSort(double A[], int n) { for (int i = 0; i < n - 1; i++) for (int j = 0; j < n - i; j++) compare_exchange(A[j], A[j + 1]); } Алгоритм пузырьковой сортировки в прямом виде достаточно сложен для распараллеливания – сравнение пар значений 42 упорядочиваемого набора данных происходит строго последовательно. В связи с этим для параллельного применения используется не сам этот алгоритм, а его модификация, известная в литературе как метод чет- нечетной перестановки. Суть модификации состоит в том, что в алгоритм сортировки вводятся два разных правила выполнения итераций метода: в зависимости от четности или нечетности номера итерации сортировки для обработки выбираются элементы с четными или нечетными индексами соответственно, сравнение выделяемых значений всегда осуществляется с их правыми соседними элементами. Таким образом, на всех нечетных итерациях сравниваются пары (при четном n), а на четных итерациях обрабатываются элементы После n -кратного повторения итераций сортировки исходный набор данных оказывается упорядоченным. // Последовательный алгоритм чет-нечетной перестановки void OddEvenSort(double A[], int n) { for (int i = 1; i < n; i++) { if (i % 2 == 1) { // нечетная итерация for (int j = 0; j < n/2 - 2; j++) compare_exchange(A[2*j + 1], A[2*j + 2]); if (n % 2 == 1) // сравнение последней пары при нечетном n compare_exchange(A[n - 2], A[n - 1]); } else // четная итерация for (int j = 1; j < n/2 - 1; j++) compare_exchange(A[2*j], A[2*j + 1]); } } Получение параллельного варианта для метода чет-нечетной перестановки уже не представляет каких-либо затруднений – сравнения пар значений на итерациях сортировки являются независимыми и могут быть выполнены параллельно. В случае n p , когда количество процессоров меньше числа упорядочиваемых значений, процессоры содержат блоки данных размера p n / и в качестве базовой подзадачи может быть использована операция "сравнить и разделить" [2,5]. // Параллельный алгоритм чет-нечетной перестановки ParallelOddEvenSort(double A[], int n) { 43 int id = GetProcId(); // номер процесса int np = GetProcNum(); // количество процессов for (int i = 0; i < np; i++ ) { if (i % 2 == 1) { // нечетная итерация if (id % 2 == 1) { // нечетный номер процесса // Cравнение-обмен с процессом — соседом справа if (id < np - 1) compare_split_min(id + 1); } else // Cравнение-обмен с процессом — соседом слева if (id > 0) compare_split_max(id - 1); } else { // четная итерация if(id % 2 == 0) { // четный номер процесса // Cравнение-обмен с процессом — соседом справа if (id < np - 1) compare_split_min(id + 1); } else // Cравнение-обмен с процессом — соседом слева compare_split_max(id - 1); } } } Для пояснения такого параллельного способа сортировки в таблице 2.3 приведен пример упорядочения данных при 4 , 16 p n (т.е. блок значений на каждом процессоре содержит 4 / p n элемента). В первом столбце таблицы приводится номер и тип итерации метода, перечисляются пары процессов, для которых параллельно выполняются операции "сравнить и разделить". Для каждого шага сортировки показано состояние упорядочиваемого набора данных до и после выполнения итерации. Таблица 2.3. Пример сортировки данных параллельным методом чет-нечетной перестановки |