Шпаковский Г.И., Серикова Н.В. Программирование для многопроцессорных систем в стандарте MPI

251

Рис. 13.1. Организация библиотеки PETSc

Библиотеки BLAS, LAPACK и MPI уже были описаны ранее. Со-

став векторных и матричных операций будет представлен несколько

позже. Индексные ряды (Index Sets) являются обобщением рядов це-

лочисленных индексов, которые используются для определения эле-

ментов векторов и матриц при выполнении операций рассылки и

сборки. Частным случаем такого ряда является обычный список коор-

динат необходимых элементов.

На основе нижеуказанных программ строятся программы более

высокого уровня абстракции, которые в свою очередь являются базой

следующего уровня программ: программ – решателей.

Для подпространств Крылова (KSP) программно реализованы ме-

тоды сопряженных градиентов, обобщенных сопряженных градиен-

тов, би-сопряженных градиентов, методы Чебышева, Ричардсона и

другие. Программы предобуславливателей (PC) используют методы

полной и неполной LU-факторизации, методы Якоби, аддитивный

метод Шварца и другие. PETSc содержит небольшой набор средств

для рисования простых фигур (Draw). Этот набор не может конкури-

ровать с более совершенными системами графики, но его наличие в

PETSc повышает удобство использования последнего.

Программы SLES являются сердцем PETSc, поскольку обеспечи-

вают единообразный и эффективный доступ ко всем пакетам решения

SLES

(Linear Equotion Solvers)

Matrices Vectors Index Sets

Прикладные коды

PDE Solvers

SNES

(Nonlinear Equotion Solvers)

RSP

(Krolov Subspase Method)

BLAS

LAPACK MPI

TS

(Time Stepping)

PC

(Preconditionrs)

Draw

252

линейных систем, включая параллельные и последовательные, пря-

мые и итерационные. Программы SNES предназначены для решения

нелинейных систем и используют верхние уровни SLES.

Все перечисленные ранее уровни используются для решения сис-

тем ДУЧП (PDE Solvers). Компонент TS (Time Stepping) обеспечива-

ет решение дифференциальных уравнений, возникающих при дискре-

тизации времязависимых PDEs.

Более полную информацию о возможностях можно получить в

“PETSc Users Manual” [30].

Поскольку PETSc использует MPI для всех межпроцессорных об-

менов, пользователь свободен в использовании процедур MPI в своем

приложении. Однако пользователь изолирован от многих деталей пе-

редачи сообщений внутри PETSc, поскольку они невидимы внутри

параллельных объектов, таких как вектора, матрицы и методы реше-

ния систем. Кроме того, PETSc обеспечивает средства, такие как

обобщенные векторные scatters/gathers и распределенные массивы,

чтобы помочь в управлении параллельными данными.

Большинство программ PETSc начинаются обращением:

PetscInitialize(int *argc, char ***argv, char *file, char *help);

которое инициализирует PETSc и MPI. Аргументы argc и argv являют-

ся аргументами командной строки и поставляются во всех C и C++

программах. Файл аргументов в качестве опции указывает на альтер-

нативное имя файла опций, который размещается по умолчанию в

пользовательском директории. Конечный аргумент help есть знаковая

строка, которая будет печататься, если программа выполняется с оп-

цией –help.

PetscInitialize автоматически вызывает MPI_Init, если он не был

инициализирован ранее. В определенных обстоятельствах, когда нуж-

но инициализировать MPI (или инициализация производится некото-

рой другой библиотекой) необходимо сначала вызвать

MPI_Init, за-

тем выполнить обращение к

PetscInitialize. По умолчанию

PetscInitialize добавляет к MPI_COMM_WORLD коммуникатор с

фиксированным именем

PETSC_COMM_WORLD.

В большинстве случаев пользователь может использовать комму-

никатор

PETSC_COMM_WORLD, чтобы указать все процессы в

данном приложении, а

PETSC_COMM_SELF указывает одиночный

процесс.

Пользователи, которые хотят использовать процедуры PETSc

только для подмножества процессоров внутри большой работы, или

253

если необходимо использовать главный процесс для координации ра-

боты подчиненных PETSc

процессов, должны описать альтернатив-

ный коммуникатор для

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);

Процедуры создания объектов являются коллективными над всеми

процессами в коммуникаторе, все процессы в коммуникаторе обязаны

обращаться к процедурам создания объектов. Кроме того, если ис-

пользуется последовательность коллективных процедур, они обязаны

вызываться в том же самом порядке на каждом процессе.

Опишем некоторые основные процедуры библиотеки, которые да-

дут представление о работе с новыми структурами данных.

Вектора. PETSc обеспечивает два основных векторных типа: по-

следовательный и параллельный (основанный на MPI). Для создания

последовательного вектора с m компонентами служит команда

VecCreateSeq(PETSC_COMM_SELF,int m,Vec *x).

254

Чтобы создать параллельный вектор необходимо указать число

компонент, которые распределяются между процессами. Команда

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 указывает число,

на единицу большее, чем номер последнего принадлежащего процес-

сору элемента вектора.

Основные векторные операции представлены в табл.13.1.

Таблица 13.1

Основные векторные операции

255

№ Название функции Операция

1 VecAXPY(Scalar *a,Vec x, Vec y);

x

ayy ∗+=

2 VecAYPX(Scalar *a,Vec x, Vec y);

ya

x

y ∗+=

3 VecWAXPY(Scalar *a,Vec x,Vec y, Vec w);

y

x

aw +∗=

4 VecAXPBY(Scalar *a,Scalar *,Vec x,Vec y);

ybxay ∗+∗=

5 VecScale(Scalar *a, Vec x);

x

a

x

∗=

6 VecDot(Vec x, Vec y, Scalar *r);

yxr ∗

′

=

7 VecTDot(Vec x, Vec y, Scalar *r);

yxr ∗

′

=

8 VecNorm(Vec x,NormType type, double *r);

type

xr =

9 VecSum(Vec x, Scalar *r);

∑

=

i

xr

10 VecCopy(Vec x, Vec y);

x

y =

11 VecSwap(Vec x, Vec y);

yxwhilexy ==

12 VecPointwiseMult(Vec x,Vec y, Vec w);

iii

yxw ∗=

13 VecPointwiseDivide(Vec x,Vec y, Vec w);

iii

yxw /=

14 VecMDot(int n,Vec x, Vec *y,Scalar *r);

][][ iyxir ∗

′

=

15 VecMTDot(int n,Vec x, Vec *y,Scalar *r);

][][ iyxir ∗

′

=

16 VecMAXPY(int n, Scalar *a,Vec y, Vec *x);

∑

∗+=

i

ixayy ][

17 VecMax(Vec x, int *idx, double *r);

i

xr max=

18 VecMin(Vec x, int *idx, double *r);

i

xr min=

19 VecAbs(Vec x);

ii

xx =

20 VecReciprocal(Vec x);

ii

xx /1=

21 VecShift(Scalar *s,Vec x);

ii

xsx +=

Матрицы. Использование матриц и векторов в PETSc сходно.

Пользователь может создавать новую параллельную или последова-

тельную матрицу

А, которая имеет М глобальных строк и N глобаль-

ных столбцов с помощью процедуры:

MatCreate(MPI_Comm comm, int m, int n, int M, int N, Mat *A);

256

где распределение элементов матрицы по процессам может быть опи-

сано на этапе выполнения, или описано для каждого номера процесса

через локальные номера строк и столбцов, используя m и n. Номера

строк в локальном процессе можно определить командой:

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);

XaYY ∗+=

2 MatMult(Mat A,Vec x, Vec y);

xAy ∗=

3 MatMultAdd(Mat A,Vec x, Vec y,Vec z);

xAyz ∗+=

4 MatMultTrans(Mat A,Vec x, Vec y);

xAy

T

∗=

5 MatMultTransAdd(Mat A,Vec x, Vec y,Vec z);

xAyz

T

∗+=

6 MatNorm(Mat A,NormType type, double *r);

type

Ar =

7 MatDiagonalScale(Mat A,Vec l,Vec r);

)()( rdiagAldiagA ∗∗=

8 MatScale(Scalar *a,Mat A);

AaA ∗=

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);

)(Adiagx =

12 MatTranspose(Mat A,Mat* B);

T

A

B

=

13 MatZeroEntries(Mat A);

0=A

257

14 MatShift(Scalar *a,Mat Y);

IaYY ∗+=

Решение систем линейных уравнений. После создания матрицы

и векторов, которые определяют линейную систему

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);

Пользователь сначала создает 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 файл).

258

Пользователь может вводить управляющие данные во время

исполнения, используя опции базы данных. Например, команда

OptionsGetInt(PETSC_NULL,"-n", &n, &flg);

проверяет, обеспечил ли пользователь опцию командной строки, что-

бы установить значение

n, которое определяет размер задачи. Если

это сделано, соответственно устанавливается переменная

n, в против-

ном случае

n не изменяется.

Пример 1. Умножение матрицы на вектор

Рассмотрим параллельную программу умножения матрицы на век-

тор. Умножаем матрицу

А, распределенную полосами из строк по

процессам, на параллельный вектор

b, распределенный по процессам.

Результирующий вектор

u будет параллельный и распределен по про-

цессам.

static char help[] = "умножение матрицы на вектор \n\

Входные параметры:\n\ -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);

259

/* заполняем матрицу A некоторыми значениями */

for (I=0; I<n; I++)

for (J=0; J<n; J++)

{ zz=I*n+J;

ierr = MatSetValues(A,1,&I,1,&J,&zz,INSERT_VALUES);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 для векторов u, b. Создаем параллельный

вектор, описывая только его глобальные размеры. При использовании

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 = 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 <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;

260

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<n; I++)

ierr = PetscMalloc(n*sizeof(double),&B[I]);CHKERRQ(ierr);

/* заполняем матрицу некоторыми значениями по столбцам */

for (I=0; I<n; I++)

for (J=0; J<n; J++)

B[I][J]=I+1;

for (J=0; J<n; J++) /* вывод на экран */

{

for (I=0; I<n; I++)

ierr = PetscPrintf(PETSC_COMM_WORLD," %f \n", B[I][J]);

CHKERRQ(ierr);

ierr = PetscPrintf(PETSC_COMM_WORLD," \n"); CHKERRQ(ierr);

}

/* выделяем память для матрицы С */

ierr = PetscMalloc(n*sizeof(double *),&C);CHKERRQ(ierr);

for (I=0; I<n; I++)

ierr = PetscMalloc(n/size*sizeof(double),&C[I]);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);

/* заполняем матрицу A некоторыми значениями */

for (I=0; I<n; I++)

for (J=0; J<n; J++)

{ zz=I+1;

ierr = MatSetValues(A,1,&I,1,&J,&zz,INSERT_VALUES);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);

Шпаковский Г.И., Серикова Н.В. Программирование для многопроцессорных систем в стандарте MPI

Подождите немного. Документ загружается.