6.3 KSP: решатели СЛАУ

Объект KSP является сердцем PETSc, поскольку он обеспечивает унифицированный и эффективный доступ ко всем средствам для решения систем линейных алгебраических уравнений (СЛАУ), включая параллельные и последовательные, прямые и итерационные методы. KSP предназначен для решения систем вида:

Ax  = b.
(6.1)

где A – матричное представление линейного оператора, b – вектор правых частей, x – вектор решения.

Для того чтобы создать решатель линейной системы создаем объект KSP:

    KSPCreate(MPI_Comm comm, KSP *ksp);

Здесь comm есть MPI коммуникатор, ksp – вновь сформированный контекст решателя.

Перед решением линейной системы с помощью KSP пользователь обязан вызвать следующую процедуру, чтобы задать матрицы, связанные с линейной системой:

    KSPSetOperators(KSP ksp, Mat Amat, Mat Pmat, MatStructure flag);

где Amat – матрица линейной системы, Pmat – матрица переобусловливания (то есть матрица, из которой должен быть сконструирован переобусловливатель). Часто в качестве матрицы предобуславливания задается матрица Amat, которая определяет линейную систему. Аргумент flag может быть использован, чтобы избежать ненужной работы при многократном решении линейных систем того же размера и тем же методом переобусловливания. Когда решается только одна система, этот флаг игнорируется. Пользователь может установить flag следующим образом:

Если имеется сомнение относительно структуры матрицы, нужно использовать флаг DIFFERENT_NONZERO_PATTERN.

И наконец, чтобы решить СЛАУ вызываем метод

    KSPSolve(KSP ksp, Vec b, Vec x);

где b и x соответственно обозначают правосторонний вектор и вектор решения.

После расчета можем получить количество проделанных итераций, вызвав функцию

    KSPGetIterationNumber(KSP ksp, int *its);

Параметр its содержит либо номер итерации, на которой сходимость была успешно достигнута, или номер итерации, на которой была обнаружена расходимость или прерывание.

Когда контекст больше не нужен, он удаляется командой

    KSPDestroy(KSP ksp);

По умолчанию KSP использует метод GMRES с предобуславливателем ILU для случая последовательного решения и block Jacobi – в случае параллельного запуска. Также доступны и многие другие методы решения. Для того чтобы вручную задать метод решения вызываем функцию

    KSPSetType(KSP ksp,KSPType method);

где method – метод, может быть одним из: KSPRICHARDSON, KSPCHEBYCHEV, KSPCG, KSPGMRES, KSPTCQMR, KSPBCGS, KSPCGS, KSPTFQMR, KSPCR, KSPLSQR, KSPBICG или KSPPREONLY (Таблица 6.3). Метод KSP также можно установить опцией командной строки -ksp_type, задав одно из следующих значений: richardson, chebychev, cg, gmres, tcqmr, bcgs, cgs, tfqmr, cr, lsqr, bicg или preonly.

Для методов Richardson, Chebychev и GMRES, можно задать некоторые их специфические свойства

    KSPRichardsonSetScale(KSP ksp, double damping_factor); 
    KSPChebychevSetEigenvalues(KSP ksp, double emax, double emin); 
    KSPGMRESSetRestart(KSP ksp, int max_steps);

Значениями параметров по умолчанию являются: damping_factor = 1.0, emax = 0.01, emin = 100.0, и max_steps = 30. Эти параметры можно задать с использованием следующих опций: -ksp_gmres_restart <n> и -ksp_richardson_scale <factor>.

По умолчанию для ортогонализации матрицы Hessenberg’a в методе GMRES используется ортогонализация Грамма-Шмидта. Метод ортогонализации можно поменять, задав опцию командной строки -ksp_gmres_irorthog или вызвать функцию

    KSPGMRESSetOrthogonalization(KSP ksp, KSPGMRESModifiedGramSchmidtOrthogonalization);

которая задает модифицированную ортоганализацию Грамма-Шмидта. Также можно задать немодифицированный (классический ) метод Грамма-Шмидта

    KSPGMRESSetOrthogonalization(KSP ksp, KSPGMRESUnmodifiedGramSchmidtOrthogonalization);

или опцией -ksp_gmres_unmodifiedgramschmidt. Заметим, что этот алгоритм численно неустойчив, но имеет несколько лучшие скоростные характеристики.

В KSP, по умолчанию, в качестве начального значения берет вектор, получаемый обнулением начального значения. Чтобы использовать ненулевое начальное приближение, нужно вызвать

    KSPSetInitialGuessNonzero(KSP ksp);

Чтобы решить линейную систему прямым методом (в настоящее время поддерживается только для последовательных матриц), можно задать следующие опции -pc_type lu -ksp_type preonly.





Метод KSPType Название опции



Richardson KSPRICHARDSON richardson



Chebychev KSPCHEBYCHEV chebychev



Conjugate Gradient KSPCG cg



BiConjugate Gradient KSPBICG bicg



Generalized Minimal Residual KSPGMRES gmres



BiCGSTAB KSPBCGS bcgs



Conjugate Gradient Squared KSPCGS cgs



Transpose-Free Quasi-Minimal Residual (1) KSPTFQMR tfqmr



Transpose-Free Quasi-Minimal Residual (2) KSPTCQMR tcqmr



Conjugate Residual KSPCR cr



Least Squares Method KSPLSQR lsqr



Shell for no KSP method KSPPREONLY preonly




Таблица 6.3: Методы решения СЛАУ в KSP

Порядок сходимости методов подпространств Крылова для линейных систем сильно зависит от их спектра. Для изменения спектра и, отсюда, для ускорения сходимости итерационных методов как правило используется переобусловливание. Переобусловливание может быть применено к системе (6.1) путем преобразования

   -1    -1           -1
(M L AM R  )(MRx ) = ML  b,
(6.2)

гдеML и MR – предобуславливатели. Если ML = I, то получим предобуславливание справа:

                 -1
r ≡ b - Ax = b- M R MRx.
(6.3)

При MR = I – предобуславливание слева:

r  ≡ M -1b- M - 1Ax = M - 1r.
 L     L      L        L
(6.4)

По умолчанию все реализации KSP используют левое переобусловливание. Правое переобусловливание может быть задано опцией -ksp_right_pc или вызовом процедуры

    KSPSetPreconditionerSide(KSP ksp, PCSide PC_RIGHT);

Контроль ошибки. Для определения сходимости используется l2 норма невязки, для ее оценки применяются следующие три величины: rtol, atol и dtol. На k-ой итерации, для оценки сходимости, используется следующее неравенство

||rk||2 < max (rtol* ||b||2,atol),
(6.5)

где rk = b - Axk . Процесс считается расходящимся при

||rk|| > dtol*||b||
    2          2
(6.6)

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

    KSPSetTolerances(KSP ksp, double rtol, double atol, double dtol, int maxits);

Здесь при задании параметра как PETSC_DEFAULT берется значение по умолчанию: rtol = 10-5, atol = 10-50, dtol = 105 и maxits = 105. Эти параметры можно задать и опциями командной строки -ksp_rtol <rtol>, -ksp_atol <atol>, -ksp_divtol <dtol> и -ksp_max_it <its>.

По умолчанию KSP решатели не выводят на экран информацию о проделанных итерациях. Для вывода информации на экран можно указать опцию командной строки -ksp_monitor. Также можно вызвать функцию

    KSPSetMonitor(KSP ksp, int (*mon)(KSP ksp, int it, double rnorm, void *ctx), void *ctx, int (*mondestroy)(void *));

В функции mon аргументами являются количество итераций, it, и l2 норма невязки, rnorm.

Предобуславливатели. Методы подпространств Крылова обычно используется с переобусловливателем. Опцией командной строки можно задать нужный тип предобуславливателя, -pc_type <methodname> или

    PCSetType(PC pc, PCType method);

В таблице 6.4 представлены основные методы переобусловливания, поддерживаемые в PETSc. Каждый предобусловливатель имеет ряд связанных с ним опций, которые могут быть установлены процедурами и опциями командной строки. Все имена таких процедур и команд имеют форму PC<TYPE>Option и -pc_<type>_option [value].





Метод PCType Название опции



Jacobi PCJACOBI jacobi



Block Jacobi PCBJACOBI bjacobi



SOR (and SSOR) PCSOR sor



SOR with Eisenstat trick PCEISENSTAT eisenstat



Incomplete Cholesky PCICC icc



Incomplete LU PCILU ilu



Additive Schwarz PCASM asm



Linear solver PCKSP ksp



Combination of preconditionersPCCOMPOSITE composite



LU PCLU lu



Cholesky PCCholesky cholesky



No preconditioning PCNONE none



Shell for user-defined PC PCSHELL shell




Таблица 6.4: Предобуславливатели в PETSc