PETSc Tutorial Satish Balay, Kris Buschelman, Bill Gropp,
Dinesh Kaushik, Matt Knepley,
Lois Curfman McInnes, Barry Smith, Hong Zhang
Mathematics and Computer Science Division
Argonne National Laboratory
http://www.mcs.anl.gov/petsc
Intended for use with version 2.1.0 of PETScNumerical Software Libraries for
the Scalable Solution of PDEs
PETSc Tutorial Satish Balay, Kris Buschelman, Bill Gropp,
Dinesh Kaushik, Matt Knepley,
Lois Curfman McInnes, Barry Smith, Hong Zhang
Mathematics and Computer Science Division
Argonne National Laboratory
http://www.mcs.anl.gov/petsc
Intended for use with version 2.1.0 of PETSc
Numerical Software Libraries for
the Scalable Solution of PDEs
The Role of PETSc
Developing parallel, non-trivial PDE solvers that deliver high performance is still difficult and requires months (or even years) of concentrated effort.
PETSc is a toolkit that can ease these difficulties and reduce the development time, but it is not a black-box PDE solver nor a silver bullet.
What is PETSc?
A freely available and supported research code
Available via http://www.mcs.anl.gov/petsc
Free for everyone, including industrial users
Hyperlinked documentation and manual pages for all routines
Many tutorial-style examples
Support via email: petsc-maint@mcs.anl.gov
Usable from Fortran 77/90, C, and C++
Portable to any parallel system supporting MPI, including
Tightly coupled systems
Cray T3E, SGI Origin, IBM SP, HP 9000, Sun Enterprise
Loosely coupled systems, e.g., networks of workstations
Compaq, HP, IBM, SGI, Sun
PCs running Linux or Windows
PETSc history
Begun in September 1991
Now: over 8,500 downloads since 1995 (versions 2.0 and 2.1)
PETSc funding and support
Department of Energy: MICS Program, DOE2000, SciDAC
National Science Foundation, Multidisciplinary Challenge Program, CISE
PETSc Concepts
How to specify the mathematics of the problem
Data objects
vectors, matrices
How to solve the problem
Solvers
linear, nonlinear, and time stepping (ODE) solvers
Parallel computing complications
Parallel data layout
structured and unstructured meshes
Structure of PETSc
Computation and Communication Kernels
MPI, MPI-IO, BLAS, LAPACK Profiling Interface PETSc PDE Application Codes Object-Oriented
Matrices, Vectors, Indices Grid
Management Linear Solvers
Preconditioners + Krylov Methods Nonlinear Solvers,
Unconstrained Minimization ODE Integrators Visualization Interface
PETSc Numerical Components
Compressed
Sparse Row
(AIJ) Blocked Compressed
Sparse Row
(BAIJ) Block
Diagonal
(BDIAG) Dense Other Indices Block Indices Stride Other Index Sets Vectors Line Search Trust Region Newton-based Methods Other Nonlinear Solvers Additive
Schwartz Block
Jacobi Jacobi ILU ICC LU
(Sequential only) Others Preconditioners Euler Backward
Euler Pseudo Time
Stepping Other Time Steppers GMRES CG CGS Bi-CG-STAB TFQMR Richardson Chebychev Other Krylov Subspace Methods Matrices Distributed Arrays Matrix-free
What is not in PETSc?
But PETSc does interface to external software that provides some of this functionality. Discretizations
Unstructured mesh generation and refinement tools
Load balancing tools
Sophisticated visualization capabilities
Flow of Control for PDE Solution
PETSc code User code Application
Initialization Function
Evaluation Jacobian
Evaluation Post-
Processing PC KSP PETSc Main Routine Linear Solvers (SLES) Nonlinear Solvers (SNES) Timestepping Solvers (TS)
Levels of Abstraction in Mathematical Software
PETSc
emphasis Application-specific interface
Programmer manipulates objects associated with the application
High-level mathematics interface
Programmer manipulates mathematical objects, such as PDEs and boundary conditions
Algorithmic and discrete mathematics interface
Programmer manipulates mathematical objects (sparse matrices, nonlinear equations), algorithmic objects (solvers) and discrete geometry (meshes)
Low-level computational kernels
e.g., BLAS-type operations
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
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.
Hello World
#include “petsc.h” int main( int argc, char *argv[] ) { PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);
PetscPrintf(PETSC_COMM_WORLD,“Hello World\n”);
PetscFinalize(); return 0; }
Data Objects
Vectors (Vec)
focus: field data arising in nonlinear PDEs
Matrices (Mat)
focus: linear operators arising in nonlinear PDEs (i.e., Jacobians) tutorial outline: data objects beginner beginner intermediate advanced intermediate Object creation
Object assembly
Setting options
Viewing
User-defined customizations
Vectors
data objects: vectors beginner proc 3 proc 2 proc 0 proc 4 proc 1 What are PETSc vectors?
Fundamental objects for storing field solutions, right-hand sides, etc.
Each process locally owns a subvector of contiguously numbered global indices
Create vectors via
VecCreate(...,Vec *)
MPI_Comm - processors that share the vector
number of elements local to this processor
or total number of elements
VecSetType(Vec,VecType)
Where VecType is
VEC_SEQ, VEC_MPI, or VEC_SHARED
Vector Assembly
data objects: vectors beginner VecSetValues(Vec,…)
number of entries to insert/add
indices of entries
values to add
mode: [INSERT_VALUES,ADD_VALUES]
VecAssemblyBegin(Vec)
VecAssemblyEnd(Vec)
Parallel Matrix and Vector Assembly
data objects: vectors and matrices beginner 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
Matrices
data objects: matrices beginner What are PETSc matrices?
Fundamental objects for storing linear operators (e.g., Jacobians)
Create matrices via
MatCreate(…,Mat *)
MPI_Comm - processors that share the matrix
number of local/global rows and columns
MatSetType(Mat,MatType)
where MatType is one of
default sparse AIJ: MPIAIJ, SEQAIJ
block sparse AIJ (for multi-component PDEs): MPIAIJ, SEQAIJ
symmetric block sparse AIJ: MPISBAIJ, SAEQSBAIJ
block diagonal: MPIBDIAG, SEQBDIAG
dense: MPIDENSE, SEQDENSE
matrix-free
etc.
Matrix Assembly
data objects: matrices beginner 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)
Parallel Matrix Distribution
Each process locally owns a submatrix of contiguously numbered global rows. proc 0 } proc 3: locally owned rows proc 3 proc 2 proc 1 proc 4 data objects: matrices beginner MatGetOwnershipRange(Mat A, int *rstart, int *rend)
rstart: first locally owned row of global matrix
rend -1: last locally owned row of global matrix
Viewers
tutorial outline: viewers beginner beginner intermediate Printing information about solver and data objects
Visualization of field and matrix data
Binary output of vector and matrix data
Comments