Программирование для многопроцессорных систем в стандарте MPI - Шпаковский Г.И., Серикова Н.В.
.pdfCALL BLACS_GRIDEXIT( ICTXT ) CALL BLACS_EXIT(0)
9998 FORMAT( 2X, A5, ' : ', I6 )
9999 FORMAT( 2X, 60A )
500 continue stop
end
subroutine pmatgen(a,DESCA,np,nq,b,DESCB,nprow,npcol,myrow,mycol) integer n, i, j, DESCA(*), DESCB(*), nprow, npcol, myrow, mycol
double precision a(*), b(*) nb = DESCA(5)
!Заполнение локальных частей матриц A и B,
!матрица A формируется по алгоритму A(I,J) = I, a
!матрица B(I,J) = 1./J здесь имеются в виду глобальные индексы. k = 0
do 250 j = 1,nq jc = (j-1)/nb
jb = mod(j-1,nb)
jn = mycol*nb + jc*npcol*nb + jb + 1 do 240 i = 1,np
ic = (i-1)/nb
ib = mod(i-1,nb)
in = myrow*nb + ic*nprow*nb + ib + 1 k = k + 1
a(k) = dfloat(in)
b(k) = 1.D+0/dfloat(jn) 240 continue
250 continue return
end
13.2. БИБЛИОТЕКА PETSC
Библиотека PETSc [30] облегчает разработку крупномасштабных приложений на языках Fortran, C и C++, является мощным средством для численного решения дифференциальных уравнений в частных производных (ДУЧП) и связанных с этим проблем на быстродействующих компьютерах. PETSc использует MPI стандарт для взаимодействия процессов. Библиотека PETSc включает различные компоненты (подобно классам в С++), каждый компонент имеет дело с частным семейством объектов (например, векторами) и операциями, которые нужно выполнять над этими объектами. Объекты и операции в PETSc определены на основе долгого опыта научных вычислений. Некоторые основные объекты PETSc представлены на рис. 13.1.
251
|
|
|
|
|
|
Прикладные коды |
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
PDE Solvers |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
TS |
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
(Time Stepping) |
|||||
|
|
|
SNES |
|
|
|
|
|
|
|
|
||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
(Nonlinear Equotion Solvers) |
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
SLES |
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
(Linear Equotion Solvers) |
||||||||
|
|
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
RSP |
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
PC |
|
|
|
Draw |
||||||
|
|
(Krolov Subspase Method) |
|
|
(Preconditionrs) |
|
|||||||||||||
|
|
|
|
|
|||||||||||||||
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Matrices |
|
|
Vectors |
|
|
|
Index Sets |
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
BLAS |
|
LAPACK |
|
MPI |
|
|
|
|
|
Рис. 13.1. Организация библиотеки PETSc
Библиотеки BLAS, LAPACK и MPI уже были описаны ранее. Состав векторных и матричных операций будет представлен несколько позже. Индексные ряды (Index Sets) являются обобщением рядов целочисленных индексов, которые используются для определения элементов векторов и матриц при выполнении операций рассылки и сборки. Частным случаем такого ряда является обычный список координат необходимых элементов.
На основе нижеуказанных программ строятся программы более высокого уровня абстракции, которые в свою очередь являются базой следующего уровня программ: программ – решателей.
Для подпространств Крылова (KSP) программно реализованы методы сопряженных градиентов, обобщенных сопряженных градиентов, би-сопряженных градиентов, методы Чебышева, Ричардсона и другие. Программы предобуславливателей (PC) используют методы полной и неполной LU-факторизации, методы Якоби, аддитивный метод Шварца и другие. PETSc содержит небольшой набор средств для рисования простых фигур (Draw). Этот набор не может конкурировать с более совершенными системами графики, но его наличие в PETSc повышает удобство использования последнего.
252
Программы SLES являются сердцем PETSc, поскольку обеспечивают единообразный и эффективный доступ ко всем пакетам решения линейных систем, включая параллельные и последовательные, прямые и итерационные. Программы 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 указывает одиночный процесс.
253
Пользователи, которые хотят использовать процедуры PETSc только для подмножества процессоров внутри большой работы, или если необходимо использовать главный процесс для координации работы подчиненных 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);
Процедуры создания объектов являются коллективными над всеми процессами в коммуникаторе, все процессы в коммуникаторе обязаны обращаться к процедурам создания объектов. Кроме того, если используется последовательность коллективных процедур, они обязаны вызываться в том же самом порядке на каждом процессе.
Опишем некоторые основные процедуры библиотеки, которые дадут представление о работе с новыми структурами данных.
254
Вектора. 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 указывает число,
255
на единицу большее, чем номер последнего принадлежащего процессору элемента вектора.
Основные векторные операции представлены в табл.13.1.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Таблица 13.1 |
|
|
Основные векторные операции |
|||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
№ |
Название функции |
|
|
|
|
|
|
|
|
|
|
|
|
|
Операция |
|
1 |
VecAXPY(Scalar *a,Vec x, Vec y); |
y = y + a x |
|
|||||||||||||
2 |
VecAYPX(Scalar *a,Vec x, Vec y); |
y = x + a y |
|
|||||||||||||
3 |
VecWAXPY(Scalar *a,Vec x,Vec y, Vec w); |
w = a x + y |
|
|||||||||||||
4 |
VecAXPBY(Scalar *a,Scalar *,Vec x,Vec y); |
y = a x +b y |
|
|||||||||||||
5 |
VecScale(Scalar *a, Vec x); |
x = a x |
|
|||||||||||||
6 |
VecDot(Vec x, Vec y, Scalar *r); |
r = x′ y |
|
|||||||||||||
7 |
VecTDot(Vec x, Vec y, Scalar *r); |
r = x′ y |
|
|||||||||||||
8 |
VecNorm(Vec x,NormType type, double *r); |
r = |
|
|
|
x |
|
|
|
type |
|
|||||
|
|
|
|
|
||||||||||||
|
|
|
|
|
|
|
||||||||||
9 |
VecSum(Vec x, Scalar *r); |
r = ∑xi |
|
|||||||||||||
10 |
VecCopy(Vec x, Vec y); |
y = x |
|
|||||||||||||
11 |
VecSwap(Vec x, Vec y); |
y = x while x = y |
|
|||||||||||||
12 |
VecPointwiseMult(Vec x,Vec y, Vec w); |
wi |
= xi yi |
|
||||||||||||
13 |
VecPointwiseDivide(Vec x,Vec y, Vec w); |
wi |
= xi / yi |
|
||||||||||||
14 |
VecMDot(int n,Vec x, Vec *y,Scalar *r); |
r[i] = x′ y[i] |
|
|||||||||||||
15 |
VecMTDot(int n,Vec x, Vec *y,Scalar *r); |
r[i] = x′ y[i] |
|
|||||||||||||
16 |
VecMAXPY(int n, Scalar *a,Vec y, Vec *x); |
y = y + ∑ai x[i] |
|
|||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
i |
|
17 |
VecMax(Vec x, int *idx, double *r); |
r = max xi |
|
|||||||||||||
18 |
VecMin(Vec x, int *idx, double *r); |
r = min xi |
|
|||||||||||||
19 |
VecAbs(Vec x); |
xi = |
|
xi |
|
|
|
|||||||||
|
|
|
||||||||||||||
|
|
|
|
|
||||||||||||
20 |
VecReciprocal(Vec x); |
xi |
=1/ xi |
|
||||||||||||
21 |
VecShift(Scalar *s,Vec x); |
xi |
= s + xi |
|
||||||||||||
Матрицы. Использование матриц и векторов в PETSc сходно. Пользователь может создавать новую параллельную или последовательную матрицу А, которая имеет М глобальных строк и N глобальных столбцов с помощью процедуры:
MatCreate(MPI_Comm comm, int m, int n, int M, int N, Mat *A);
где распределение элементов матрицы по процессам может быть описано на этапе выполнения, или описано для каждого номера процесса через локальные номера строк и столбцов, используя m и n. Номера строк в локальном процессе можно определить командой:
256
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); |
Y = Y + a X |
|
||||||||
2 |
MatMult(Mat A,Vec x, Vec y); |
y = A x |
|
||||||||
3 |
MatMultAdd(Mat A,Vec x, Vec y,Vec z); |
z = y + A x |
|
||||||||
4 |
MatMultTrans(Mat A,Vec x, Vec y); |
y = AT x |
|
||||||||
5 |
MatMultTransAdd(Mat A,Vec x, Vec y,Vec z); |
z = y + AT x |
|
||||||||
6 |
MatNorm(Mat A,NormType type, double *r); |
r = |
|
|
|
A |
|
|
|
type |
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|||||
7 |
MatDiagonalScale(Mat A,Vec l,Vec r); |
A = diag(l) A diag(r) |
|
||||||||
8 |
MatScale(Scalar *a,Mat A); |
A = a A |
|
||||||||
9 |
MatConvert(Mat A,MatType type,Mat *B); |
B = A |
|
||||||||
10 |
MatCopy(Mat A,Mat B,MatStructure); |
B = A |
|
||||||||
11 |
MatGetDiagonal(Mat A,Vec x); |
x = diag(A) |
|
||||||||
12 |
MatTranspose(Mat A,Mat* B); |
B = AT |
|
||||||||
13 |
MatZeroEntries(Mat A); |
A = 0 |
|
||||||||
14 |
MatShift(Scalar *a,Mat Y); |
Y = Y + a I |
|
||||||||
Решение систем линейных уравнений. После создания матрицы и векторов, которые определяют линейную систему 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);
257
Пользователь сначала создает 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, распределенный по процессам.
258
Результирующий вектор 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);
/* заполняем матрицу 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 . */
259
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; |
|
|
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;
260