PETSc

1
PETSc
Material adapted from a tutorial by

• Satish Balay
• Bill Gropp
• Lois Curfman McInnes
• Barry Smith
• www-fp.mcs.anl.gov/petsc/docs/tutorials/index.htm
• Mathematics and Computer Science Division
• Argonne National Laboratory
• http//www.mcs.anl.gov/petsc

2
PETSc Philosophy

• Writing hand-parallelized (message-passing or
  distributed shared memory) application codes from
  scratch is extremely difficult and time
  consuming.
• Scalable parallelizing compilers for real
  application codes are very far in the future.
• We can ease the development of parallel
  application codes by developing general-purpose,
  parallel numerical PDE libraries.

3
PETSc Concepts in the Tutorial

• How to specify the mathematics of the problem
• Data objects
• vectors, matrices
• How to solve the problem
• Solvers
• linear, nonlinear, and timestepping (ODE) solvers
• Parallel computing complications
• Parallel data layout
• structured and unstructured meshes

4
Tutorial Outline

• Getting started
• sample results
• programming paradigm
• Data objects
• vectors (e.g., field variables)
• matrices (e.g., sparse Jacobians)
• Viewers
• object information
• visualization
• Solvers
• linear
• nonlinear
• timestepping (and ODEs)

• Data layout and ghost values
• structured mesh problems
• unstructured mesh problems
• partitioning and coloring
• Putting it all together
• a complete example
• Debugging and error handling
• Profiling and performance tuning
• -log_summary
• stages and preloading
• user-defined events
• Extensibility issues

5
Tutorial Approach
From the perspective of an application programmer

• Advanced
• user-defined customization of algorithms and data
  structures
• Developer
• advanced customizations, intended primarily for
  use by library developers

• Beginner
• basic functionality, intended for use by most
  programmers
• Intermediate
• selecting options, performance evaluation and
  tuning

6
Incremental Application Improvement

• Beginner
• Get the application up and walking
• Intermediate
• Experiment with options
• Determine opportunities for improvement
• Advanced
• Extend algorithms and/or data structures as
  needed
• Developer
• Consider interface and efficiency issues for
  integration and interoperability of multiple
  toolkits

7
The PETSc Programming Model

• Goals
• Portable, runs everywhere
• Performance
• Scalable parallelism
• Approach
• Distributed memory, shared-nothing
• Requires only a compiler (single node or
  processor)
• Access to data on remote machines through MPI
• Can still exploit compiler discovered
  parallelism on each node (e.g., SMP)
• Hide within parallel objects the details of the
  communication
• User orchestrates communication at a higher
  abstract level than message passing

tutorial introduction
8
Collectivity

• MPI communicators (MPI_Comm) specify collectivity
  (processors involved in a computation)
• All PETSc creation routines for solver and data
  objects are collective with respect to a
  communicator, e.g., VecCreate(MPI_Comm comm, int
  m, int M, Vec x)
• Some operations are collective, while others are
  not, e.g.,
• collective VecNorm( )
• not collective VecGetLocalSize()
• If a sequence of collective routines is used,
  they must be called in the same order on each
  processor

tutorial introduction
9
PDE Application Codes
tutorial introduction
10
PETSc Numerical Components
tutorial introduction
11
Flow of Control for PDE Solution
Main Routine
Timestepping Solvers (TS)
Nonlinear Solvers (SNES)
Linear Solvers (SLES)
PETSc
PC
KSP
Application Initialization
Function Evaluation
Jacobian Evaluation
Post- Processing
PETSc code
User code
tutorial introduction
12
True Portability

• Tightly coupled systems
• Cray T3D/T3E
• SGI/Origin
• IBM SP
• Convex Exemplar
• Loosely coupled systems, e.g., networks of
  workstations
• Sun OS, Solaris 2.5
• IBM AIX
• DEC Alpha
• HP

• Linux
• Freebsd
• Windows NT/95

tutorial introduction
13
Data Objects

• Vectors (Vec)
• focus field data arising in nonlinear PDEs
• Matrices (Mat)
• focus linear operators arising in nonlinear PDEs
  (i.e., Jacobians)

• Object creation
• Object assembly
• Setting options
• Viewing
• User-defined customizations

beginner
beginner
intermediate
intermediate
advanced
tutorial outline data objects
14
Vectors

• Fundamental objects for storing field solutions,
  right-hand sides, etc.
• VecCreateMPI(...,Vec )
• MPI_Comm - processors that share the vector
• number of elements local to this processor
• total number of elements
• Each process locally owns a subvector of
  contiguously numbered global indices

proc 0
proc 1
proc 2
proc 3
proc 4
data objects vectors
beginner
15
Vectors Example
include "vec.h" int main(int argc,char argv)
Vec x,y,w / vectors /
Vec z / array of
vectors / double norm,v,v1,v2int n
20,ierr Scalar one 1.0,two 2.0,three
3.0,dots3,dot PetscInitialize(argc,argv,(c
har)0,help) ierr OptionsGetInt(PETSC_NULL,"n
",n,PETSC_NULL) CHKERRA(ierr) ierr
VecCreate(PETSC_COMM_WORLD,PETSC_DECIDE,n,x)
CHKERRA(ierr) ierr VecSetFromOptions(x)CHKER
RA(ierr) ierr VecDuplicate(x,y)CHKERRA(ierr
) ierr VecDuplicate(x,w)CHKERRA(ierr)
ierr VecSet(one,x)CHKERRA(ierr) ierr
VecSet(two,y)CHKERRA(ierr) ierr
VecSet(one,z0)CHKERRA(ierr) ierr
VecSet(two,z1)CHKERRA(ierr) ierr
VecSet(three,z2)CHKERRA(ierr) ierr
VecDot(x,x,dot)CHKERRA(ierr) ierr
VecMDot(3,x,z,dots)CHKERRA(ierr)
16
Vectors Example
ierr PetscPrintf(PETSC_COMM_WORLD,"Vector
length d\n",(int)dot) CHKERRA(ierr)
ierr VecScale(two,x)CHKERRA(ierr) ierr
VecNorm(x,NORM_2,norm)CHKERRA(ierr) v
norm-2.0sqrt((double)n) if (v gt -1.e-10 v lt
1.e-10) v 0.0 ierr PetscPrintf(PETSC_COMM_
WORLD,"VecScale g\n",v) CHKERRA(ierr)
ierr VecDestroy(x)CHKERRA(ierr) ierr
VecDestroy(y)CHKERRA(ierr) ierr
VecDestroy(w)CHKERRA(ierr) ierr
VecDestroyVecs(z,3)CHKERRA(ierr)
PetscFinalize() return 0
17
Vector Assembly

• VecSetValues(Vec,)
• number of entries to insert/add
• indices of entries
• values to add
• mode INSERT_VALUES,ADD_VALUES
• VecAssemblyBegin(Vec)
• VecAssemblyEnd(Vec)

data objects vectors
beginner
18
Parallel Matrix and Vector Assembly

• Processors may generate any entries in vectors
  and matrices
• Entries need not be generated on the processor on
  which they ultimately will be stored
• PETSc automatically moves data during the
  assembly process if necessary

data objects vectors and matrices
beginner
19
Selected Vector Operations
20
Sparse Matrices

• Fundamental objects for storing linear operators
  (e.g., Jacobians)
• MatCreateMPIAIJ(,Mat )
• MPI_Comm - processors that share the matrix
• number of local rows and columns
• number of global rows and columns
• optional storage pre-allocation information

data objects matrices
beginner
21
Parallel Matrix Distribution
Each process locally owns a submatrix of
contiguously numbered global rows.
proc 0

• MatGetOwnershipRange(Mat A, int rstart, int
  rend)
• rstart first locally owned row of global
  matrix
• rend -1 last locally owned row of global matrix

data objects matrices
beginner
22
Matrix Assembly

• MatSetValues(Mat,)
• number of rows to insert/add
• indices of rows and columns
• number of columns to insert/add
• values to add
• mode INSERT_VALUES,ADD_VALUES
• MatAssemblyBegin(Mat)
• MatAssemblyEnd(Mat)

data objects matrices
beginner
23
Viewers

• Printing information about solver and data
  objects
• Visualization of field and matrix data
• Binary output of vector and matrix data

beginner
beginner
intermediate
tutorial outline viewers
24
Viewer Concepts

• Information about PETSc objects
• runtime choices for solvers, nonzero info for
  matrices, etc.
• Data for later use in restarts or external tools
• vector fields, matrix contents
• various formats (ASCII, binary)
• Visualization
• simple x-window graphics
• vector fields
• matrix sparsity structure

beginner
viewers
25
Viewing Vector Fields

• VecView(Vec x,Viewer v)
• Default viewers
• ASCII (sequential) VIEWER_STDOUT_SELF
• ASCII (parallel) VIEWER_STDOUT_WORLD
• X-windows VIEWER_DRAW_WORLD
• Default ASCII formats
• VIEWER_FORMAT_ASCII_DEFAULT
• VIEWER_FORMAT_ASCII_MATLAB
• VIEWER_FORMAT_ASCII_COMMON
• VIEWER_FORMAT_ASCII_INFO
• etc.

Solution components, using runtime option
-snes_vecmonitor
velocity v
velocity u
temperature T
vorticity z
viewers
beginner
26
Viewing Matrix Data

• MatView(Mat A, Viewer v)
• Runtime options available after matrix assembly
• -mat_view_info
• info about matrix assembly
• -mat_view_draw
• sparsity structure
• -mat_view
• data in ASCII
• etc.

viewers
beginner
27
Linear Solvers

• Goal Support the solution of linear systems,
• Axb,
• particularly for sparse, parallel problems
  arising
• within PDE-based models
• User provides
• Code to evaluate A, b

solverslinear
beginner
28
Linear PDE Solution
Main Routine
PETSc
Linear Solvers (SLES)
Solve Ax b
PC
KSP
Application Initialization
Evaluation of A and b
Post- Processing
PETSc code
User code
solverslinear
beginner
29
Linear Solvers (SLES)
SLES Scalable Linear Equations Solvers

• Application code interface
• Choosing the solver
• Setting algorithmic options
• Viewing the solver
• Determining and monitoring convergence
• Providing a different preconditioner matrix
• Matrix-free solvers
• User-defined customizations

beginner
beginner
intermediate
intermediate
intermediate
intermediate
advanced
tutorial outline solvers linear
advanced
30
Linear Solvers in PETSc 2.0
Preconditioners (PC)
Krylov Methods (KSP)

• Block Jacobi
• Overlapping Additive Schwarz
• ICC, ILU via BlockSolve95
• ILU(k), LU (sequential only)
• etc.

• Conjugate Gradient
• GMRES
• CG-Squared
• Bi-CG-stab
• Transpose-free QMR
• etc.

solverslinear
beginner
31
Basic Linear Solver Code (C/C)
SLES sles / linear solver
context / Mat A /
matrix / Vec x, b /
solution, RHS vectors / int n, its
/ problem dimension, number of
iterations / MatCreate(MPI_COMM_WORLD,n,n,A)
/ assemble matrix / VecCreate(MPI_COMM_WORLD,n,
x) VecDuplicate(x,b)
/ assemble RHS vector
/ SLESCreate(MPI_COMM_WORLD,sles)
SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTE
RN) SLESSetFromOptions(sles) SLESSolve(sles,b,x,
its)
solverslinear
beginner
32
Basic Linear Solver Code (Fortran)
SLES sles Mat A Vec
x, b integer n, its, ierr call
MatCreate(MPI_COMM_WORLD,n,n,A,ierr) call
VecCreate(MPI_COMM_WORLD,n,x,ierr) call
VecDuplicate(x,b,ierr) call
SLESCreate(MPI_COMM_WORLD,sles,ierr) call
SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTER
N,ierr) call SLESSetFromOptions(sles,ierr) call
SLESSolve(sles,b,x,its,ierr)
C then assemble matrix and right-hand-side
vector
solverslinear
beginner
33
Customization Options

• Procedural Interface
• Provides a great deal of control on a
  usage-by-usage basis inside a single code
• Gives full flexibility inside an application
• Command Line Interface
• Applies same rule to all queries via a database
• Enables the user to have complete control at
  runtime, with no extra coding

solverslinear
intermediate
34
Setting Solver Options within Code

• SLESGetKSP(SLES sles,KSP ksp)
• KSPSetType(KSP ksp,KSPType type)
• KSPSetTolerances(KSP ksp,double rtol,double
  atol,double dtol, int maxits)
• ...
• SLESGetPC(SLES sles,PC pc)
• PCSetType(PC pc,PCType)
• PCASMSetOverlap(PC pc,int overlap)
• ...

solverslinear
intermediate
35
Setting Solver Options at Runtime

• -ksp_type cg,gmres,bcgs,tfqmr,
• -pc_type lu,ilu,jacobi,sor,asm,
• -ksp_max_it ltmax_itersgt
• -ksp_gmres_restart ltrestartgt
• -pc_asm_overlap ltoverlapgt
• -pc_asm_type basic,restrict,interpolate,none
• etc ...

1
2
solverslinear
36
SLES Review of Basic Usage

• SLESCreate( ) - Create SLES context
• SLESSetOperators( ) - Set linear operators
• SLESSetFromOptions( ) - Set runtime solver
  options for SLES, KSP,PC
• SLESSolve( ) - Run linear solver
• SLESView( ) - View solver options
  actually used at runtime (alternative
  -sles_view)
• SLESDestroy( ) - Destroy solver

solverslinear
beginner
37
SLES Review of Selected Preconditioner Options
1
2
And many more options...
solvers linear preconditioners
38
SLES Review of Selected Krylov Method Options
1
2
And many more options...
solvers linear Krylov methods
39
SLES Runtime Script Example
solverslinear
intermediate
40
Viewing SLES Runtime Options
solverslinear
intermediate
41
SLES Example Programs
Location www.mcs.anl.gov/petsc/src/sles/examples
/tutorials/

• ex1.c, ex1f.F - basic uniprocessor codes
• ex2.c, ex2f.F - basic parallel codes
• ex11.c - using complex numbers
• ex4.c - using different linear system
  and
• preconditioner
  matrices
• ex9.c - repeatedly solving different
  linear systems
• ex15.c - setting a user-defined
  preconditioner

1
2
3
And many more examples ...
solverslinear