Поскольку PETSc использует модель передачи сообщений для параллельного программирования и пользуется MPI для межпроцессорного взаимодействия, пользователь может свободно употреблять в своем коде процедуры MPI там, где это необходимо. Однако по умолчанию пользователь огражден от многих деталей обмена сообщениями внутри PETSc, поскольку они скрыты в таких параллельных объектах, как векторы, матрицы и решатели. К тому же, PETSc продоставляет такие инструменты, как обобщенную сборку/рассылку векторов и распределенные массивы, чтобы облегчить управление параллельными данными. Помните, что пользователь должен определить коммуникатор перед созданием любого объекта PETSc (такого как вектор, матрица или решатель), чтобы указать процессоры, на которые будет распределен объект. Например, как упоминалось выше, командами для создания матрицы, вектора и линейного решателя являются:
VecCreate (MPI_Comm comm ,Vec *x);
SLESCreate (MPI_Comm comm ,SLES *sles);
/* Запуск программы: mpirun -np <procs> ex2 [-help]
[all PETSc options] */
static char help[] =
"Параллельное решение линейной системы через SLES.\n\
Входные параметры включают:\n\
-random_exact_sol : использовать случайный вектор\
точного решения\n\
-view_exact_sol : вывод вектора точного решения в stdout\n\
-m <mesh_x> : количество узлов сетки по x\n\
-n <mesh_n> : количество узлов сетки по y\n\n";
/*T
Идея: параллельный пример, основанный на SLES;
Идея: SLES^Laplacian, 2d
Идея: Laplacian, 2d
Процессоры: n
T*/
/*
Включаем "petscsles.h", чтобы использовать решатели SLES.
Этот файл автоматически включает:
petsc.h - основные процедуры PETSc petscvec.h - векторы
petscsys.h - системные процедуры petscmat.h - матрицы
petscis.h - индексные множества petscksp.h - методы
подпространств Крылова
petscviewer.h - просмотрщики petscpc.h - предобработчики
*/
#include "petscsles.h"
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args) {
Vec x,b,u; /* приближенное решение, RHS,
точное решение */
Mat A; /* матрица линейной системы */
SLES sles; /* контекст линейного решателя */
PetscRandom rctx; /* контекст генератора случайных чисел */
PetscReal norm; /* погрешность решения */
int i,j,I,J,Istart,Iend,ierr,m = 8,n = 7,its;
PetscTruth flg;
PetscScalar v,one = 1.0,neg_one = -1.0;
KSP ksp;
PetscInitialize(&argc, &args, (char*)0, help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,PETSC_NULL);
CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
CHKERRQ(ierr);
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Вычисляем матрицу и правосторонний вектор, определяющий
линейную систему Ax = b.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
/*
Создаем параллельную матрицу, определяя только ее общие размеры.
При использовании MatCreate() формат матрицы может быть задан
во время выполнения.
Аналогично параллельное разделение матрицы определяется PETSc
во время выполнения.
Для проблем существенной размерности, предварительное
распределение памяти для матрицы влияет на
производительность.
Поскольку предварительное распределение невозможно при
использовании общей процедуры создания матрицы MatCreate(),
рекомендуется для практических задач вместо этого использовать
процедуры создания матриц специальных форматов, например:
MatCreateMPIAIJ() - параллельная AIJ
(упакованная разреженная строка)
MatCreateMPIBAIJ() - параллельная блочная AIJ
См. главу о матрицах в руководстве пользователя
для детальной информации.
*/
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,
m*n,&A);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
/*
В настоящее время все параллельные форматы матриц PETSc
разделяются среди процессоров на непрерывные последовательности
строк. Определим, какие строки матрицы получены локально
*/
ierr = MatGetOwnershipRange(A,\&Istart,\&Iend);CHKERRQ(ierr);
/*
Параллельно устанавливаем элементы матрицы для 2-D
пятиточечного шаблона:
- Каждый процессор должен вставить только элементы, которыми
он владеет локально (любые нелокальные элементы будут
пересланы соответствующему процессору во время
сборки матрицы).
- Всегда указывайте глобальные номера строк и столбцов
элементов матрицы.
Здесь используется не совсем стандартная нумерация,
которая перебирает сначала все неизвестные для x = h, затем
для x = 2h и т.д.; поэтому Вы видите J = I +- n вместо
J = I +- m, как Вы ожидали. Более стандартная нумерация сначала
проходит все переменные для y = h, затем y = 2h и т.д.
*/
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);
}
/*
Собираем матрицу, используя двухшаговый процесс:
MatAssemblyBegin(), MatAssemblyEnd()}
Можно выполнять вычисления, пока сообщения передаются,
если поместить код между двумя этими операторами.
*/
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*
Создаем параллельные векторы.
- Первый вектор формируется с нуля, а затем делаются дубликаты,
если нужно.
- При использовании VecCreate(), VecSetSizes и
VecSetFromOptions() в этом примере, определяются только
глобальные размеры вектора; параллельное разделение
определяется во время выполнения.
- При решении линейной системы и векторы, и матрицы должны
соответственно разделяться. PETSc автоматически генерирует
правильно разделенные матрицы и векторы, если процедуры
MatCreate() и VecCreate() используются с одним и тем же
коммуникатором.
- Пользователь может самостоятельно определить локальные
размеры векторов и матриц, если требуется более сложное
разделение (заменив аргумент PETSC_DECIDE в процедуре
VecSetSizes() ниже).
*/
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
/*
Установим точное решение; затем вычислим правосторонний
вектор. По умолчанию мы используем точное решение вектора
со всеми элементами 1.0; однако использование
опции времени выполнения -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);
/*
Установим операторы. Здесь матрица, определяющая линейную
систему, служит также матрицей предобработчиков.
*/
ierr = SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);
CHKERRQ(ierr);
/*
Установим значения линейного решателя по умолчанию для
этой задачи (необязательно):
- Извлекая контексты KSP и PC из контекста SLES, затем можно
непосредственно вызывать любые процедуры KSP и PC
для установки различных опций.
- Следующие два оператора опциональны; все параметры
могут быть альтернативно заданы во время выполнения
через SLESSetFromOptions().
Все значения по умолчанию можно переопределить
во время выполнения, как показано ниже.
*/
ierr = SLESGetKSP(sles,&ksp);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,
PETSC_DEFAULT,PETSC_DEFAULT); CHKERRQ(ierr);
/*
Установим опции времени выполнения, например,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
Эти опции могут переопределить значения, указаннаые выше
до тех пор, пока SLESSetFromOptions() будет вызван после любой
процедуры настройки.
*/
ierr = SLESSetFromOptions(sles);CHKERRQ(ierr);
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Решаем линейную систему
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
ierr = SLESSolve(sles,b,x,&its);CHKERRQ(ierr);
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Проверяем решение и очищаем его
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
/*
Проверяем ошибку
*/
ierr = VecAXPY(&neg_one,u,x);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
/*
Масштабируем ошибку
*/
/* norm *= sqrt(1.0/((m+1)*(n+1))); */
/*
Выводим информацию о сходимости. 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(b);CHKERRQ(ierr);
ierr = VecDestroy(x);CHKERRQ(ierr);
ierr = MatDestroy(A);CHKERRQ(ierr);
/*
Всегда вызывайте PetscFinalize() перед выходом из программы.
Эта процедура
- финализирует библиотеки PETSc, а также MPI
- выдает итоговую и диагностическую информацию, если указаны
определенные опции времени выполнения
(например, -log_summary).
*/
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
}