Программирование для многопроцессорных систем в стандарте MPI - Шпаковский Г.И., Серикова Н.В.. Программирование для многопроцессорных систем в стандарте MPI -. Организация вычислений в многопроцессорных системах
 Скачать 1.61 Mb.
|
|
PETSC_COMM_WORLD обращением: PetscSetCommWorld(MPI_Comm comm); перед вызовом PetscInitialize, но после обращения к MPI_Init. PetscSetCommWorld в большинстве случаев вызывается однажды на процесс. Все PETSc процедуры возвращают целое значение, указывающее, имела ли место ошибка на протяжении вызова. Код ошибки устанав- ливается в ненулевое значение, если ошибка была обнаружена, в про- тивном случае устанавливается нуль. Все PETSc программы должны вызывать PetscFinalize как их фи- нальное предложение: PetscFinalize(). Эта процедура вызывает MPI_Finalize, если MPIбыл запущен вы- зовом PetscInitialize. Если MPIбыл запущен внешним по отношению к PETSc образом, например, пользователем, то пользователь несет от- ветственность за вызов MPI_Finalize. Пользователь обязан описывать коммуникатор при создании лю- бого объекта PETSc (вектора, матрицы, метода решения), чтобы ука- зать процессы, на которые распределяется объект. Например, некото- рые команды для создания матриц, векторов и методов решения сис- тем линейных уравнений таковы: MatCreate(MPI_Comm comm, int M, int N, Mat *A); VecCreate(MPI_Comm comm, int m, int M, Vec *x); SLESCreate(MPI_Comm comm, SLES *sles); Процедуры создания объектов являются коллективными над всеми процессами в коммуникаторе, все процессы в коммуникаторе обязаны обращаться к процедурам создания объектов. Кроме того, если ис- пользуется последовательность коллективных процедур, они обязаны вызываться в том же самом порядке на каждом процессе. Опишем некоторые основные процедуры библиотеки, которые да- дут представление о работе с новыми структурами данных. 255 Вектора. PETSc обеспечивает два основных векторных типа: по- следовательный и параллельный (основанный на MPI). Для создания последовательного вектора с m компонентами служит команда VecCreateSeq(PETSC_COMM_SELF,int m,Vec *x). Чтобы создать параллельный вектор необходимо указать число компонент, которые распределяются между процессами. Команда VecCreateMPI(MPI_Comm comm, int m, int M,Vec *x); создает вектор, распределенный между всеми процессами коммуника- тора comm, m указывает число компонент вектора в локальном про- цессе, M – общее число компонент вектора. Или локальную или гло- бальную размерность, но не обе одновременно, можно установить равными PETSC_DECIDE, чтобы указать, что PETSc определит их на этапе выполнения. Параллельный или последовательный вектор x с глобальным раз- мером M можно создать так: VecCreate(MPI_Comm comm,int m,int M,Vec *x); VecSetFromOptions(). В этом случае автоматически генерируется тип вектора (последова- тельный или параллельный). Дополнительные векторы того же типа могут быть сформированы с помощью вызова: VecDuplicate(Vec old,Vec *new); Команды VecSet(Scalar *value,Vec x); VecSetValues(Vec x,int n,int *indices, Scalar *values,INSERT_VALUES); позволяют заполнить векторы некоторыми значениями: или устанав- ливают все компоненты вектора равными некоторому скалярному значению, или присваивают различные значение каждому компонен- ту. Scalar определяется как double в C/C++ (или соответственно double precision в языке Фортран) для версий PETSc, которая не была скомпилирована для использования с комплексными числами. Для параллельных векторов, которые распределяются по процес- сорам в соответствии с их номерами, процедура: VecGetOwnershipRange (Vec vec, int *low, int *high) позволяет определить локальную для данного процессора часть век- тора. Здесь аргумент low указывает первый компонент вектора, при- надлежащий данному процессору, а аргумент high указывает число, 256 на единицу большее, чем номер последнего принадлежащего процес- сору элемента вектора. Основные векторные операции представлены в табл.13.1. Таблица 13.1 Основные векторные операции № Название функции Операция 1 VecAXPY(Scalar *a,Vec x, Vec y); x a y y ∗ + = 2 VecAYPX(Scalar *a,Vec x, Vec y); y a x y ∗ + = 3 VecWAXPY(Scalar *a,Vec x,Vec y, Vec w); y x a w + ∗ = 4 VecAXPBY(Scalar *a,Scalar *,Vec x,Vec y); y b x a y ∗ + ∗ = 5 VecScale(Scalar *a, Vec x); x a x ∗ = 6 VecDot(Vec x, Vec y, Scalar *r); y x r ∗ ′ = 7 VecTDot(Vec x, Vec y, Scalar *r); y x r ∗ ′ = 8 VecNorm(Vec x,NormType type, double *r); type x r = 9 VecSum(Vec x, Scalar *r); ∑ = i x r 10 VecCopy(Vec x, Vec y); x y = 11 VecSwap(Vec x, Vec y); y x while x y = = 12 VecPointwiseMult(Vec x,Vec y, Vec w); i i i y x w ∗ = 13 VecPointwiseDivide(Vec x,Vec y, Vec w); i i i y x w / = 14 VecMDot(int n,Vec x, Vec *y,Scalar *r); ] [ ] [ i y x i r ∗ ′ = 15 VecMTDot(int n,Vec x, Vec *y,Scalar *r); ] [ ] [ i y x i r ∗ ′ = 16 VecMAXPY(int n, Scalar *a,Vec y, Vec *x); ∑ ∗ + = i i i x a y y ] [ 17 VecMax(Vec x, int *idx, double *r); i x r max = 18 VecMin(Vec x, int *idx, double *r); i x r min = 19 VecAbs(Vec x); i i x x = 20 VecReciprocal(Vec x); i i x x / 1 = 21 VecShift(Scalar *s,Vec x); i i x s x + = Матрицы. Использование матриц и векторов в PETSc сходно. Пользователь может создавать новую параллельную или последова- тельную матрицу А, которая имеет М глобальных строк и N глобаль- ных столбцов с помощью процедуры: MatCreate(MPI_Comm comm, int m, int n, int M, int N, Mat *A); где распределение элементов матрицы по процессам может быть опи- сано на этапе выполнения, или описано для каждого номера процесса через локальные номера строк и столбцов, используя m и n. Номера строк в локальном процессе можно определить командой: 257 MatGetOwnershipRange(Mat A, int *first_row, int *last_row). Значения элементов матрицы можно установить командой: MatSetValues(Mat A, int m, int *im, int n, int *in, Scalar *values, INSERT_VALUES); После того, как все элементы помещены в матрицу, она должна быть обработана парой команд: MatAssemblyBegin(Mat A,MAT_FINAL_ASSEMBLY); MatAssemblyEnd(Mat A,MAT_FINAL_ASSEMBLY); Основные матричные операции представлены в табл. 13.2. Таблица 13.2 Основные матричные операции № Название функции Операция 1 MatAXPY(Scalar *a,Mat X, Mat Y); X a Y Y ∗ + = 2 MatMult(Mat A,Vec x, Vec y); x A y ∗ = 3 MatMultAdd(Mat A,Vec x, Vec y,Vec z); x A y z ∗ + = 4 MatMultTrans(Mat A,Vec x, Vec y); x A y T ∗ = 5 MatMultTransAdd(Mat A,Vec x, Vec y,Vec z); x A y z T ∗ + = 6 MatNorm(Mat A,NormType type, double *r); type A r = 7 MatDiagonalScale(Mat A,Vec l,Vec r); ) ( ) ( r diag A l diag A ∗ ∗ = 8 MatScale(Scalar *a,Mat A); A a A ∗ = 9 MatConvert(Mat A,MatType type,Mat *B); A B = 10 MatCopy(Mat A,Mat B,MatStructure); A B = 11 MatGetDiagonal(Mat A,Vec x); ) ( A diag x = 12 MatTranspose(Mat A,Mat* B); T A B = 13 MatZeroEntries(Mat A); 0 = A 14 MatShift(Scalar *a,Mat Y); I a Y Y ∗ + = Решение систем линейных уравнений. После создания матрицы и векторов, которые определяют линейную систему Ax=b , пользова- тель может решить эту систему с помощью SLES, применив последо- вательность команд: SLESCreate(MPI_Comm comm,SLES *sles); SLESSetOperators(SLES sles,Mat A,Mat PrecA,MatStructure flag); SLESSetFromOptions(SLES sles); SLESSolve(SLES sles,Vec b,Vec x,int *its); SLESDestroy(SLES sles); 258 Пользователь сначала создает SLESконтекст и устанавливает опе- раторы, связанные с системой, затем устанавливает различные опции для настройки решения, решает линейную систему и, наконец, разру- шает SLES контекст. Команда SLESSetFromOptions() разрешает пользователю настраивать метод линейного решения на этапе выпол- нения, используя определенный набор опций. Используя этот набор, пользователь может не только избрать итеративный метод и предоб- работку, но и предписать устойчивость сходимости, установить раз- личные мониторинговые процедуры и т.д. Решение систем нелинейных уравнений. Большинство ДУЧП задач – нелинейны. PETSc обеспечивает интерфейс для решения не- линейных задач вызовом SNESаналогичноSLES. 13.3. ПРИМЕРЫ Рассмотрим несколько простейших задач с использованием биб- лиотеки PETSc. Задачи иллюстрируют общие принципы создания программ и использования новых структур данных. Все программы должны включать заголовочные файлы для PETSc: #include "petscsles.h", где petscsles.h есть include файл для SLES компонента. Каждая про- грамма PETSc обязана описать include файл, который соответствует наиболее высокому уровню объектов PETSc, необходимых внутри программы; все необходимые include файлы нижнего уровня включа- ются автоматически внутри файлов верхнего уровня. Например, petscsles.h включает petscmat.h (матрицы), petscvec.h (вектора), и petsc.h (основной PETSc файл). Пользователь может вводить управляющие данные во время ис- полнения, используя опции базы данных. Например, команда OptionsGetInt(PETSC_NULL,"-n", &n, &flg); проверяет, обеспечил ли пользователь опцию командной строки, что- бы установить значение n, которое определяет размер задачи. Если это сделано, соответственно устанавливается переменная n, в против- ном случае n не изменяется. Пример 1. Умножение матрицы на вектор Рассмотрим параллельную программу умножения матрицы на век- тор. Умножаем матрицу А, распределенную полосами из строк по процессам, на параллельный вектор b, распределенный по процессам. 259 Результирующий вектор u будет параллельный и распределен по про- цессам. static char help[] = "умножение матрицы на вектор \n\ Входные параметры:\n\ -n #include "petscmat.h" int main(int argc,char **args) { Vec b,u; /* вектора */ Mat A; /* матрица */ int I,J,Istart,Iend,ierr,n = 8; Scalar zz,one = 1.0; PetscInitialize(&argc,&args,(char *)0,help); /* инициализация Petsc */ /* задаем порядок матрицы или по умолчанию */ ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr); /* Создаем параллельную структуру типа Mat для матрицы А. Матрица будет разделена полосами из строк по процессам. Создаем параллельную матрицу, описывая только ее глобальные размеры. При использовании MatCreate() формат распределения элементов матрицы по процессам будет определен на этапе исполнения: PETSC_DECIDE . */ ierr = MatCreate (PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, n,n, &A); CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); /* Определяем, какие строки матрицы размещены локально. Istart,Iend номера строк матрицы А в процессе */ ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); ierr = PetscPrintf(PETSC_COMM_WORLD," %d %d\n",Istart,Iend); CHKERRQ(ierr); /* заполняем матрицу A некоторыми значениями */ for (I=0; I } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* вывод на экран */ ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* создаем структуру Vect для векторов u, b. Создаем параллельный вектор, описывая только его глобальные размеры. При использовании VecCreate() формат распределения элементов вектора по процессам будет определен на этапе исполнения: PETSC_DECIDE . */ 260 ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,&u); CHKERRQ(ierr); ierr = VecSetFromOptions(u);CHKERRQ(ierr); ierr = VecDuplicate(u,&b);CHKERRQ(ierr); ierr = VecSet(&one,b);CHKERRQ(ierr); /* заполняем вектор b единицами */ /* вывод на экран вектора b */ ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* параллельное умножение */ ierr = MatMult(A,b,u);CHKERRQ(ierr); /* вывод на экран результирующего вектора u */ ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); ierr = PetscFinalize();CHKERRQ(ierr); return 0; } Пример 2. Умножение матриц Обобщим предыдущую программу для умножения матрицы на матрицу. Умножаем параллельную матрицу А, распределенную по- лосами строк по процессам на последовательную матрицу В в каж- дом процессе. Результирующая матрица С будет распределена поло- сами столбцов по процессам. static char help[] = "умножение матрицы на матрицу.\n\ Входные параметры:\n\ -n #include "petscmat.h" int main(int argc,char **args) { Vec b,u; /* вектора */ Mat A; /* матрица */ int I,J,ierr,n = 8,*cc,size; Scalar zz,*ci,**B,**C; PetscInitialize(&argc,&args,(char *)0,help); /* инициализация Petsc */ /* задаем порядок матрицы или по умолчанию */ ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr); /* определяем количество процессов приложения */ ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); /* выделяем память для матрицы В */ ierr = PetscMalloc(n*sizeof(double *),&B);CHKERRQ(ierr); for (I=0; I 261 for (J=0; J } /* выделяем память для матрицы С */ ierr = PetscMalloc(n*sizeof(double *),&C);CHKERRQ(ierr); for (I=0; I /* заполняем матрицу A некоторыми значениями */ for (I=0; I } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* вывод на экран */ ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* создаем структуру Vect для столбцов матрицы В и строк матрицы C. Создаем параллельный вектор, описывая только его глобальные размеры. При использова- нии VecCreate() формат распределения элементов вектора по процессам будет оп- ределен на этапе исполнения: PETSC_DECIDE . */ ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,&u); CHKERRQ(ierr); ierr = VecSetFromOptions(u);CHKERRQ(ierr); ierr = VecDuplicate(u,&b);CHKERRQ(ierr); /* заполняем массив порядка индексов */ ierr = PetscMalloc(n*sizeof(int),&cc);CHKERRQ(ierr); for (J=0; J 262 /*основной цикл умножения матрицы А на строку-вектор матрицы В */ for (I=0; I /* параллельное умножение */ ierr = MatMult(A,u,b);CHKERRQ(ierr); /* заполнение результатов в матрицу С */ ierr = VecGetArray(b,&ci);CHKERRQ(ierr); for (J=0; J /* вывод результирующей матрицы С */ for (I=0; I } ierr = PetscFinalize();CHKERRQ(ierr); return 0; } Пример 3. Умножение матриц (вариант с рассылкой) Если необходимо иметь результирующую матрицу во всех про- цессах, можно использовать механизм scatter/gather. Рассмотрим, как это организовано в библиотеке PETSc. Умножаем параллельную мат- рицу А , распределенную полосами строк по процессам на последова- тельную матрицу В в каждом процессе. Результирующая последова- тельная матрица С будет получена всеми процессами. static char help[] = "умножение матрицы на матрицу.\n\ Входные параметры:\n\ -n #include "petscmat.h" int main(int argc,char **args) { Vec b,u,c; /* вектора */ Mat A; /* матрица */ int I,J,ierr,n = 8,*cc; Scalar zz,**B,*ci,**C; IS is; VecScatter ctx; /* описание переменной для организации рассылки */ PetscInitialize(&argc,&args,(char *)0,help); /* инициализация Petsc */ /* задаем порядок матрицы или по умолчанию */ 263 ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL); CHKERRQ(ierr); /* выделяем память для матрицы В */ ierr = PetscMalloc(n*sizeof(double *),&B);CHKERRQ(ierr); for (I=0; I /* вывод на экран матрицы В */ for (J=0; J } /* Создаем параллельную структуру типа Mat для матрицы А. Матрица будет разделена полосами из строк по процессам. Создаем параллельную матрицу, описывая только ее глобальные размеры. При использовании MatCreate() формат распределения элементов матрицы по процессам будет определен на этапе ис- полнения: PETSC_DECIDE . */ ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE, n,n,&A); CHKERRQ(ierr); ierr = MatSetFromOptions(A);CHKERRQ(ierr); /* заполняем матрицу A некоторыми значениями */ for (I=0; I CHKERRQ(ierr); } ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); /* вывод на экран */ ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); /* создаем структуру Vect для столбцов матрицы В и строк матрицы C. Создаем параллельный вектор, описывая только его глобальные размеры. При использова- нии VecCreate() формат распределения элементов вектора по процессам будет определен на этапе исполнения: PETSC_DECIDE . */ 264 ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,&u); CHKERRQ(ierr); ierr = VecSetFromOptions(u);CHKERRQ(ierr); ierr = VecDuplicate(u,&b);CHKERRQ(ierr); /* создаем структуры для организации рассылки */ ierr = VecCreateSeq(PETSC_COMM_SELF,n,&c);CHKERRQ(ierr); ierr = ISCreateStride(PETSC_COMM_SELF,n,0,1,&is);CHKERRQ(ierr); ierr = VecScatterCreate(b,is,c,is,&ctx);CHKERRQ(ierr); /* заполняем массив порядка индексов */ ierr = PetscMalloc(n*sizeof(int),&cc);CHKERRQ(ierr); for (J=0; J /* параллельное умножение */ ierr = MatMult(A,u,b);CHKERRQ(ierr); /* рассылка вектора всем процессам */ ierr = VecScatterBegin(b,c,INSERT_VALUES,SCATTER_FORWARD,ctx); CHKERRQ(ierr); ierr = VecScatterEnd(b,c,INSERT_VALUES,SCATTER_FORWARD,ctx); CHKERRQ(ierr); /* заполнение результатов в матрицу С */ ierr = VecGetArray(c,&ci);CHKERRQ(ierr); for (J=0; J } /* вывод результирующей матрицы С */ for (I=0; I } ierr = PetscFinalize();CHKERRQ(ierr); return 0; } |