Чтобы помочь пользователю начать освоение PETSc немедленно, следует
начать с простого однопроцессорного примера,
который решает одномерную задачу Лапласа с конечными разностями.
Этот последовательный код, который находится в файле
${PETSC_DIR}/src/sles/examples/tutorials/ex1.c,
иллюстрирует решение линейной системы с помощью библиотек SLES, интерфейсов
предобработчиков, методов подпространств Крылова и прямых линейных
решателей PETSc. По мере движения по коду выделяется несколько наиболее
важных частей примера:
/* Program usage: mpirun ex1 [-help] [all PETSc options] */
static char help[] = "Решает тридиагональную линейную систему
через SLES.\n\n";
/*T
Концепция: SLES^решение системы линейных уравнений
Процессоры: 1
T*/
/*
Включаем "petscsles.h" чтобы использовать решатели SLES .
Отметьте, что этот файл автоматически включает:
petsc.h - основные процедуры PETSc petscvec.h - векторы
petscsys.h - системные процедуры petscmat.h - матрицы
petscis.h - индексные множества petscksp.h - методы
подпространств Крылова
petscviewer.h - просмотрщики petscpc.h - предобработчики
Замечание:
Соответствующий параллельный пример находится в файле ex23.c
*/
#include "petscsles.h"
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
Vec x, b, u; /* приближенное решение, RHS, точное решение */
Mat A; /* матрица линейной системы */
SLES sles; /* контекст линейного решателя */
PC pc; /* контекст предобработчика */
KSP ksp; /* контекст метода подпространств Крылова */
PetscReal norm; /* допустимая ошибка решения */
int ierr,i,n = 10,col[3],its,size;
PetscScalar neg_one = -1.0,one = 1.0,value[3];
PetscInitialize(&argc,&args,(char *)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);
CHKERRQ(ierr);
if (size != 1) SETERRQ(1,"Это однопроцессорный пример!");
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
CHKERRQ(ierr);
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Вычисляет матрицу и правосторонний вектор, который определяет
линейную систему Ax = b.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
/*
Создает векторы. Заметьте, что первый вектор формируется с нуля,
а затем размножается, если нужно.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) x, "Решение");CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
/*
Создает матрицу. При использовании MatCreate() формат матрицы
можно определять в процессе выполнения.
Замечание по производительности:
Для задач существенного объема предварительное распределение
памяти для матрицы критично для достижения хорошей
производительности. Поскольку предварительное распределение
невозможно средствами общей процедуры создания матрицы MatCreate(),
рекомендуется для конкретных задач использовать
специализированные процедуры, например:
MatCreateSeqAIJ() - последовательная AIJ
(упакованная разреженная строка),
MatCreateSeqBAIJ() - блочная AIJ
См. главу о матрицах руководства пользователя
для более полной информации.
*/
ierr = MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&A);
CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
/*
Собирает матрицу
*/
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=1; i<n-1; i++)
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
CHKERRQ(ierr);
i = n - 1; col[0] = n - 2; col[1] = n - 1;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
CHKERRQ(ierr);
i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*
Установим точное решение; затем вычислим правосторонний вектор.
*/
ierr = VecSet(&one,u);CHKERRQ(ierr);
ierr = MatMult(A,u,b);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 = SLESGetPC(sles,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-7,
PETSC_DEFAULT,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);
/*
Смотреть информацию от решателя; вместо этого можно
использовать опцию
-sles_view для вывода этой информации на экран после
завершения SLESSolve().
*/
ierr = SLESView(sles,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Проверяем решение и очищаем его
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
/*
Проверяем ошибки
*/
ierr = VecAXPY(&neg_one,u,x);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,
"Погрешность %A, Итераций %d\n",norm,its); CHKERRQ(ierr);
/*
Освобождаем рабочее пространство. Все объекты PETSc
должны быть уничтожены, как только они станут ненужными.
*/
ierr = VecDestroy(x);CHKERRQ(ierr);
ierr = VecDestroy(u);CHKERRQ(ierr);
ierr = VecDestroy(b);CHKERRQ(ierr);
ierr = MatDestroy(A);CHKERRQ(ierr);
ierr = SLESDestroy(sles);CHKERRQ(ierr);
/* Всегда вызываем PetscFinalize() перед выходом из программы.
Эта процедура
- финализирует библиотеки PETSc и MPI
- предоставляет общую и диагностическую информацию,
если указаны определенные опции времени выполнения
(например, -log_summary).
*/
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
}