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

261

/* создаем структуру 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<n; J++)

cc[J]=J;

/*основной цикл умножения матрицы А на строку-вектор матрицы В */

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

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

ierr = VecSetValues(u,n,cc,B[I],INSERT_VALUES);CHKERRQ(ierr);

/* параллельное умножение */

ierr = MatMult(A,u,b);CHKERRQ(ierr);

/* заполнение результатов в матрицу С */

ierr = VecGetArray(b,&ci);CHKERRQ(ierr);

for (J=0; J<n/size; J++) C[I][J]=ci[J];

ierr = VecRestoreArray(b,&ci);CHKERRQ(ierr);

}

/* вывод результирующей матрицы С */

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

{

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

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

CHKERRQ(ierr);

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

}

ierr = PetscFinalize();CHKERRQ(ierr);

return 0;

}

Пример 3. Умножение матриц (вариант с рассылкой)

Если необходимо иметь результирующую матрицу во всех про-

цессах, можно использовать механизм scatter/gather. Рассмотрим, как

это организовано в библиотеке PETSc. Умножаем параллельную мат-

рицу

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

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

В в каждом процессе. Результирующая последова-

тельная матрица

С будет получена всеми процессами.

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

262

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

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

ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);

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

/* выделяем память для результирующей матрицы С */

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

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

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

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

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

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

{ B[I][J]=I+1; C[I][J]=0; }

/* вывод на экран матрицы В */

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

}

/* Создаем параллельную структуру типа 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;

263

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 для столбцов матрицы В и строк матрицы 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 = 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<n; J++)

cc[J]=J;

/* основной цикл умножения матрицы А на строку-вектор матрицы В */

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

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

ierr = VecSetValues(u,n,cc,B[I],INSERT_VALUES);CHKERRQ(ierr);

/* параллельное умножение */

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

C[I][J]=ci[J];

ierr = VecRestoreArray(c,&ci);CHKERRQ(ierr);

}

/* вывод результирующей матрицы С */

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

264

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

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

CHKERRQ(ierr);

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

}

ierr = PetscFinalize();CHKERRQ(ierr);

return 0;

}

Пример 4. Решение системы линейных уравнений

Следующий пример иллюстрирует параллельное решение систе-

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

фейс для инициации программы, создания векторов и матриц и реше-

ния линейной системы является тем же для однопроцессорного и мно-

гопроцессорного примеров. Главное различие состоит в том, что каж-

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

раллельном случае.

static char help[] = "Solves a linear system in parallel with SLES.\n\

Input parameters include:\n\

-random_exact_sol : use a random exact solution vector\n\

-view_exact_sol : write exact solution vector to stdout\n\

-m <mesh_x> : number of mesh points in x-direction\n\

-n <mesh_n> : number of mesh points in y-direction\n\n";

#include "petscsles.h"

int main(int argc,char **args)

{

Vec x,b,u; /* приближенное решение, RHS, точное решение */

Mat A; /* матрица линейной системы */

SLES sles; /* метод решения линейной системы */

PetscRandom rctx; /* генератор случайных чисел */

double norm; /* норма ошибки решения */

int i,j,I,J,Istart,Iend,ierr,m = 8,n = 7,its;

PetscTruth flg;

Scalar v,one = 1.0,neg_one = -1.0;

KSP ksp;

PetscInitialize(&argc,&args,(char *)0,help);

ierr = PetscOptionsGetInt(PETSC_NULL,"-",&m,PETSC_NULL);CHKERRQ(ierr);

ierr = PetscOptionsGetInt(PETSC_NULL,"-",&n,PETSC_NULL);CHKERRQ(ierr);

265

/* Организуем матрицу и правосторонний вектор, который определяет линейную

систему, Ax = b. Создаем параллельную матрицу, описывая только ее глобальные

размеры. Когда используется MatCreate(), формат матрицы может быть описан на

этапе исполнения. */

ierr=MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,

m*n,m*n,&A); CHKERRQ(ierr);

ierr = MatSetFromOptions(A);CHKERRQ(ierr);

/* Матрица разделена крупными кусками строк по процессам. Определяем, какие

строки матрицы размещены локально.*/

ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);

/* Размещаем матричные элементы. Каждому процессу нужно разместить только

те элементы, которые принадлежат ему локально (все нелокальные элементы

будут посланы в соответствующий процессор во время сборки матрицы). */

for (I=Istart; I<Iend; I++)

{ v = -1.0; i = I/n; j = I – i*n;

if (i>0)

{ J = I – n;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);

CHKERRQ(ierr);

}

if (i<m-1)

{ J = I + n;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);

CHKERRQ(ierr);

}

if (j>0)

{ J = I – 1;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);

CHKERRQ(ierr);

}

if (j<n-1)

{ J = I + 1;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);

CHKERRQ(ierr);

}

v = 4.0;

ierr = MatSetValues(A,1,&I,1,&I,&v,INSERT_VALUES);

CHKERRQ(ierr);

}

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

/* Создаем параллельные векторы. Формируем вектор и затем дублируем, если

необходимо. Когда в этом примере используется VecCreate() и

266

VecSetFromOptions(), указывается только глобальный размер вектора; парал-

лельное разбиение определяется на этапе исполнения. Когда решается линейная

система, векторы и матрицы обязаны быть разбиты соответственно. PETSc авто-

матически генерирует соответствующее разбиение матриц и векторов, когда

MatCreate() и VecCreate() используются с тем же самым коммуникатором. Поль-

зователь может альтернативно описать размеры локальных векторов и матриц,

когда необходимо более сложное разбиение (путем замещения аргумента

PETSC_DECIDE в VecCreate()). */

ierr = VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,&u);

CHKERRQ(ierr);

ierr = VecSetFromOptions(u);CHKERRQ(ierr);

ierr = VecDuplicate(u,&b);CHKERRQ(ierr);

ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

/* Установим точное решение, затем вычислим правосторонний вектор. По умол-

чанию используем точное решение вектора для случая, когда все элементы

вектора равны единице. Альтернативно, используя опцию – random_sol, фор-

мируем решение вектора со случайными компонентами. */

ierr = PetscOptionsHasName(PETSC_NULL,"-random_exact_sol",&flg);

CHKERRQ(ierr);

if (flg)

{ ierr = PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);

CHKERRQ(ierr);

ierr = VecSetRandom(rctx,u);CHKERRQ(ierr);

ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);

}

else ierr = VecSet(&one,u);CHKERRQ(ierr);

ierr = MatMult(A,u,b);CHKERRQ(ierr);

/* Выводим точное решение вектора, если необходимо */

ierr = PetscOptionsHasName(PETSC_NULL,"-view_exact_sol",&flg);

CHKERRQ(ierr);

if (flg)

{ ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);

CHKERRQ(ierr);

}

/* Создаем метод решения системы линейных уравнений */

ierr = SLESCreate(PETSC_COMM_WORLD,&sles);CHKERRQ(ierr);

/* Устанавливаем операторы. Здесь матрица, которая определяет линейную

систему, служит как preconditioning матрица. */

ierr = SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);

CHKERRQ(ierr);

/* Устанавливается метод решения */

ierr = SLESGetKSP(sles,&ksp);CHKERRQ(ierr);

ierr = KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,

267

PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

ierr = SLESSetFromOptions(sles);CHKERRQ(ierr);

/* Решается линейная система */

ierr = SLESSolve(sles,b,x,&its);CHKERRQ(ierr);

/* Проверяется решение norm *= sqrt(1.0/((m+1)*(n+1))); */

ierr = VecAXPY(&neg_one,u,x);CHKERRQ(ierr);

ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);

/* Печатается информация о сходимости. PetscPrintf() создает одно предложение

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

ляется PetscFPrintf(), которая печатает в файл. */

ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations d\n",

norm,its); CHKERRQ(ierr);

/* Освобождается рабочее пространство. Все PETSc объекты должны быть разру-

шены, когда они больше не нужны. */

ierr = SLESDestroy(sles);CHKERRQ(ierr);ierr = VecDestroy(u);CHKERRQ(ierr);

ierr = VecDestroy(x);CHKERRQ(ierr);ierr = VecDestroy(b);CHKERRQ(ierr);

ierr = MatDestroy(A);CHKERRQ(ierr);

ierr = PetscFinalize();CHKERRQ(ierr);

return 0;

}

КОНТРОЛЬНЫЕ ВОПРОСЫ К ГЛАВЕ 13

Контрольные вопросы к 13.1

1. Для чего предназначена библиотека ScaLAPACK?

2. Для чего предназначена библиотека PETSc?

3. Опишите уровень операций в библиотеки BLAS?

4. Для чего предназначена библиотека BLAS?

5. Для чего предназначена библиотека LAPACK?

6. Каково назначение библиотеки BLAСS?

7. Объясните назначение параметров функции vTRSD2D.

8. Какие четыре шага надо сделать, чтобы выполнить программу ScaLAPACK?

9. Опишите каждый из этих четырех шагов.

10. Объясните на примере программы умножения матриц особенности работы

библиотеки ScaLAPACK.

Контрольные вопросы к 13.2

1. Какие задачи позволяет решать библиотека PETSc?

2. Нужно ли вызывать MPI_Init при написании PETSc программы?

3. В каком случае необходимо создание коммуникатора и как это сделать?

4. Какие типы векторов существуют в библиотеке? Как создать вектор того или

иного типа?

5. В каком случае нужно создавать параллельные векторы?

6. Можно ли распределить компоненты параллельного вектора по процессам не

равными частями?

268

7.

Как определить номера строк в локальном процессе при создании параллель-

ной матрицы?

8. Какие матричные операции библиотеки можно использовать для параллель-

ного решения задачи умножения матрицы на матрицу?

9. Необходим ли MPICH для инсталляции PETSc?

10. Как выполнить PETSc-программу на 4 процессах, задавая различные размер-

ности задачи?

269

ПРИЛОЖЕНИЯ

Приложение 1. КОНСТАНТЫ ДЛЯ ЯЗЫКОВ С И FORTRAN

Здесь приведены константы, определенные в файлах mpi.h (для C) и mpif.h

(для языка Fortran).

/* возвращаемые коды (как для C, так и для Fortran) */

MPI_SUCCESS

MPI_ERR_BUFFER

MPI_ERR_COUNT

MPI_ERR_TYPE

MPI_ERR_TAG

MPI_ERR_COMM

MPI_ERR_RANK

MPI_ERR_REQUEST

MPI_ERR_ROOT

MPI_ERR_GROUP

MPI_ERR_OP

MPI_ERR_TOPOLOGY

MPI_ERR_DIMS

MPI_ERR_ARG

MPI_ERR_UNKNOWN

MPI_ERR_TRUNCATE

MPI_ERR_OTHER

MPI_ERR_INTERN

MPI_ERR_PENDING

MPI_ERR_IN_STATUS

MPI_ERR_LASTCODE

/* константы (как для C, так и для Fortran) */

MPI_BOTTOM

MPI_PROC_NULL

MPI_ANY_SOURCE

MPI_ANY_TAG

MPI_UNDEFINED

MPI_BSEND_OVERHEAD

MPI_KEYVAL_INVALID

/* размер статуса и резервируемые индексные значения (Fortran) */

MPI_STATUS_SIZE

MPI_SOURCE

MPI_TAG

MPI_ERROR

/* описатели обработчика ошибок (C и Fortran) */

MPI_ERRORS_ARE_FATAL

MPI_ERRORS_RETURN

/* максимальные размеры строк */

270

MPI_MAX_PROCESSOR_NAME

MPI_MAX_ERROR_STRING

/* элементарные типы данных (C) */

MPI_CHAR

MPI_SHORT

MPI_INT

MPI_LONG

MPI_UNSIGNED_CHAR

MPI_UNSIGNED_SHORT

MPI_UNSIGNED

MPI_UNSIGNED_LONG

MPI_FLOAT

MPI_DOUBLE

MPI_LONG_DOUBLE

MPI_BYTE

MPI_PACKED

/* элементарные типы данных(Fortran) */

MPI_INTEGER

MPI_REAL

MPI_DOUBLE_PRECISION

MPI_COMPLEX

MPI_LOGICAL

MPI_CHARACTER

MPI_BYTE

MPI_PACKED

/* типы данных для функций редукции (C) */

MPI_FLOAT_INT

MPI_DOUBLE_INT

MPI_LONG_INT

MPI_2INT

MPI_SHORT_INT

MPI_LONG_DOUBLE_INT

/* типы данных для функций редукции (Fortran) */

MPI_2REAL

MPI_2DOUBLE_PRECISION

MPI_2INTEGER

/* дополнительные типы данных (Fortran) */

MPI_INTEGER1

MPI_INTEGER2

MPI_INTEGER4

MPI_REAL2

MPI_REAL4

MPI_REAL8

etc.

/* дополнительные типы данных (C) */

MPI_LONG_LONG_INT

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

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