Introduction to petsc
1 / 67

Introduction to PETSc - PowerPoint PPT Presentation

  • Uploaded on

Introduction to PETSc. Numerical Software Libraries for the Scalable Solution of PDEs. Satish Balay, Bill Gropp, Lois C. McInnes, Barry Smith Mathematics and Computer Science Division Argonne National Laboratory full tutorial as presented at the SIAM Parallel Processing

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Introduction to PETSc' - genera

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Introduction to petsc

Introduction to PETSc

Numerical Software Libraries for

the Scalable Solution of PDEs

Satish Balay, Bill Gropp, Lois C. McInnes, Barry Smith

Mathematics and Computer Science Division

Argonne National Laboratory

full tutorial as presented at

the SIAM Parallel Processing

Conference, March 1999,

available via


  • Writing hand-parallelized 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.

  • Caveats

    • 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 reduce the development time, but it is not a black-box PDE solver nor a silver bullet.

Petsc concepts as illustrated via a complete nonlinear pde example flow in a driven cavity
PETSc ConceptsAs illustrated via a complete nonlinear PDE example: flow in a driven cavity

  • 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

  • Debugging, profiling, extensibility

Tutorial approach









Tutorial Approach

From the perspective of an application programmer:

  • Beginner

    • basic functionality, intended for use by most programmers

  • Intermediate

    • selecting options, performance evaluation and tuning

  • Advanced

    • user-defined customization of algorithms and data structures

  • Developer

    • advanced customizations, intended primarily for use by library developers

Incremental application improvement
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

Component interactions for numerical pdes











Component Interactions for Numerical PDEs



Cfd on an unstructured grid
CFD on an Unstructured Grid

  • 3D incompressible Euler

  • Tetrahedral grid

  • Up to 11 million unknowns

  • Based on a legacy NASA code, FUN3d, developed by W. K. Anderson

  • Fully implicit steady-state

  • Primary PETSc tools: nonlinear solvers (SNES) and vector scatters (VecScatter)

Results courtesy of collaborators Dinesh Kaushik and David Keyes, Old Dominion Univ., partially funded by NSF and ASCI level 2 grant

Scalability for fixed size problem

Missing data point

Scalability for Fixed-Size Problem



Scalability study
Scalability Study

600 MHz T3E, 2.8M vertices

Multiphase flow
Multiphase Flow

  • Oil reservoir simulation: fully implicit, time-dependent

  • First, and only, fully implicit, parallel compositional simulator

  • 3D EOS model (8 DoF per cell)

  • Structured Cartesian mesh

  • Over 4 million cell blocks, 32 million DoF

  • Primary PETSc tools: linear solvers (SLES)

    • restarted GMRES with Block Jacobi preconditioning

    • Point-block ILU(0) on each processor

  • Over 10.6 gigaflops sustained performance on 128 nodes of an IBM SP. 90+ percent parallel efficiency

Results courtesy of collaborators Peng Wang and Jason Abate, Univ. of Texas at Austin, partially funded by DOE ER FE/MICS

Pc and sp comparison
PC and SP Comparison

179,000 unknowns (22,375 cell blocks)

  • PC: Fast ethernet (100 Megabits/second) network of 300 Mhz Pentium PCs with 66 Mhz bus

  • SP: 128 node IBM SP with 160 MHz Power2super processors and 2 memory cards

The petsc programming model
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

Pde application codes

PETSc PDE Application Codes

ODE Integrators


Nonlinear Solvers,

Unconstrained Minimization


Linear Solvers

Preconditioners + Krylov Methods


Matrices, Vectors, Indices



Profiling Interface

Computation and Communication Kernels


PDE Application Codes

Petsc numerical components

Nonlinear Solvers

Time Steppers

Newton-based Methods





Pseudo Time



Line Search

Trust Region

Krylov Subspace Methods


















(Sequential only)




Sparse Row


Blocked Compressed

Sparse Row








Index Sets


Block Indices



PETSc Numerical Components

Mesh definitions for our purposes
Mesh Definitions: For Our Purposes

  • Structured

    • Determine neighbor relationships purely from logical I, J, K coordinates

    • PETSc support via DA

  • Semi-Structured

    • No direct PETSc support; could work with tools such as SAMRAI, Overture, Kelp

  • Unstructured

    • Do not explicitly use logical I, J, K coordinates

    • PETSc support via lower-level VecScatter utilities

One is always free to manage the mesh data as if it is


A freely available and supported research code
A Freely Available and Supported Research Code

  • Available via

  • Usable in C, C++, and Fortran77/90 (with minor limitations in Fortran 77/90 due to their syntax)

  • Users manual in Postscript and HTML formats

  • Hyperlinked manual pages for all routines

  • Many tutorial-style examples

  • Support via email: [email protected]

True portability
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

Petsc history
PETSc History

  • September 1991: begun by Bill and Barry

  • May 1992: first release of PETSc 1.0

  • January 1993: joined by Lois

  • September 1994: begun PETSc 2.0

  • June 1995: first release of version 2.0

  • October 1995: joined by Satish

  • Now: over 4,000 downloads of version 2.0

  • Department of Energy, MICS Program

  • National Science Foundation, Multidisciplinary Challenge Program, CISE

PETSc Funding and Support

A complete example driven cavity model

Velocity-vorticity formulation, with flow driven by lid and/or bouyancy

Finite difference discretization with 4 DoF per mesh point

Primary PETSc tools: SNES, DA

Solution Components

velocity: v

velocity: u


temperature: T

vorticity: z

A Complete Example:Driven Cavity Model

Example code: petsc/src/snes/examples/tutorials/ex8.c

Application code author: D. E. Keyes

Driven cavity solution approach
Driven Cavity Solution Approach and/or bouyancy


Main Routine


Nonlinear Solvers (SNES)


F(u) = 0

Linear Solvers (SLES)














User code

PETSc code

Driven cavity program
Driven Cavity Program and/or bouyancy

  • PartA: Data objects (Vec and Mat), Solvers (SNES)

  • Part B: Parallel data layout(DA)

  • Part C: Nonlinear function evaluation

    • ghost point updates

    • local function computation

  • Part D: Jacobian evaluation

    • default colored finite differencing approximation

  • Experimentation

Collectivity and/or bouyancy

  • 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

Vectors and/or bouyancy

proc 0

  • 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 1

proc 2

proc 3

proc 4

Sparse matrices

proc 0 and/or bouyancy

proc 1

proc 2

proc 3

proc 4

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

  • Each process locally owns a submatrix of contiguously numbered global rows.

Parallel matrix and vector assembly
Parallel Matrix and Vector Assembly and/or bouyancy

  • Processes may generate any entries in vectors and matrices

  • Entries need not be generated by the process on which they ultimately will be stored

  • PETSc automatically moves data during assembly if necessary

  • Vector example:

  • VecSetValues(Vec,…)

    • number of entries to insert/add

    • indices of entries

    • values to add


  • VecAssemblyBegin(Vec)

  • VecAssemblyEnd(Vec)

Matrix assembly
Matrix Assembly and/or bouyancy

  • MatSetValues(Mat,…)

    • number of rows to insert/add

    • indices of rows and columns

    • number of columns to insert/add

    • values to add


  • MatAssemblyBegin(Mat)

  • MatAssemblyEnd(Mat)

Blocked sparse matrices
Blocked Sparse Matrices and/or bouyancy

  • For multi-component PDEs

  • MatCreateMPIBAIJ(…,Mat *)

    • MPI_Comm - processes that share the matrix

    • block size

    • number of local rows and columns

    • number of global rows and columns

    • optional storage pre-allocation information

  • 3D compressible Euler code

  • Block size 5

  • IBM Power2

Solvers usage concepts

Linear ( and/or bouyancy SLES)

Nonlinear (SNES)

Timestepping (TS)

Context variables

Solver options

Callback routines


important concepts

Solvers: Usage Concepts

Solver Classes

Usage Concepts

Nonlinear solvers snes
Nonlinear Solvers (SNES) and/or bouyancy

Goal: For problems arising from PDEs, support the general solution of F(u) = 0

  • Newton-based methods, including

    • Line search strategies

    • Trust region approaches

    • Pseudo-transient continuation

    • Matrix-free variants

  • User provides:

    • Code to evaluate F(u)

    • Code to evaluate Jacobian of F(u) (optional)

      • or use sparse finite difference approximation

      • or use automatic differentiation (coming soon!)

  • User can customize all phases of solution

Context variables
Context Variables and/or bouyancy

  • Are the key to solver organization

  • Contain the complete state of an algorithm, including

    • parameters (e.g., convergence tolerance)

    • functions that run the algorithm (e.g., convergence monitoring routine)

    • information about the current state (e.g., iteration number)

Creating the snes context
Creating the SNES Context and/or bouyancy

  • C/C++ version


  • Fortran version


  • Provides an identical user interface for all linear solvers

    • uniprocessor and parallel

    • real and complex numbers

Basic nonlinear solver code c c
Basic Nonlinear Solver Code (C/C++) and/or bouyancy

SNES snes; /* nonlinear solver context */

Mat J; /* Jacobian matrix */

Vec x, F; /* solution, residual vectors */

int n, its; /* problem dimension, number of iterations */

ApplicationCtx usercontext; /* user-defined application context */










Basic nonlinear solver code fortran
Basic Nonlinear Solver Code (Fortran) and/or bouyancy

SNES snes

Mat J

Vec x, F

int n, its


call MatCreate(MPI_COMM_WORLD,n,n,J,ierr)

call VecCreate(MPI_COMM_WORLD,n,x,ierr)

call VecDuplicate(x,F,ierr)


call SNESSetFunction(snes,F,EvaluateFunction,PETSC_NULL,ierr)

call SNESSetJacobian(snes,J,EvaluateJacobian,PETSC_NULL,ierr)

call SNESSetFromOptions(snes,ierr)

call SNESSolve(snes,x,its,ierr)

Solvers based on callbacks

important concept and/or bouyancy

Solvers Based on Callbacks

  • User provides routines to perform actions that the library requires. For example,

    • SNESSetFunction(SNES,...)

      • uservector - vector to store function values

      • userfunction - name of the user’s function

      • usercontext - pointer to private data for the user’s function

  • Now, whenever the library needs to evaluate the user’s nonlinear function, the solver may call the application code directly with its own local state.

  • usercontext: serves as an application context object. Data are handled through such opaque objects; the library never sees irrelevant application data

Sample application context driven cavity problem
Sample Application Context: and/or bouyancy Driven Cavity Problem

typedef struct {

/* - - - - - - - - - - - - - - - - basic application data - - - - - - - - - - - - - - - - - */

double lid_velocity, prandtl, grashof; /* problem parameters */

int mx, my; /* discretization parameters */

int mc; /* number of DoF per node */

int draw_contours; /* flag - drawing contours */

/* - - - - - - - - - - - - - - - - - - - - parallel data - - - - - - - - - - - - - - - - - - - - - */

MPI_Comm comm; /* communicator */

DA da; /* distributed array */

Vec localF, localX; /* local ghosted vectors */

} AppCtx;

Sample function evaluation code driven cavity problem
Sample Function Evaluation Code: and/or bouyancy Driven Cavity Problem

UserComputeFunction(SNES snes, Vec X, Vec F, void *ptr)


AppCtx *user = (AppCtx *) ptr; /* user-defined application context */

int istart, iend, jstart, jend; /* local starting and ending grid points */

Scalar *f; /* local vector data */


/* Communicate nonlocal ghost point data */

VecGetArray( F, &f );

/* Compute local function components; insert into f[ ] */

VecRestoreArray( F, &f );


return 0;


Sample local computational loops driven cavity problem
Sample Local Computational Loops: and/or bouyancy Driven Cavity Problem

#define U(i) 4*(i)

#define V(i) 4*(i)+1

#define Omega(i) 4*(i)+2

#define Temp(i) 4*(i)+3


for ( j = jstart; j<jend; j++ ) {

row = (j - gys) * gxm + istart - gxs - 1;

for ( i = istart; i<iend; i++ ) {

row++; u = x[U(row)];

uxx = (two * u - x[ U (row-1)] - x [U (row+1) ] ) * hydhx;

uyy = (two * u - x[ U (row-gxm)] - x [ U (row+gxm) ] ) * hxdhy;

f [U(row)] = uxx + uyy - p5 * (x [ (Omega (row+gxm)] - x [Omega (row-gxm)] ) * hx;

} }


Customization options
Customization Options and/or bouyancy

  • 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

Setting solver options within code
Setting Solver Options within Code and/or bouyancy

  • SNESSetTolerances(SNES snes, double atol, double rtol, double stol, int maxits, int maxf)

  • SNESSetType(SNES snes,SNESType type)

  • etc.

Uniform access to all linear and nonlinear solvers
Uniform access to all linear and nonlinear solvers and/or bouyancy

  • -ksp_type [cg,gmres,bcgs,tfqmr,…]

  • -pc_type [lu,ilu,jacobi,sor,asm,…]

  • -snes_type [ls,tr,…]

  • -snes_line_search <line search method>

  • -sles_ls <parameters>

  • -snes_rtol <relative_tol>

  • etc...



Runtime script example
Runtime Script Example and/or bouyancy

Customization via callbacks setting a user defined line search routine
Customization via Callbacks: and/or bouyancy Setting a user-defined line search routine

SNESSetLineSearch(SNES snes,int(*ls)(…),void *lsctx)

where available line search routines ls( ) include:

  • SNESCubicLineSearch( ) - cubic line search

  • SNESQuadraticLineSearch( ) - quadratic line search

  • SNESNoLineSearch( ) - full Newton step

  • SNESNoLineSearchNoNorms( ) - full Newton step but calculates no norms (faster in parallel, useful when using a fixed number of Newton iterations instead of usual convergence testing)

  • YourOwnFavoriteLineSearchRoutine( ) - whatever you like



Snes review of basic usage
SNES: Review of Basic Usage and/or bouyancy

  • SNESCreate( ) - Create SNES context

  • SNESSetFunction( ) - Set function eval. routine

  • SNESSetJacobian( ) - Set Jacobian eval. routine

  • SNESSetFromOptions( ) - Set runtime solver options for [SNES,SLES,KSP,PC]

  • SNESSolve( ) - Run nonlinear solver

  • SNESView( ) - View solver options actually used at runtime (alternative: -snes_view)

  • SNESDestroy( ) - Destroy solver

Snes review of selected options
SNES: Review of Selected Options and/or bouyancy



And many more options...

Snes example programs
SNES: Example Programs and/or bouyancy


  • ex1.c, ex1f.F - basic uniprocessor codes

  • ex4.c, ex4f.F - uniprocessor nonlinear PDE (1 DoF per node)

  • ex5.c, ex5f.F - parallel nonlinear PDE (1 DoF per node)

  • ex8.c - parallel nonlinear PDE (multiple DoF per node, driven cavity problem)



  • And many more examples ...

Data layout and ghost values usage concepts

Structured and/or bouyancy

DA objects


VecScatter objects

Geometric data

Data structure creation

Ghost point updates

Local numerical computation

important concepts

Data Layout and Ghost Values : Usage Concepts

Managing field data layout and required ghost values is the key to high performance of most PDE-based parallel programs.

Mesh Types

Usage Concepts

Ghost values

Local node and/or bouyancy

Ghost node

Ghost Values

Ghost values: To evaluate a local function f(x) , each process requires its local portion of the vector x as well as its ghost values -- or bordering portions of x that are owned by neighboring processes.

Communication and physical discretization
Communication and Physical Discretization and/or bouyancy







Data Structure


Ghost Point

Data Structures

Ghost Point






Loops over



DACreate( )

DAGlobalToLocal( )

structured meshes







Loops over


VecScatter( )

VecScatterCreate( )

unstructured meshes


Global and local representations

5 and/or bouyancy







Local node

Ghost node

Global and Local Representations

Global: each process stores a unique local set of vertices (and each vertex is owned by exactly one process)

Local: each process stores a unique local set of vertices as well as ghost nodes from neighboring processes

Distributed arrays
Distributed Arrays and/or bouyancy

Data layout and ghost values

Proc 10

Proc 10

Proc 0

Proc 1

Proc 0

Proc 1





Logically regular meshes
Logically Regular Meshes and/or bouyancy

  • DA - Distributed Array: object containing info about vector layout across processes and communication of ghost values

  • Form a DA (1, 2, or 3-dim)

    • DACreate2d(….,DA *)


      • number of grid points in x- and y-directions

      • processors in x- and y-directions

      • degrees of freedom per node

      • stencil width

      • ...

  • Update ghostpoints

    • DAGlobalToLocalBegin(DA,…)

    • DAGlobalToLocalEnd(DA,…)

Vectors and das
Vectors and DAs and/or bouyancy

  • The DA object contains information about the data layout and ghost values, but not the actual field data, which is contained in PETSc vectors

  • Global vector: parallel

    • each process stores a unique local portion

    • DACreateGlobalVector(DA da,Vec *gvec);

  • Local work vector: sequential

    • each processor stores its local portion plus ghost values

    • DACreateLocalVector(DA da,Vec *lvec);

    • uses “natural” local numbering of indices (0,1,…nlocal-1)

Updating the local representation
Updating the Local Representation and/or bouyancy

Two-step process that enables overlapping computation and communication

  • DAGlobalToLocalBegin(DA,

    • Vec global_vec,


    • Vec local_vec);

  • DAGlobalToLocal End(DA,…)

Profiling and/or bouyancy

  • Integrated monitoring of

    • time

    • floating-point performance

    • memory usage

    • communication

  • All PETSc events are logged if compiled with -DPETSC_LOG (default); can also profile application code segments

  • Print summary data with option: -log_summary

  • See supplementary handout with summary data

Driven cavity running the program 1
Driven Cavity: and/or bouyancy Running the program (1)

Matrix-free Jacobian approximation with no preconditioning (via -snes_mf) … does not use explicit Jacobian evaluation

  • 1 processor: (thermally-driven flow)

    • mpirun -np 1 ex8 -snes_mf -snes_monitor -grashof 1000.0 -lidvelocity 0.0

  • 2 processors, view DA (and pausing for mouse input):

    • mpirun -np 2 ex8 -snes_mf -snes_monitor -da_view_draw -draw_pause -1

  • View contour plots of converging iterates

    • mpirun ex8 -snes_mf -snes_monitor -snes_vecmonitor

Driven cavity running the program 2
Driven Cavity: and/or bouyancy Running the program (2)

  • Use MatFDColoring for sparse finite difference Jacobian approximation; view SNES options used at runtime

    • mpirun ex8 -snes_view -mat_view_info

  • Set trust region Newton method instead of default line search

    • mpirun ex8 -snes_type tr -snes_view -snes_monitor

  • Set transpose-free QMR as Krylov method; set relative KSP convergence tolerance to be .01

    • mpirun ex8 -ksp_type tfqmr -ksp_rtol .01 -snes_monitor

Sample error traceback
Sample Error Traceback and/or bouyancy

Breakdown in ILU factorization due to a zero pivot

Summary and/or bouyancy

  • Using callbacks to set up the problems for nonlinear solvers

  • Managing data layout and ghost point communication with DAs (for structured meshes) and VecScatters (for unstructured meshes)

  • Evaluating parallel functions and Jacobians

  • Consistent profiling and error handling

Extensibility issues
Extensibility Issues and/or bouyancy

  • Most PETSc objects are designed to allow one to “drop in” a new implementation with a new set of data structures (similar to implementing a new class in C++).

  • Heavily commented example codes include

    • Krylov methods: petsc/src/sles/ksp/impls/cg

    • preconditioners: petsc/src/sles/pc/impls/jacobi

  • Feel free to discuss more details with us in person.

Caveats revisited
Caveats Revisited and/or bouyancy

  • 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.

  • Users are invited to interact directly with us regarding correctness or performance issues by writing to [email protected]

Where is petsc headed next
Where is PETSc headed next? and/or bouyancy

  • As parallel programming models continue to evolve, so will PETSc

  • New algorithms

  • Extended tools for parallel mesh/vector manipulation

  • Numerical software interoperability issues, e.g., in collaboration with

    • Argonne colleagues in the ALICE project


    • DOE colleagues in the Common Component Architecture (CCA) Forum


References and/or bouyancy

  • PETSc Users Manual, by Balay, Gropp, McInnes, and Smith

  • Using MPI, by Gropp, Lusk, and Skjellum

  • MPI Homepage:

  • Domain Decomposition, by Smith, Bjorstad, and Gropp

Petsc at nersc
PETSc at NERSC and/or bouyancy

  • Installed on CRAY T3E-900:

  • Info available via:

  • Access via modules package

    • module load petsc

  • Sample makefile

    • $PETSC_DIR/src/sles/examples/tutorials/makefile

  • Sample example program

    • $PETSC_DIR/src/sles/examples/tutorials/ex2.c

      • this example is discussed in detail in Part I of the PETSc Users Manual, available via