SNES: Nonlinear Solvers#
Note
This chapter is being cleaned up by Jed Brown. Contributions are welcome.
The solution of largescale nonlinear problems pervades many facets of
computational science and demands robust and flexible solution
strategies. The SNES
library of PETSc provides a powerful suite of
datastructureneutral numerical routines for such problems. Built on
top of the linear solvers and data structures discussed in preceding
chapters, SNES
enables the user to easily customize the nonlinear
solvers according to the application at hand. Also, the SNES
interface is identical for the uniprocess and parallel cases; the only
difference in the parallel version is that each process typically forms
only its local contribution to various matrices and vectors.
The SNES
class includes methods for solving systems of nonlinear
equations of the form
where \(\mathbf{F}: \, \Re^n \to \Re^n\). Newtonlike methods provide the
core of the package, including both line search and trust region
techniques. A suite of nonlinear Krylov methods and methods based upon
problem decomposition are also included. The solvers are discussed
further in The Nonlinear Solvers. Following the PETSc design
philosophy, the interfaces to the various solvers are all virtually
identical. In addition, the SNES
software is completely flexible, so
that the user can at runtime change any facet of the solution process.
PETSc’s default method for solving the nonlinear equation is Newton’s method. The general form of the \(n\)dimensional Newton’s method for solving (3) is
where \(\mathbf{x}_0\) is an initial approximation to the solution and \(\mathbf{J}(\mathbf{x}_k) = \mathbf{F}'(\mathbf{x}_k)\), the Jacobian, is nonsingular at each iteration. In practice, the Newton iteration (4) is implemented by the following two steps:
Other defectcorrection algorithms can be implemented by using different choices for \(J(\mathbf{x}_k)\).
Basic SNES Usage#
In the simplest usage of the nonlinear solvers, the user must merely provide a C, C++, or Fortran routine to evaluate the nonlinear function (3). The corresponding Jacobian matrix can be approximated with finite differences. For codes that are typically more efficient and accurate, the user can provide a routine to compute the Jacobian; details regarding these applicationprovided routines are discussed below. To provide an overview of the use of the nonlinear solvers, browse the concrete example in ex1.c or skip ahead to the discussion.
Listing: src/snes/tutorials/ex1.c
static char help[] = "Newton's method for a twovariable system, sequential.\n\n";
/*
Include "petscsnes.h" so that we can use SNES solvers. Note that this
file automatically includes:
petscsys.h  base PETSc routines petscvec.h  vectors
petscmat.h  matrices
petscis.h  index sets petscksp.h  Krylov subspace methods
petscviewer.h  viewers petscpc.h  preconditioners
petscksp.h  linear solvers
*/
/*F
This examples solves either
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{x^2_0 + x_0 x_1  3}{x_0 x_1 + x^2_1  6}
\end{equation}
or if the {\tt hard} options is given
\begin{equation}
F\genfrac{(}{)}{0pt}{}{x_0}{x_1} = \genfrac{(}{)}{0pt}{}{\sin(3 x_0) + x_0}{x_1}
\end{equation}
F*/
#include <petscsnes.h>
/*
Userdefined routines
*/
extern PetscErrorCode FormJacobian1(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode FormFunction1(SNES,Vec,Vec,void*);
extern PetscErrorCode FormJacobian2(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode FormFunction2(SNES,Vec,Vec,void*);
int main(int argc,char **argv)
{
SNES snes; /* nonlinear solver context */
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
Vec x,r; /* solution, residual vectors */
Mat J; /* Jacobian matrix */
PetscMPIInt size;
PetscScalar pfive = .5,*xx;
PetscBool flg;
PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Example is only for sequential runs");
/*                                  
Create nonlinear solver context
                                  */
PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
/*                                   
Create matrix and vector data structures; set corresponding routines
                                   */
/*
Create vectors for solution and nonlinear function
*/
PetscCall(VecCreate(PETSC_COMM_WORLD,&x));
PetscCall(VecSetSizes(x,PETSC_DECIDE,2));
PetscCall(VecSetFromOptions(x));
PetscCall(VecDuplicate(x,&r));
/*
Create Jacobian matrix data structure
*/
PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,2,2));
PetscCall(MatSetFromOptions(J));
PetscCall(MatSetUp(J));
PetscCall(PetscOptionsHasName(NULL,NULL,"hard",&flg));
if (!flg) {
/*
Set function evaluation routine and vector.
*/
PetscCall(SNESSetFunction(snes,r,FormFunction1,NULL));
/*
Set Jacobian matrix data structure and Jacobian evaluation routine
*/
PetscCall(SNESSetJacobian(snes,J,J,FormJacobian1,NULL));
} else {
PetscCall(SNESSetFunction(snes,r,FormFunction2,NULL));
PetscCall(SNESSetJacobian(snes,J,J,FormJacobian2,NULL));
}
/*                                  
Customize nonlinear solver; set runtime options
                                   */
/*
Set linear solver defaults for this problem. By extracting the
KSP and PC contexts from the SNES context, we can then
directly call any KSP and PC routines to set various options.
*/
PetscCall(SNESGetKSP(snes,&ksp));
PetscCall(KSPGetPC(ksp,&pc));
PetscCall(PCSetType(pc,PCNONE));
PetscCall(KSPSetTolerances(ksp,1.e4,PETSC_DEFAULT,PETSC_DEFAULT,20));
/*
Set SNES/KSP/KSP/PC runtime options, e.g.,
snes_view snes_monitor ksp_type <ksp> pc_type <pc>
These options will override those specified above as long as
SNESSetFromOptions() is called _after_ any other customization
routines.
*/
PetscCall(SNESSetFromOptions(snes));
/*                                  
Evaluate initial guess; then solve nonlinear system
                                   */
if (!flg) {
PetscCall(VecSet(x,pfive));
} else {
PetscCall(VecGetArray(x,&xx));
xx[0] = 2.0; xx[1] = 3.0;
PetscCall(VecRestoreArray(x,&xx));
}
/*
Note: The user should initialize the vector, x, with the initial guess
for the nonlinear solver prior to calling SNESSolve(). In particular,
to employ an initial guess of zero, the user should explicitly set
this vector to zero by calling VecSet().
*/
PetscCall(SNESSolve(snes,NULL,x));
if (flg) {
Vec f;
PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD));
PetscCall(SNESGetFunction(snes,&f,0,0));
PetscCall(VecView(r,PETSC_VIEWER_STDOUT_WORLD));
}
/*                                  
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
                                   */
PetscCall(VecDestroy(&x)); PetscCall(VecDestroy(&r));
PetscCall(MatDestroy(&J)); PetscCall(SNESDestroy(&snes));
PetscCall(PetscFinalize());
return 0;
}
/*  */
/*
FormFunction1  Evaluates nonlinear function, F(x).
Input Parameters:
. snes  the SNES context
. x  input vector
. ctx  optional userdefined context
Output Parameter:
. f  function vector
*/
PetscErrorCode FormFunction1(SNES snes,Vec x,Vec f,void *ctx)
{
const PetscScalar *xx;
PetscScalar *ff;
/*
Get pointers to vector data.
 For default PETSc vectors, VecGetArray() returns a pointer to
the data array. Otherwise, the routine is implementation dependent.
 You MUST call VecRestoreArray() when you no longer need access to
the array.
*/
PetscCall(VecGetArrayRead(x,&xx));
PetscCall(VecGetArray(f,&ff));
/* Compute function */
ff[0] = xx[0]*xx[0] + xx[0]*xx[1]  3.0;
ff[1] = xx[0]*xx[1] + xx[1]*xx[1]  6.0;
/* Restore vectors */
PetscCall(VecRestoreArrayRead(x,&xx));
PetscCall(VecRestoreArray(f,&ff));
return 0;
}
/*  */
/*
FormJacobian1  Evaluates Jacobian matrix.
Input Parameters:
. snes  the SNES context
. x  input vector
. dummy  optional userdefined context (not used here)
Output Parameters:
. jac  Jacobian matrix
. B  optionally different preconditioning matrix
. flag  flag indicating matrix structure
*/
PetscErrorCode FormJacobian1(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
{
const PetscScalar *xx;
PetscScalar A[4];
PetscInt idx[2] = {0,1};
/*
Get pointer to vector data
*/
PetscCall(VecGetArrayRead(x,&xx));
/*
Compute Jacobian entries and insert into matrix.
 Since this is such a small problem, we set all entries for
the matrix at once.
*/
A[0] = 2.0*xx[0] + xx[1]; A[1] = xx[0];
A[2] = xx[1]; A[3] = xx[0] + 2.0*xx[1];
PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES));
/*
Restore vector
*/
PetscCall(VecRestoreArrayRead(x,&xx));
/*
Assemble matrix
*/
PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
if (jac != B) {
PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
}
return 0;
}
/*  */
PetscErrorCode FormFunction2(SNES snes,Vec x,Vec f,void *dummy)
{
const PetscScalar *xx;
PetscScalar *ff;
/*
Get pointers to vector data.
 For default PETSc vectors, VecGetArray() returns a pointer to
the data array. Otherwise, the routine is implementation dependent.
 You MUST call VecRestoreArray() when you no longer need access to
the array.
*/
PetscCall(VecGetArrayRead(x,&xx));
PetscCall(VecGetArray(f,&ff));
/*
Compute function
*/
ff[0] = PetscSinScalar(3.0*xx[0]) + xx[0];
ff[1] = xx[1];
/*
Restore vectors
*/
PetscCall(VecRestoreArrayRead(x,&xx));
PetscCall(VecRestoreArray(f,&ff));
return 0;
}
/*  */
PetscErrorCode FormJacobian2(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
{
const PetscScalar *xx;
PetscScalar A[4];
PetscInt idx[2] = {0,1};
/*
Get pointer to vector data
*/
PetscCall(VecGetArrayRead(x,&xx));
/*
Compute Jacobian entries and insert into matrix.
 Since this is such a small problem, we set all entries for
the matrix at once.
*/
A[0] = 3.0*PetscCosScalar(3.0*xx[0]) + 1.0; A[1] = 0.0;
A[2] = 0.0; A[3] = 1.0;
PetscCall(MatSetValues(B,2,idx,2,idx,A,INSERT_VALUES));
/*
Restore vector
*/
PetscCall(VecRestoreArrayRead(x,&xx));
/*
Assemble matrix
*/
PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
if (jac != B) {
PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
}
return 0;
}
To create a SNES
solver, one must first call SNESCreate()
as
follows:
SNESCreate(MPI_Comm comm,SNES *snes);
The user must then set routines for evaluating the residual function (3) and its associated Jacobian matrix, as discussed in the following sections.
To choose a nonlinear solution method, the user can either call
SNESSetType(SNES snes,SNESType method);
or use the option snes_type <method>
, where details regarding the
available methods are presented in The Nonlinear Solvers. The
application code can take complete control of the linear and nonlinear
techniques used in the Newtonlike method by calling
SNESSetFromOptions(snes);
This routine provides an interface to the PETSc options database, so
that at runtime the user can select a particular nonlinear solver, set
various parameters and customized routines (e.g., specialized line
search variants), prescribe the convergence tolerance, and set
monitoring routines. With this routine the user can also control all
linear solver options in the KSP
, and PC
modules, as discussed
in KSP: Linear System Solvers.
After having set these routines and options, the user solves the problem by calling
where x
should be initialized to the initial guess before calling and contains the solution on return.
In particular, to employ an initial guess of
zero, the user should explicitly set this vector to zero by calling
VecZeroEntries(x)
. Finally, after solving the nonlinear system (or several
systems), the user should destroy the SNES
context with
SNESDestroy(SNES *snes);
Nonlinear Function Evaluation#
When solving a system of nonlinear equations, the user must provide a a residual function (3), which is set using
SNESSetFunction(SNES snes,Vec f,PetscErrorCode (*FormFunction)(SNES snes,Vec x,Vec f,void *ctx),void *ctx);
The argument f
is an optional vector for storing the solution; pass NULL
to have the SNES
allocate it for you.
The argument ctx
is an optional userdefined context, which can
store any private, applicationspecific data required by the function
evaluation routine; NULL
should be used if such information is not
needed. In C and C++, a userdefined context is merely a structure in
which various objects can be stashed; in Fortran a user context can be
an integer array that contains both parameters and pointers to PETSc
objects.
SNES Tutorial ex5
and
SNES Tutorial ex5f
give examples of userdefined application contexts in C and Fortran,
respectively.
Jacobian Evaluation#
The user must also specify a routine to form some approximation of the
Jacobian matrix, A
, at the current iterate, x
, as is typically
done with
SNESSetJacobian(SNES snes,Mat Amat,Mat Pmat,PetscErrorCode (*FormJacobian)(SNES snes,Vec x,Mat A,Mat B,void *ctx),void *ctx);
The arguments of the routine FormJacobian()
are the current iterate,
x
; the (approximate) Jacobian matrix, Amat
; the matrix from
which the preconditioner is constructed, Pmat
(which is usually the
same as Amat
); and an optional userdefined Jacobian context,
ctx
, for applicationspecific data. Note that the SNES
solvers
are all datastructure neutral, so the full range of PETSc matrix
formats (including “matrixfree” methods) can be used.
Matrices discusses information regarding
available matrix formats and options, while MatrixFree Methods focuses on matrixfree methods in
SNES
. We briefly touch on a few details of matrix usage that are
particularly important for efficient use of the nonlinear solvers.
A common usage paradigm is to assemble the problem Jacobian in the
preconditioner storage B
, rather than A
. In the case where they
are identical, as in many simulations, this makes no difference.
However, it allows us to check the analytic Jacobian we construct in
FormJacobian()
by passing the snes_mf_operator
flag. This
causes PETSc to approximate the Jacobian using finite differencing of
the function evaluation (discussed in Finite Difference Jacobian Approximations),
and the analytic Jacobian becomes merely the preconditioner. Even if the
analytic Jacobian is incorrect, it is likely that the finite difference
approximation will converge, and thus this is an excellent method to
verify the analytic Jacobian. Moreover, if the analytic Jacobian is
incomplete (some terms are missing or approximate),
snes_mf_operator
may be used to obtain the exact solution, where
the Jacobian approximation has been transferred to the preconditioner.
One such approximate Jacobian comes from “Picard linearization” which writes the nonlinear system as
where \(\mathbf{A}(\mathbf{x})\) usually contains the lowerderivative parts of the equation. For example, the nonlinear diffusion problem
would be linearized as
Usually this linearization is simpler to implement than Newton and the
linear problems are somewhat easier to solve. In addition to using
snes_mf_operator
with this approximation to the Jacobian, the
Picard iterative procedure can be performed by defining \(\mathbf{J}(\mathbf{x})\)
to be \(\mathbf{A}(\mathbf{x})\). Sometimes this iteration exhibits better global
convergence than Newton linearization.
During successive calls to FormJacobian()
, the user can either
insert new matrix contexts or reuse old ones, depending on the
application requirements. For many sparse matrix formats, reusing the
old space (and merely changing the matrix elements) is more efficient;
however, if the matrix structure completely changes, creating an
entirely new matrix context may be preferable. Upon subsequent calls to
the FormJacobian()
routine, the user may wish to reinitialize the
matrix entries to zero by calling MatZeroEntries()
. See
Other Matrix Operations for details on the reuse of the matrix
context.
The directory $PETSC_DIR/src/snes/tutorials
provides a variety of
examples.
The Nonlinear Solvers#
As summarized in Table PETSc Nonlinear Solvers, SNES
includes
several Newtonlike nonlinear solvers based on line search techniques
and trust region methods. Also provided are several nonlinear Krylov
methods, as well as nonlinear methods involving decompositions of the
problem.
Each solver may have associated with it a set of options, which can be
set with routines and options database commands provided for this
purpose. A complete list can be found by consulting the manual pages or
by running a program with the help
option; we discuss just a few in
the sections below.
Method 
SNESType 
Options Name 
Default Line Search 

Line Search Newton 


Trust region Newton 

— 

Nonlinear Richardson 


Nonlinear CG 


Nonlinear GMRES 


QuasiNewton 


Full Approximation Scheme 

— 

Nonlinear ASM 

– 

ASPIN 


Nonlinear GaussSeidel 

– 

Anderson Mixing 

– 

Newton with constraints (1) 


Newton with constraints (2) 


Multistage Smoothers 

– 

Composite 

– 

Linear solve only 

– 

Python Shell 


– 
Shell (userdefined) 

– 
Line Search Newton#
The method SNESNEWTONLS
(snes_type newtonls
) provides a
line search Newton method for solving systems of nonlinear equations. By
default, this technique employs cubic backtracking
[DennisJrS83]. Alternative line search techniques are
listed in Table PETSc Line Search Methods.
Line Search 
SNESLineSearchType 
Options Name 

Backtracking 


(damped) step 


L2norm Minimization 


Critical point 


Shell 

Every SNES
has a line search context of type SNESLineSearch
that
may be retrieved using
SNESGetLineSearch(SNES snes,SNESLineSearch *ls);.
There are several default options for the line searches. The order of
polynomial approximation may be set with snes_linesearch_order
or
SNESLineSearchSetOrder(SNESLineSearch ls, PetscInt order);
for instance, 2 for quadratic or 3 for cubic. Sometimes, it may not be
necessary to monitor the progress of the nonlinear iteration. In this
case, snes_linesearch_norms
or
SNESLineSearchSetComputeNorms(SNESLineSearch ls,PetscBool norms);
may be used to turn off function, step, and solution norm computation at the end of the linesearch.
The default line search for the line search Newton method,
SNESLINESEARCHBT
involves several parameters, which are set to
defaults that are reasonable for many applications. The user can
override the defaults by using the following options:
snes_linesearch_alpha <alpha>
snes_linesearch_maxstep <max>
snes_linesearch_minlambda <tol>
Besides the backtracking linesearch, there are SNESLINESEARCHL2
,
which uses a polynomial secant minimization of \(F(x)_2\), and
SNESLINESEARCHCP
, which minimizes \(F(x) \cdot Y\) where
\(Y\) is the search direction. These are both potentially iterative
line searches, which may be used to find a betterfitted steplength in
the case where a single secant search is not sufficient. The number of
iterations may be set with snes_linesearch_max_it
. In addition, the
convergence criteria of the iterative line searches may be set using
function tolerances snes_linesearch_rtol
and
snes_linesearch_atol
, and steplength tolerance
snes_linesearch_ltol
.
Custom line search types may either be defined using
SNESLineSearchShell
, or by creating a custom user line search type
in the model of the preexisting ones and register it using
SNESLineSearchRegister(const char sname[],PetscErrorCode (*function)(SNESLineSearch));.
Trust Region Methods#
The trust region method in SNES
for solving systems of nonlinear
equations, SNESNEWTONTR
(snes_type newtontr
), is taken from the
MINPACK project [MoreSGH84]. Several parameters can be
set to control the variation of the trust region size during the
solution process. In particular, the user can control the initial trust
region radius, computed by
by setting \(\Delta_0\) via the option snes_tr_delta0 <delta0>
.
Nonlinear Krylov Methods#
A number of nonlinear Krylov methods are provided, including Nonlinear Richardson, conjugate gradient, GMRES, and Anderson Mixing. These methods are described individually below. They are all instrumental to PETSc’s nonlinear preconditioning.
Nonlinear Richardson. The nonlinear Richardson iteration merely takes the form of a line searchdamped fixedpoint iteration of the form
where the default linesearch is SNESLINESEARCHL2
. This simple solver
is mostly useful as a nonlinear smoother, or to provide line search
stabilization to an inner method.
Nonlinear Conjugate Gradients. Nonlinear CG is equivalent to linear
CG, but with the steplength determined by line search
(SNESLINESEARCHCP
by default). Five variants (FletcherReed,
HestenesSteifel, PolakRibierePolyak, DaiYuan, and Conjugate Descent)
are implemented in PETSc and may be chosen using
SNESNCGSetType(SNES snes, SNESNCGType btype);
Anderson Mixing and Nonlinear GMRES Methods. Nonlinear GMRES and Anderson Mixing methods combine the last \(m\) iterates, plus a new fixedpoint iteration iterate, into a residualminimizing new iterate.
QuasiNewton Methods#
QuasiNewton methods store iterative rankone updates to the Jacobian
instead of computing it directly. Three limitedmemory quasiNewton
methods are provided, LBFGS, which are described in
Table PETSc quasiNewton solvers. These all are encapsulated under
snes_type qn
and may be changed with snes_qn_type
. The default
is LBFGS, which provides symmetric updates to an approximate Jacobian.
This iteration is similar to the line search Newton methods.
QN Method 

Options Name 
Default Line Search 

LBFGS 



“Good” Broyden 



“Bad” Broyden 


One may also control the form of the initial Jacobian approximation with
SNESQNSetScaleType(SNES snes, SNESQNScaleType stype);
and the restart type with
SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype);
The Full Approximation Scheme#
The Full Approximation Scheme is a nonlinear multigrid correction. At
each level, there is a recursive cycle control SNES
instance, and
either one or two nonlinear solvers as smoothers (up and down). Problems
set up using the SNES
DMDA
interface are automatically
coarsened. FAS differs slightly from PCMG
, in that the hierarchy is
constructed recursively. However, much of the interface is a onetoone
map. We describe the “get” operations here, and it can be assumed that
each has a corresponding “set” operation. For instance, the number of
levels in the hierarchy may be retrieved using
SNESFASGetLevels(SNES snes, PetscInt *levels);
There are four SNESFAS
cycle types, SNES_FAS_MULTIPLICATIVE
,
SNES_FAS_ADDITIVE
, SNES_FAS_FULL
, and SNES_FAS_KASKADE
. The
type may be set with
SNESFASSetType(SNES snes,SNESFASType fastype);.
and the cycle type, 1 for V, 2 for W, may be set with
SNESFASSetCycles(SNES snes, PetscInt cycles);.
Much like the interface to PCMG
described in Multigrid Preconditioners, there are interfaces to recover the
various levels’ cycles and smoothers. The level smoothers may be
accessed with
SNESFASGetSmoother(SNES snes, PetscInt level, SNES *smooth);
SNESFASGetSmootherUp(SNES snes, PetscInt level, SNES *smooth);
SNESFASGetSmootherDown(SNES snes, PetscInt level, SNES *smooth);
and the level cycles with
SNESFASGetCycleSNES(SNES snes,PetscInt level,SNES *lsnes);.
Also akin to PCMG
, the restriction and prolongation at a level may
be acquired with
SNESFASGetInterpolation(SNES snes, PetscInt level, Mat *mat);
SNESFASGetRestriction(SNES snes, PetscInt level, Mat *mat);
In addition, FAS requires special restriction for solutionlike variables, called injection. This may be set with
SNESFASGetInjection(SNES snes, PetscInt level, Mat *mat);.
The coarse solve context may be acquired with
SNESFASGetCoarseSolve(SNES snes, SNES *smooth);
Nonlinear Additive Schwarz#
Nonlinear Additive Schwarz methods (NASM) take a number of local nonlinear subproblems, solves them independently in parallel, and combines those solutions into a new approximate solution.
SNESNASMSetSubdomains(SNES snes,PetscInt n,SNES subsnes[],VecScatter iscatter[],VecScatter oscatter[],VecScatter gscatter[]);
allows for the user to create these local subdomains. Problems set up
using the SNES
DMDA
interface are automatically decomposed. To
begin, the type of subdomain updates to the whole solution are limited
to two types borrowed from PCASM
: PC_ASM_BASIC
, in which the
overlapping updates added. PC_ASM_RESTRICT
updates in a
nonoverlapping fashion. This may be set with
SNESNASMSetType(SNES snes,PCASMType type);.
SNESASPIN
is a helper SNES
type that sets up a nonlinearly
preconditioned Newton’s method using NASM as the preconditioner.
General Options#
This section discusses options and routines that apply to all SNES
solvers and problem classes. In particular, we focus on convergence
tests, monitoring routines, and tools for checking derivative
computations.
Convergence Tests#
Convergence of the nonlinear solvers can be detected in a variety of
ways; the user can even specify a customized test, as discussed below.
Most of the nonlinear solvers use SNESConvergenceTestDefault()
,
however, SNESNEWTONTR
uses a methodspecific additional convergence
test as well. The convergence tests involves several parameters, which
are set by default to values that should be reasonable for a wide range
of problems. The user can customize the parameters to the problem at
hand by using some of the following routines and options.
One method of convergence testing is to declare convergence when the
norm of the change in the solution between successive iterations is less
than some tolerance, stol
. Convergence can also be determined based
on the norm of the function. Such a test can use either the absolute
size of the norm, atol
, or its relative decrease, rtol
, from an
initial guess. The following routine sets these parameters, which are
used in many of the default SNES
convergence tests:
This routine also sets the maximum numbers of allowable nonlinear
iterations, its
, and function evaluations, fcts
. The
corresponding options database commands for setting these parameters are:
snes_atol <atol>
snes_rtol <rtol>
snes_stol <stol>
snes_max_it <its>
snes_max_funcs <fcts>
A related routine is SNESGetTolerances()
.
Convergence tests for trust regions methods often use an additional
parameter that indicates the minimum allowable trust region radius. The
user can set this parameter with the option snes_trtol <trtol>
or
with the routine
SNESSetTrustRegionTolerance(SNES snes,PetscReal trtol);
Users can set their own customized convergence tests in SNES
by
using the command
SNESSetConvergenceTest(SNES snes,PetscErrorCode (*test)(SNES snes,PetscInt it,PetscReal xnorm, PetscReal gnorm,PetscReal f,SNESConvergedReason reason, void *cctx),void *cctx,PetscErrorCode (*destroy)(void *cctx));
The final argument of the convergence test routine, cctx
, denotes an
optional userdefined context for private data. When solving systems of
nonlinear equations, the arguments xnorm
, gnorm
, and f
are
the current iterate norm, current step norm, and function norm,
respectively. SNESConvergedReason
should be set positive for
convergence and negative for divergence. See include/petscsnes.h
for
a list of values for SNESConvergedReason
.
Convergence Monitoring#
By default the SNES
solvers run silently without displaying
information about the iterations. The user can initiate monitoring with
the command
SNESMonitorSet(SNES snes,PetscErrorCode (*mon)(SNES,PetscInt its,PetscReal norm,void* mctx),void *mctx,PetscErrorCode (*monitordestroy)(void**));
The routine, mon
, indicates a userdefined monitoring routine, where
its
and mctx
respectively denote the iteration number and an
optional userdefined context for private data for the monitor routine.
The argument norm
is the function norm.
The routine set by SNESMonitorSet()
is called once after every
successful step computation within the nonlinear solver. Hence, the user
can employ this routine for any applicationspecific computations that
should be done after the solution update. The option snes_monitor
activates the default SNES
monitor routine,
SNESMonitorDefault()
, while snes_monitor_lg_residualnorm
draws
a simple line graph of the residual norm’s convergence.
One can cancel hardwired monitoring routines for SNES
at runtime
with snes_monitor_cancel
.
As the Newton method converges so that the residual norm is small, say
\(10^{10}\), many of the final digits printed with the
snes_monitor
option are meaningless. Worse, they are different on
different machines; due to different roundoff rules used by, say, the
IBM RS6000 and the Sun SPARC. This makes testing between different
machines difficult. The option snes_monitor_short
causes PETSc to
print fewer of the digits of the residual norm as it gets smaller; thus
on most of the machines it will always print the same numbers making
crossprocess testing easier.
The routines
SNESGetSolution(SNES snes,Vec *x);
SNESGetFunction(SNES snes,Vec *r,void *ctx,int(**func)(SNES,Vec,Vec,void*));
return the solution vector and function vector from a SNES
context.
These routines are useful, for instance, if the convergence test
requires some property of the solution or function other than those
passed with routine arguments.
Checking Accuracy of Derivatives#
Since handcoding routines for Jacobian matrix evaluation can be error
prone, SNES
provides easytouse support for checking these matrices
against finite difference versions. In the simplest form of comparison,
users can employ the option snes_test_jacobian
to compare the
matrices at several points. Although not exhaustive, this test will
generally catch obvious problems. One can compare the elements of the
two matrices by using the option snes_test_jacobian_view
, which
causes the two matrices to be printed to the screen.
Another means for verifying the correctness of a code for Jacobian
computation is running the problem with either the finite difference or
matrixfree variant, snes_fd
or snes_mf
; see Finite Difference Jacobian Approximations or MatrixFree Methods.
If a
problem converges well with these matrix approximations but not with a
userprovided routine, the problem probably lies with the handcoded
matrix. See the note in Jacobian Evaluation about
assembling your Jabobian in the “preconditioner” slot Pmat
.
The correctness of user provided MATSHELL
Jacobians in general can be
checked with MatShellTestMultTranspose()
and MatShellTestMult()
.
The correctness of user provided MATSHELL
Jacobians via TSSetRHSJacobian()
can be checked with TSRHSJacobianTestTranspose()
and TSRHSJacobianTest()
that check the correction of the matrixtranspose vector product and the
matrixproduct. From the command line, these can be checked with
ts_rhs_jacobian_test_mult_transpose
mat_shell_test_mult_transpose_view
ts_rhs_jacobian_test_mult
mat_shell_test_mult_view
Inexact Newtonlike Methods#
Since exact solution of the linear Newton systems within (4)
at each iteration can be costly, modifications
are often introduced that significantly reduce these expenses and yet
retain the rapid convergence of Newton’s method. Inexact or truncated
Newton techniques approximately solve the linear systems using an
iterative scheme. In comparison with using direct methods for solving
the Newton systems, iterative methods have the virtue of requiring
little space for matrix storage and potentially saving significant
computational work. Within the class of inexact Newton methods, of
particular interest are NewtonKrylov methods, where the subsidiary
iterative technique for solving the Newton system is chosen from the
class of Krylov subspace projection methods. Note that at runtime the
user can set any of the linear solver options discussed in KSP: Linear System Solvers,
such as ksp_type <ksp_method>
and
pc_type <pc_method>
, to set the Krylov subspace and preconditioner
methods.
Two levels of iterations occur for the inexact techniques, where during each global or outer Newton iteration a sequence of subsidiary inner iterations of a linear solver is performed. Appropriate control of the accuracy to which the subsidiary iterative method solves the Newton system at each global iteration is critical, since these inner iterations determine the asymptotic convergence rate for inexact Newton techniques. While the Newton systems must be solved well enough to retain fast local convergence of the Newton’s iterates, use of excessive inner iterations, particularly when \(\ \mathbf{x}_k  \mathbf{x}_* \\) is large, is neither necessary nor economical. Thus, the number of required inner iterations typically increases as the Newton process progresses, so that the truncated iterates approach the true Newton iterates.
A sequence of nonnegative numbers \(\{\eta_k\}\) can be used to indicate the variable convergence criterion. In this case, when solving a system of nonlinear equations, the update step of the Newton process remains unchanged, and direct solution of the linear system is replaced by iteration on the system until the residuals
satisfy
Here \(\mathbf{x}_0\) is an initial approximation of the solution, and \(\ \cdot \\) denotes an arbitrary norm in \(\Re^n\) .
By default a constant relative convergence tolerance is used for solving
the subsidiary linear systems within the Newtonlike methods of
SNES
. When solving a system of nonlinear equations, one can instead
employ the techniques of Eisenstat and Walker [EW96]
to compute \(\eta_k\) at each step of the nonlinear solver by using
the option snes_ksp_ew
. In addition, by adding one’s own
KSP
convergence test (see Convergence Tests), one can easily create one’s own,
problemdependent, inner convergence tests.
MatrixFree Methods#
The SNES
class fully supports matrixfree methods. The matrices
specified in the Jacobian evaluation routine need not be conventional
matrices; instead, they can point to the data required to implement a
particular matrixfree method. The matrixfree variant is allowed only
when the linear systems are solved by an iterative method in combination
with no preconditioning (PCNONE
or pc_type
none
), a
userprovided preconditioner matrix, or a userprovided preconditioner
shell (PCSHELL
, discussed in Preconditioners); that
is, obviously matrixfree methods cannot be used with a direct solver,
approximate factorization, or other preconditioner which requires access
to explicit matrix entries.
The user can create a matrixfree context for use within SNES
with
the routine
MatCreateSNESMF(SNES snes,Mat *mat);
This routine creates the data structures needed for the matrixvector
products that arise within Krylov space iterative
methods [BS90].
The default SNES
matrixfree approximations can also be invoked with the command
snes_mf
. Or, one can retain the userprovided Jacobian
preconditioner, but replace the userprovided Jacobian matrix with the
default matrix free variant with the option snes_mf_operator
.
MatCreateSNESMF()
uses
MatCreateMFFD(Vec x, Mat *mat);
which can also be used directly for users who need a matrixfree matrix but are not using SNES
.
The user can set one parameter to control the Jacobianvector product approximation with the command
MatMFFDSetFunctionError(Mat mat,PetscReal rerror);
The parameter rerror
should be set to the square root of the
relative error in the function evaluations, \(e_{rel}\); the default
is the square root of machine epsilon (about \(10^{8}\) in double
precision), which assumes that the functions are evaluated to full
floatingpoint precision accuracy. This parameter can also be set from
the options database with mat_mffd_err <err>
In addition, PETSc provides ways to register new routines to compute
the differencing parameter (\(h\)); see the manual page for
MatMFFDSetType()
and MatMFFDRegister()
. We currently provide two
default routines accessible via mat_mffd_type <ds or wp>
. For
the default approach there is one “tuning” parameter, set with
MatMFFDDSSetUmin(Mat mat,PetscReal umin);
This parameter, umin
(or \(u_{min}\)), is a bit involved; its
default is \(10^{6}\) . Its command line form is mat_mffd_umin <umin>
.
The Jacobianvector product is approximated via the formula
where \(h\) is computed via
This approach is taken from Brown and Saad [BS90]. The second approach, taken from Walker and Pernice, [PW98], computes \(h\) via
This has no tunable parameters, but note that inside the nonlinear solve for the entire linear iterative process \(u\) does not change hence \(\sqrt{1 + u}\) need be computed only once. This information may be set with the options
MatMFFDWPSetComputeNormU(Mat mat,PetscBool );
or mat_mffd_compute_normu <true or false>
. This information is used
to eliminate the redundant computation of these parameters, therefore
reducing the number of collective operations and improving the
efficiency of the application code. This takes place automatically for the PETSc GMRES solver with left preconditioning.
It is also possible to monitor the differencing parameters h that are computed via the routines
MatMFFDSetHHistory(Mat,PetscScalar *,int);
MatMFFDResetHHistory(Mat,PetscScalar *,int);
MatMFFDGetH(Mat,PetscScalar *);
We include an explicit example of using matrixfree methods in ex3.c.
Note that by using the option snes_mf
one can
easily convert any SNES
code to use a matrixfree NewtonKrylov
method without a preconditioner. As shown in this example,
SNESSetFromOptions()
must be called after SNESSetJacobian()
to
enable runtime switching between the userspecified Jacobian and the
default SNES
matrixfree form.
Listing: src/snes/tutorials/ex3.c
static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
This example employs a userdefined monitoring routine and optionally a userdefined\n\
routine to check candidate iterates produced by line search routines.\n\
The command line options include:\n\
pre_check_iterates : activate checking of iterates\n\
post_check_iterates : activate checking of iterates\n\
check_tol <tol>: set tolerance for iterate checking\n\
user_precond : activate a (trivial) userdefined preconditioner\n\n";
/*
Include "petscdm.h" so that we can use data management objects (DMs)
Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
Include "petscsnes.h" so that we can use SNES solvers. Note that this
file automatically includes:
petscsys.h  base PETSc routines
petscvec.h  vectors
petscmat.h  matrices
petscis.h  index sets
petscksp.h  Krylov subspace methods
petscviewer.h  viewers
petscpc.h  preconditioners
petscksp.h  linear solvers
*/
#include <petscdm.h>
#include <petscdmda.h>
#include <petscsnes.h>
/*
Userdefined routines.
*/
PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
PetscErrorCode FormInitialGuess(Vec);
PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
/*
Userdefined application context
*/
typedef struct {
DM da; /* distributed array */
Vec F; /* righthandside of PDE */
PetscMPIInt rank; /* rank of processor */
PetscMPIInt size; /* size of communicator */
PetscReal h; /* mesh spacing */
PetscBool sjerr; /* if or not to test jacobian domain error */
} ApplicationCtx;
/*
Userdefined context for monitoring
*/
typedef struct {
PetscViewer viewer;
} MonitorCtx;
/*
Userdefined context for checking candidate iterates that are
determined by line search methods
*/
typedef struct {
Vec last_step; /* previous iterate */
PetscReal tolerance; /* tolerance for changes between successive iterates */
ApplicationCtx *user;
} StepCheckCtx;
typedef struct {
PetscInt its0; /* num of prevous outer KSP iterations */
} SetSubKSPCtx;
int main(int argc,char **argv)
{
SNES snes; /* SNES context */
SNESLineSearch linesearch; /* SNESLineSearch context */
Mat J; /* Jacobian matrix */
ApplicationCtx ctx; /* userdefined context */
Vec x,r,U,F; /* vectors */
MonitorCtx monP; /* monitoring context */
StepCheckCtx checkP; /* stepchecking context */
SetSubKSPCtx checkP1;
PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
PetscScalar xp,*FF,*UU,none = 1.0;
PetscInt its,N = 5,i,maxit,maxf,xs,xm;
PetscReal abstol,rtol,stol,norm;
PetscBool flg,viewinitial = PETSC_FALSE;
PetscCall(PetscInitialize(&argc,&argv,(char*)0,help));
PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank));
PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size));
PetscCall(PetscOptionsGetInt(NULL,NULL,"n",&N,NULL));
ctx.h = 1.0/(N1);
ctx.sjerr = PETSC_FALSE;
PetscCall(PetscOptionsGetBool(NULL,NULL,"test_jacobian_domain_error",&ctx.sjerr,NULL));
PetscCall(PetscOptionsGetBool(NULL,NULL,"view_initial",&viewinitial,NULL));
/*                                  
Create nonlinear solver context
                                  */
PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes));
/*                                   
Create vector data structures; set function evaluation routine
                                   */
/*
Create distributed array (DMDA) to manage parallel grid and vectors
*/
PetscCall(DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da));
PetscCall(DMSetFromOptions(ctx.da));
PetscCall(DMSetUp(ctx.da));
/*
Extract global and local vectors from DMDA; then duplicate for remaining
vectors that are the same types
*/
PetscCall(DMCreateGlobalVector(ctx.da,&x));
PetscCall(VecDuplicate(x,&r));
PetscCall(VecDuplicate(x,&F)); ctx.F = F;
PetscCall(VecDuplicate(x,&U));
/*
Set function evaluation routine and vector. Whenever the nonlinear
solver needs to compute the nonlinear function, it will call this
routine.
 Note that the final routine argument is the userdefined
context that provides applicationspecific data for the
function evaluation routine.
*/
PetscCall(SNESSetFunction(snes,r,FormFunction,&ctx));
/*                                   
Create matrix data structure; set Jacobian evaluation routine
                                   */
PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N));
PetscCall(MatSetFromOptions(J));
PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
PetscCall(MatMPIAIJSetPreallocation(J,3,NULL,3,NULL));
/*
Set Jacobian matrix data structure and default Jacobian evaluation
routine. Whenever the nonlinear solver needs to compute the
Jacobian matrix, it will call this routine.
 Note that the final routine argument is the userdefined
context that provides applicationspecific data for the
Jacobian evaluation routine.
*/
PetscCall(SNESSetJacobian(snes,J,J,FormJacobian,&ctx));
/*
Optionally allow userprovided preconditioner
*/
PetscCall(PetscOptionsHasName(NULL,NULL,"user_precond",&flg));
if (flg) {
KSP ksp;
PC pc;
PetscCall(SNESGetKSP(snes,&ksp));
PetscCall(KSPGetPC(ksp,&pc));
PetscCall(PCSetType(pc,PCSHELL));
PetscCall(PCShellSetApply(pc,MatrixFreePreconditioner));
}
/*                                  
Customize nonlinear solver; set runtime options
                                   */
/*
Set an optional userdefined monitoring routine
*/
PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer));
PetscCall(SNESMonitorSet(snes,Monitor,&monP,0));
/*
Set names for some vectors to facilitate monitoring (optional)
*/
PetscCall(PetscObjectSetName((PetscObject)x,"Approximate Solution"));
PetscCall(PetscObjectSetName((PetscObject)U,"Exact Solution"));
/*
Set SNES/KSP/KSP/PC runtime options, e.g.,
snes_view snes_monitor ksp_type <ksp> pc_type <pc>
*/
PetscCall(SNESSetFromOptions(snes));
/*
Set an optional userdefined routine to check the validity of candidate
iterates that are determined by line search methods
*/
PetscCall(SNESGetLineSearch(snes, &linesearch));
PetscCall(PetscOptionsHasName(NULL,NULL,"post_check_iterates",&post_check));
if (post_check) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n"));
PetscCall(SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP));
PetscCall(VecDuplicate(x,&(checkP.last_step)));
checkP.tolerance = 1.0;
checkP.user = &ctx;
PetscCall(PetscOptionsGetReal(NULL,NULL,"check_tol",&checkP.tolerance,NULL));
}
PetscCall(PetscOptionsHasName(NULL,NULL,"post_setsubksp",&post_setsubksp));
if (post_setsubksp) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n"));
PetscCall(SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1));
}
PetscCall(PetscOptionsHasName(NULL,NULL,"pre_check_iterates",&pre_check));
if (pre_check) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n"));
PetscCall(SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP));
}
/*
Print parameters used for convergence testing (optional) ... just
to demonstrate this routine; this information is also printed with
the option snes_view
*/
PetscCall(SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf));
/*                                  
Initialize application:
Store righthandside of PDE and exact solution
                                   */
/*
Get local grid boundaries (for 1dimensional DMDA):
xs, xm  starting grid index, width of local grid (no ghost points)
*/
PetscCall(DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL));
/*
Get pointers to vector data
*/
PetscCall(DMDAVecGetArray(ctx.da,F,&FF));
PetscCall(DMDAVecGetArray(ctx.da,U,&UU));
/*
Compute local vector entries
*/
xp = ctx.h*xs;
for (i=xs; i<xs+xm; i++) {
FF[i] = 6.0*xp + PetscPowScalar(xp+1.e12,6.0); /* +1.e12 is to prevent 0^6 */
UU[i] = xp*xp*xp;
xp += ctx.h;
}
/*
Restore vectors
*/
PetscCall(DMDAVecRestoreArray(ctx.da,F,&FF));
PetscCall(DMDAVecRestoreArray(ctx.da,U,&UU));
if (viewinitial) {
PetscCall(VecView(U,0));
PetscCall(VecView(F,0));
}
/*                                  
Evaluate initial guess; then solve nonlinear system
                                   */
/*
Note: The user should initialize the vector, x, with the initial guess
for the nonlinear solver prior to calling SNESSolve(). In particular,
to employ an initial guess of zero, the user should explicitly set
this vector to zero by calling VecSet().
*/
PetscCall(FormInitialGuess(x));
PetscCall(SNESSolve(snes,NULL,x));
PetscCall(SNESGetIterationNumber(snes,&its));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its));
/*                                  
Check solution and clean up
                                   */
/*
Check the error
*/
PetscCall(VecAXPY(x,none,U));
PetscCall(VecNorm(x,NORM_2,&norm));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its));
if (ctx.sjerr) {
SNESType snestype;
PetscCall(SNESGetType(snes,&snestype));
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"SNES Type %s\n",snestype));
}
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
PetscCall(PetscViewerDestroy(&monP.viewer));
if (post_check) PetscCall(VecDestroy(&checkP.last_step));
PetscCall(VecDestroy(&x));
PetscCall(VecDestroy(&r));
PetscCall(VecDestroy(&U));
PetscCall(VecDestroy(&F));
PetscCall(MatDestroy(&J));
PetscCall(SNESDestroy(&snes));
PetscCall(DMDestroy(&ctx.da));
PetscCall(PetscFinalize());
return 0;
}
/*  */
/*
FormInitialGuess  Computes initial guess.
Input/Output Parameter:
. x  the solution vector
*/
PetscErrorCode FormInitialGuess(Vec x)
{
PetscScalar pfive = .50;
PetscFunctionBeginUser;
PetscCall(VecSet(x,pfive));
PetscFunctionReturn(0);
}
/*  */
/*
FormFunction  Evaluates nonlinear function, F(x).
Input Parameters:
. snes  the SNES context
. x  input vector
. ctx  optional userdefined context, as set by SNESSetFunction()
Output Parameter:
. f  function vector
Note:
The userdefined context can contain any applicationspecific
data needed for the function evaluation.
*/
PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
{
ApplicationCtx *user = (ApplicationCtx*) ctx;
DM da = user>da;
PetscScalar *xx,*ff,*FF,d;
PetscInt i,M,xs,xm;
Vec xlocal;
PetscFunctionBeginUser;
PetscCall(DMGetLocalVector(da,&xlocal));
/*
Scatter ghost points to local vector, using the 2step process
DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
By placing code between these two statements, computations can
be done while messages are in transition.
*/
PetscCall(DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal));
PetscCall(DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal));
/*
Get pointers to vector data.
 The vector xlocal includes ghost point; the vectors x and f do
NOT include ghost points.
 Using DMDAVecGetArray() allows accessing the values using global ordering
*/
PetscCall(DMDAVecGetArray(da,xlocal,&xx));
PetscCall(DMDAVecGetArray(da,f,&ff));
PetscCall(DMDAVecGetArray(da,user>F,&FF));
/*
Get local grid boundaries (for 1dimensional DMDA):
xs, xm  starting grid index, width of local grid (no ghost points)
*/
PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL));
PetscCall(DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL));
/*
Set function values for boundary points; define local interior grid point range:
xsi  starting interior grid index
xei  ending interior grid index
*/
if (xs == 0) { /* left boundary */
ff[0] = xx[0];
xs++;xm;
}
if (xs+xm == M) { /* right boundary */
ff[xs+xm1] = xx[xs+xm1]  1.0;
xm;
}
/*
Compute function over locally owned part of the grid (interior points only)
*/
d = 1.0/(user>h*user>h);
for (i=xs; i<xs+xm; i++) ff[i] = d*(xx[i1]  2.0*xx[i] + xx[i+1]) + xx[i]*xx[i]  FF[i];
/*
Restore vectors
*/
PetscCall(DMDAVecRestoreArray(da,xlocal,&xx));
PetscCall(DMDAVecRestoreArray(da,f,&ff));
PetscCall(DMDAVecRestoreArray(da,user>F,&FF));
PetscCall(DMRestoreLocalVector(da,&xlocal));
PetscFunctionReturn(0);
}
/*  */
/*
FormJacobian  Evaluates Jacobian matrix.
Input Parameters:
. snes  the SNES context
. x  input vector
. dummy  optional userdefined context (not used here)
Output Parameters:
. jac  Jacobian matrix
. B  optionally different preconditioning matrix
. flag  flag indicating matrix structure
*/
PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
{
ApplicationCtx *user = (ApplicationCtx*) ctx;
PetscScalar *xx,d,A[3];
PetscInt i,j[3],M,xs,xm;
DM da = user>da;
PetscFunctionBeginUser;
if (user>sjerr) {
PetscCall(SNESSetJacobianDomainError(snes));
PetscFunctionReturn(0);
}
/*
Get pointer to vector data
*/
PetscCall(DMDAVecGetArrayRead(da,x,&xx));
PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL));
/*
Get range of locally owned matrix
*/
PetscCall(DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL));
/*
Determine starting and ending local indices for interior grid points.
Set Jacobian entries for boundary points.
*/
if (xs == 0) { /* left boundary */
i = 0; A[0] = 1.0;
PetscCall(MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES));
xs++;xm;
}
if (xs+xm == M) { /* right boundary */
i = M1;
A[0] = 1.0;
PetscCall(MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES));
xm;
}
/*
Interior grid points
 Note that in this case we set all elements for a particular
row at once.
*/
d = 1.0/(user>h*user>h);
for (i=xs; i<xs+xm; i++) {
j[0] = i  1; j[1] = i; j[2] = i + 1;
A[0] = A[2] = d; A[1] = 2.0*d + 2.0*xx[i];
PetscCall(MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES));
}
/*
Assemble matrix, using the 2step process:
MatAssemblyBegin(), MatAssemblyEnd().
By placing code between these two statements, computations can be
done while messages are in transition.
Also, restore vector.
*/
PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
PetscCall(DMDAVecRestoreArrayRead(da,x,&xx));
PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
PetscFunctionReturn(0);
}
/*  */
/*
Monitor  Optional userdefined monitoring routine that views the
current iterate with an xwindow plot. Set by SNESMonitorSet().
Input Parameters:
snes  the SNES context
its  iteration number
norm  2norm function value (may be estimated)
ctx  optional userdefined context for private data for the
monitor routine, as set by SNESMonitorSet()
Note:
See the manpage for PetscViewerDrawOpen() for useful runtime options,
such as nox to deactivate all xwindow output.
*/
PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
{
MonitorCtx *monP = (MonitorCtx*) ctx;
Vec x;
PetscFunctionBeginUser;
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm));
PetscCall(SNESGetSolution(snes,&x));
PetscCall(VecView(x,monP>viewer));
PetscFunctionReturn(0);
}
/*  */
/*
PreCheck  Optional userdefined routine that checks the validity of
candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
Input Parameters:
snes  the SNES context
xcurrent  current solution
y  search direction and length
Output Parameters:
y  proposed step (search direction and length) (possibly changed)
changed_y  tells if the step has changed or not
*/
PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
{
PetscFunctionBeginUser;
*changed_y = PETSC_FALSE;
PetscFunctionReturn(0);
}
/*  */
/*
PostCheck  Optional userdefined routine that checks the validity of
candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
Input Parameters:
snes  the SNES context
ctx  optional userdefined context for private data for the
monitor routine, as set by SNESLineSearchSetPostCheck()
xcurrent  current solution
y  search direction and length
x  the new candidate iterate
Output Parameters:
y  proposed step (search direction and length) (possibly changed)
x  current iterate (possibly modified)
*/
PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
{
PetscInt i,iter,xs,xm;
StepCheckCtx *check;
ApplicationCtx *user;
PetscScalar *xa,*xa_last,tmp;
PetscReal rdiff;
DM da;
SNES snes;
PetscFunctionBeginUser;
*changed_x = PETSC_FALSE;
*changed_y = PETSC_FALSE;
PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
check = (StepCheckCtx*)ctx;
user = check>user;
PetscCall(SNESGetIterationNumber(snes,&iter));
/* iteration 1 indicates we are working on the second iteration */
if (iter > 0) {
da = user>da;
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check>tolerance));
/* Access local array data */
PetscCall(DMDAVecGetArray(da,check>last_step,&xa_last));
PetscCall(DMDAVecGetArray(da,x,&xa));
PetscCall(DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL));
/*
If we fail the userdefined check for validity of the candidate iterate,
then modify the iterate as we like. (Note that the particular modification
below is intended simply to demonstrate how to manipulate this data, not
as a meaningful or appropriate choice.)
*/
for (i=xs; i<xs+xm; i++) {
if (!PetscAbsScalar(xa[i])) rdiff = 2*check>tolerance;
else rdiff = PetscAbsScalar((xa[i]  xa_last[i])/xa[i]);
if (rdiff > check>tolerance) {
tmp = xa[i];
xa[i] = .5*(xa[i] + xa_last[i]);
*changed_x = PETSC_TRUE;
PetscCall(PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%g, x_last=%g, diff=%g, x_new=%g\n",i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i])));
}
}
PetscCall(DMDAVecRestoreArray(da,check>last_step,&xa_last));
PetscCall(DMDAVecRestoreArray(da,x,&xa));
}
PetscCall(VecCopy(x,check>last_step));
PetscFunctionReturn(0);
}
/*  */
/*
PostSetSubKSP  Optional userdefined routine that reset SubKSP options when hierarchical bjacobi PC is used
e.g,
mpiexec n 8 ./ex3 nox n 10000 ksp_type fgmres pc_type bjacobi pc_bjacobi_blocks 4 sub_ksp_type gmres sub_ksp_max_it 3 post_setsubksp sub_ksp_rtol 1.e16
Set by SNESLineSearchSetPostCheck().
Input Parameters:
linesearch  the LineSearch context
xcurrent  current solution
y  search direction and length
x  the new candidate iterate
Output Parameters:
y  proposed step (search direction and length) (possibly changed)
x  current iterate (possibly modified)
*/
PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
{
SetSubKSPCtx *check;
PetscInt iter,its,sub_its,maxit;
KSP ksp,sub_ksp,*sub_ksps;
PC pc;
PetscReal ksp_ratio;
SNES snes;
PetscFunctionBeginUser;
PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
check = (SetSubKSPCtx*)ctx;
PetscCall(SNESGetIterationNumber(snes,&iter));
PetscCall(SNESGetKSP(snes,&ksp));
PetscCall(KSPGetPC(ksp,&pc));
PetscCall(PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps));
sub_ksp = sub_ksps[0];
PetscCall(KSPGetIterationNumber(ksp,&its)); /* outer KSP iteration number */
PetscCall(KSPGetIterationNumber(sub_ksp,&sub_its)); /* inner KSP iteration number */
if (iter) {
PetscCall(PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %D %D, subksp_it %D\n",iter,check>its0,its,sub_its));
ksp_ratio = ((PetscReal)(its))/check>its0;
maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
if (maxit < 2) maxit = 2;
PetscCall(KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit));
PetscCall(PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %D\n\n",(double)ksp_ratio,maxit));
}
check>its0 = its; /* save current outer KSP iteration number */
PetscFunctionReturn(0);
}
/*  */
/*
MatrixFreePreconditioner  This routine demonstrates the use of a
userprovided preconditioner. This code implements just the null
preconditioner, which of course is not recommended for general use.
Input Parameters:
+ pc  preconditioner
 x  input vector
Output Parameter:
. y  preconditioned vector
*/
PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
{
PetscCall(VecCopy(x,y));
return 0;
}
Table Jacobian Options summarizes the various matrix situations
that SNES
supports. In particular, different linear system matrices
and preconditioning matrices are allowed, as well as both matrixfree
and applicationprovided preconditioners. If ex3.c is run with
the options snes_mf
and user_precond
then it uses a
matrixfree application of the matrixvector multiple and a user
provided matrix free Jacobian.
Matrix Use 
Conventional Matrix Formats 
Matrixfree versions 
Jacobian Matrix 
Create matrix with 
Create matrix with 
Preconditioning Matrix 
Create matrix with 
Use 
SNESSetJacobian()
or with a DM
variant such as DMDASNESSetJacobianLocal()
SNES also provides some less wellintegrated code to apply matrixfree finite differencing using an automatically computed measurement of the
noise of the functions. This can be selected with snes_mf_version 2
; it does not use MatCreateMFFD()
but has similar options that start with
snes_mf_
instead of mat_mffd_
. Note that this alternative prefix only works for version 2 differencing.
Finite Difference Jacobian Approximations#
PETSc provides some tools to help approximate the Jacobian matrices efficiently via finite differences. These tools are intended for use in certain situations where one is unable to compute Jacobian matrices analytically, and matrixfree methods do not work well without a preconditioner, due to very poor conditioning. The approximation requires several steps:
First, one colors the columns of the (not yet built) Jacobian matrix, so that columns of the same color do not share any common rows.
Next, one creates a
MatFDColoring
data structure that will be used later in actually computing the Jacobian.Finally, one tells the nonlinear solvers of
SNES
to use theSNESComputeJacobianDefaultColor()
routine to compute the Jacobians.
A code fragment that demonstrates this process is given below.
ISColoring iscoloring;
MatFDColoring fdcoloring;
MatColoring coloring;
/*
This initializes the nonzero structure of the Jacobian. This is artificial
because clearly if we had a routine to compute the Jacobian we wouldn't
need to use finite differences.
*/
FormJacobian(snes,x, &J, &J, &user);
/*
Color the matrix, i.e. determine groups of columns that share no common
rows. These columns in the Jacobian can all be computed simultaneously.
*/
MatColoringCreate(J, &coloring);
MatColoringSetType(coloring,MATCOLORINGSL);
MatColoringSetFromOptions(coloring);
MatColoringApply(coloring, &iscoloring);
MatColoringDestroy(&coloring);
/*
Create the data structure that SNESComputeJacobianDefaultColor() uses
to compute the actual Jacobians via finite differences.
*/
MatFDColoringCreate(J,iscoloring, &fdcoloring);
ISColoringDestroy(&iscoloring);
MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction, &user);
MatFDColoringSetFromOptions(fdcoloring);
/*
Tell SNES to use the routine SNESComputeJacobianDefaultColor()
to compute Jacobians.
*/
SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
Of course, we are cheating a bit. If we do not have an analytic formula
for computing the Jacobian, then how do we know what its nonzero
structure is so that it may be colored? Determining the structure is
problem dependent, but fortunately, for most structured grid problems
(the class of problems for which PETSc was originally designed) if one
knows the stencil used for the nonlinear function one can usually fairly
easily obtain an estimate of the location of nonzeros in the matrix.
This is harder in the unstructured case, but one typically knows where the nonzero entries are from the mesh topology and distribution of degrees of freedom.
If using DMPlex
(DMPlex: Unstructured Grids in PETSc) for unstructured meshes, the nonzero locations will be identified in DMCreateMatrix()
and the procedure above can be used.
Most external packages for unstructured meshes have similar functionality.
One need not necessarily use a MatColoring
object to determine a
coloring. For example, if a grid can be colored directly (without using
the associated matrix), then that coloring can be provided to
MatFDColoringCreate()
. Note that the user must always preset the
nonzero structure in the matrix regardless of which coloring routine is
used.
PETSc provides the following coloring algorithms, which can be selected using MatColoringSetType()
or via the command line argument mat_coloring_type
.
Algorithm 

Parallel 


smallestlast [MoreSGH84] 

No 

largestfirst [MoreSGH84] 

No 

incidencedegree [MoreSGH84] 

No 

JonesPlassmann [JP93] 

Yes 

Greedy 

Yes 

Natural (1 color per column) 

Yes 

Power (\(A^k\) followed by 1coloring) 

Yes 
As for the matrixfree computation of Jacobians (MatrixFree Methods), two parameters affect the accuracy of the finite difference Jacobian approximation. These are set with the command
MatFDColoringSetParameters(MatFDColoring fdcoloring,PetscReal rerror,PetscReal umin);
The parameter rerror
is the square root of the relative error in the
function evaluations, \(e_{rel}\); the default is the square root of
machine epsilon (about \(10^{8}\) in double precision), which
assumes that the functions are evaluated approximately to floatingpoint
precision accuracy. The second parameter, umin
, is a bit more
involved; its default is \(10e^{6}\) . Column \(i\) of the
Jacobian matrix (denoted by \(F_{:i}\)) is approximated by the
formula
where \(h\) is computed via:
for MATMFFD_DS
or:
for MATMFFD_WP
(default). These parameters may be set from the options
database with
mat_fd_coloring_err <err>
mat_fd_coloring_umin <umin>
mat_fd_type <htype>
Note that MatColoring
type MATCOLORINGSL
, MATCOLORINGLF
, and
MATCOLORINGID
are sequential algorithms. MATCOLORINGJP
and
MATCOLORINGGREEDY
are parallel algorithms, although in practice they
may create more colors than the sequential algorithms. If one computes
the coloring iscoloring
reasonably with a parallel algorithm or by
knowledge of the discretization, the routine MatFDColoringCreate()
is scalable. An example of this for 2D distributed arrays is given below
that uses the utility routine DMCreateColoring()
.
DMCreateColoring(da,IS_COLORING_GHOSTED, &iscoloring);
MatFDColoringCreate(J,iscoloring, &fdcoloring);
MatFDColoringSetFromOptions(fdcoloring);
ISColoringDestroy( &iscoloring);
Note that the routine MatFDColoringCreate()
currently is only
supported for the AIJ and BAIJ matrix formats.
Variational Inequalities#
SNES
can also solve variational inequalities with box constraints.
These are nonlinear algebraic systems with additional inequality
constraints on some or all of the variables:
\(Lu_i \le u_i \le Hu_i\). Some or all of the lower bounds may be
negative infinity (indicated to PETSc with SNES_VI_NINF
) and some or
all of the upper bounds may be infinity (indicated by SNES_VI_INF
).
The command
SNESVISetVariableBounds(SNES,Vec Lu,Vec Hu);
is used to indicate that one is solving a variational inequality. The
option snes_vi_monitor
turns on extra monitoring of the active set
associated with the bounds and snes_vi_type
allows selecting from
several VI solvers, the default is preferred.
Nonlinear Preconditioning#
The mathematical framework of nonlinear preconditioning is explained in detail in [BKST15].
Nonlinear preconditioning in PETSc involves the use of an inner SNES
instance to define the step for an outer SNES
instance. The inner
instance may be extracted using
SNESGetNPC(SNES snes,SNES *npc);
and passed runtime options using the npc_
prefix. Nonlinear
preconditioning comes in two flavors: left and right. The side may be
changed using snes_npc_side
or SNESSetNPCSide()
. Left nonlinear
preconditioning redefines the nonlinear function as the action of the
nonlinear preconditioner \(\mathbf{M}\);
Right nonlinear preconditioning redefines the nonlinear function as the function on the action of the nonlinear preconditioner;
which can be interpreted as putting the preconditioner into “striking distance” of the solution by outer acceleration.
In addition, basic patterns of solver composition are available with the
SNESType
SNESCOMPOSITE
. This allows for two or more SNES
instances to be combined additively or multiplicatively. By command
line, a set of SNES
types may be given by comma separated list
argument to snes_composite_sneses
. There are additive
(SNES_COMPOSITE_ADDITIVE
), additive with optimal damping
(SNES_COMPOSITE_ADDITIVEOPTIMAL
), and multiplicative
(SNES_COMPOSITE_MULTIPLICATIVE
) variants which may be set with
SNESCompositeSetType(SNES,SNESCompositeType);
New subsolvers may be added to the composite solver with
and accessed with
 BS90(1,2)
Peter N. Brown and Youcef Saad. Hybrid Krylov methods for nonlinear systems of equations. SIAM J. Sci. Stat. Comput., 11:450–481, 1990.
 BKST15
Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu. Composing scalable nonlinear algebraic solvers. SIAM Review, 57(4):535–565, 2015. http://www.mcs.anl.gov/papers/P20100112.pdf. URL: http://www.mcs.anl.gov/papers/P20100112.pdf, doi:10.1137/130936725.
 EW96
S. C. Eisenstat and H. F. Walker. Choosing the forcing terms in an inexact Newton method. SIAM J. Scientific Computing, 17:16–32, 1996.
 JP93
Mark T. Jones and Paul E. Plassmann. A parallel graph coloring heuristic. SIAM J. Sci. Comput., 14(3):654–669, 1993.
 MoreSGH84(1,2,3,4)
Jorge J. Moré, Danny C. Sorenson, Burton S. Garbow, and Kenneth E. Hillstrom. The MINPACK project. In Wayne R. Cowell, editor, Sources and Development of Mathematical Software, 88–111. 1984.
 PW98
M. Pernice and H. F. Walker. NITSOL: a Newton iterative solver for nonlinear systems. SIAM J. Sci. Stat. Comput., 19:302–318, 1998.
 DennisJrS83
J. E. Dennis Jr. and Robert B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. PrenticeHall, Inc., Englewood Cliffs, NJ, 1983.