# Presentation: 895994

## PETSc Tutorial

- the Scalable Solution of PDEs

- Kris Buschelman
- Barry Smith
- Mathematics and Computer Science Division
- Argonne National Laboratory
- _[http://www.mcs.anl.gov/petsc](http://www.mcs.anl.gov/petsc)_
*Intended for use with version 2.1.3 of PETSc*

## The PETSc Team

- Dinesh
- Kaushik

- Lois
- Curfman
- McInnes

- Kris
- Buschelman

- Satish
- Balay

- Bill
- Gropp

- Barry
- Smith

- Hong
- Zhang

- Matt
- Knepley

## Tutorial Objectives

- Introduce the Portable, Extensible Toolkit for Scientific Computation (PETSc)
- Demonstrate how to write a complete parallel implicit PDE solver using PETSc
- Introduce PETSc interfaces to other software packages
- Explain how to learn more about PETSc

## 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 (Open Source, non-GPL)
- 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
- 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
- Data objects
- vectors, matrices
- How to solve
- Solver objects
- linear, nonlinear, and time stepping (ODE) solvers
- Parallel computing complications
- Parallel data layout
- structured and unstructured meshes

## Tutorial Topics (Long Course)

**Getting started**
- motivating examples
- programming paradigm
**Data objects**
- vectors (e.g., field
- variables)
- matrices (e.g., sparse Jacobians)
**Viewers**
- object information
- visualization
**Solvers**
- Linear
**Profiling and     performance tuning**

**Solvers (cont.)**
- nonlinear
- timestepping (and ODEs)
**Data layout and ghost values**
- structured and unstructured mesh problems
**Putting it all together**
- a complete example
**Debugging and error handling**
**New features**
**Using PETSc with other software packages**

## Slide 8

## Using PETSc with Other Packages

- Linear solvers
- AMG  _[http://www.](http://www.mgnet.org/mgnet-codes-gmd.html)__[mgnet](http://www.mgnet.org/mgnet-codes-gmd.html)__[.org/](http://www.mgnet.org/mgnet-codes-gmd.html)__[mgnet](http://www.mgnet.org/mgnet-codes-gmd.html)__[-codes-](http://www.mgnet.org/mgnet-codes-gmd.html)__[gmd](http://www.mgnet.org/mgnet-codes-gmd.html)__[.html](http://www.mgnet.org/mgnet-codes-gmd.html)_
- BlockSolve95  _[http://www.mcs.anl.gov/BlockSolve95](http://www.mcs.anl.gov/BlockSolve95)_
- ILUTP  _[http://www.](http://www.cs.umn.edu/~saad/)__[cs](http://www.cs.umn.edu/~saad/)__[.](http://www.cs.umn.edu/~saad/)__[umn](http://www.cs.umn.edu/~saad/)__[.](http://www.cs.umn.edu/~saad/)__[edu](http://www.cs.umn.edu/~saad/)__[/~](http://www.cs.umn.edu/~saad/)__[saad](http://www.cs.umn.edu/~saad/)__[/](http://www.cs.umn.edu/~saad/)_
- LUSOL  _[http://www.](http://www.sbsi-sol-optimize.com/)__[sbsi](http://www.sbsi-sol-optimize.com/)__[-sol-optimize.com](http://www.sbsi-sol-optimize.com/)_
- SPAI  _[http://www.](http://www.sam.math.ethz.ch/~grote/spai)__[sam](http://www.sam.math.ethz.ch/~grote/spai)__[.math.](http://www.sam.math.ethz.ch/~grote/spai)__[ethz](http://www.sam.math.ethz.ch/~grote/spai)__[.](http://www.sam.math.ethz.ch/~grote/spai)__[ch](http://www.sam.math.ethz.ch/~grote/spai)__[/~](http://www.sam.math.ethz.ch/~grote/spai)__[grote](http://www.sam.math.ethz.ch/~grote/spai)__[/](http://www.sam.math.ethz.ch/~grote/spai)__[spai](http://www.sam.math.ethz.ch/~grote/spai)_
- SuperLU  _[http://www.](http://www.nersc.gov/~xiaoye/SuperLU)__[nersc](http://www.nersc.gov/~xiaoye/SuperLU)__[.gov/~](http://www.nersc.gov/~xiaoye/SuperLU)__[xiaoye](http://www.nersc.gov/~xiaoye/SuperLU)__[/](http://www.nersc.gov/~xiaoye/SuperLU)__[SuperLU](http://www.nersc.gov/~xiaoye/SuperLU)_
- Spooles http://www.netlib.org/linalg/spooles/spooles.2.2.html
- DSCPack _http://www.cse.psu.edu/~raghavan/_
- Hypre _http://www.llnl.gov/CASC/hypre_
- Optimization software
- TAO  _[http://www.mcs.anl.gov/](http://www.mcs.anl.gov/tao)__[tao](http://www.mcs.anl.gov/tao)_
- Veltisto  _[http://www.cs.nyu.edu/~biros/](http://www.cs.nyu.edu/~biros/veltisto)__[veltisto](http://www.cs.nyu.edu/~biros/veltisto)_

- Mesh and discretization tools
- Overture  _[http://www.](http://www.llnl.gov/CASC/Overture)__[llnl](http://www.llnl.gov/CASC/Overture)__[.gov/CASC/Overture](http://www.llnl.gov/CASC/Overture)_
- SAMRAI  _[http://www.](http://www.llnl.gov/CASC/SAMRAI)__[llnl](http://www.llnl.gov/CASC/SAMRAI)__[.gov/CASC/SAMRAI](http://www.llnl.gov/CASC/SAMRAI)_
- SUMAA3d  _[http://www.mcs.anl.gov/sumaa3d](http://www.mcs.anl.gov/sumaa3d)_
- GlobalArrays
- SciRun
- ODE solvers
- PVODE  _http://www.llnl.gov/CASC/SUNDIALS_
- Others
- Matlab  _[http://www.](http://www.mathworks.com/)__[mathworks](http://www.mathworks.com/)__[.com](http://www.mathworks.com/)_
- ParMETIS http://www.cs.umn.edu/~karypis/metis/parmetis
- ADIFOR/ADIC

## CFD on an Unstructured Mesh

- 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 Dinesh Kaushik and David Keyes, Old Dominion Univ., partially funded by NSF and ASCI level 2 grant*

## Fixed-size Parallel Scaling Results (GFlop/s)

- Dimension=11,047,096

## Inside the Parallel Scaling Results on ASCI Red

- ONERA M6 wing test case, tetrahedral grid of 2.8 million vertices (about 11 million unknowns)
- on up to 3072 ASCI Red nodes (each with dual Pentium Pro 333 MHz processors)

## Multiphase Flow

- Oil reservoir simulation: fully implicit, time-dependent
- First 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 process
- 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*

## Speedup Comparison

## Structures Simulations

- ALE3D (LLNL structures code) test problems
- Simulation with over 16 million degrees of freedom
- Run on NERSC 512 processor T3E and LLNL ASCI Blue Pacific
- Primary PETSc tools: multigrid linear solvers (SLES)

*Results courtesy of Mark Adams (Univ. of California, Berkeley)*

## ALE3D Test Problem Performance

- NERSC Cray T3E Scaled Performance15,000 DoF per processor

## Tutorial Approach

- 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

- From the perspective of an application programmer:

## 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
- Full tutorials available at http://www.mcs.anl.gov/petsc/docs/tutorials

## Structure of PETSc

## PETSc Numerical Components

- 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?

- Discretizations
- Unstructured mesh generation and refinement tools
- Load balancing tools
- Sophisticated visualization capabilities

- But PETSc does interface to external software that provides some of this functionality.

## 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)

## Flow of Control for PDE Solution

- KSP

- PETSc

- Linear Solvers (SLES)

- Application
- Initialization

- Function
- Evaluation

- Jacobian
- Evaluation

- Post-
- Processing

- Other Tools

- SAMRAI

- Overture

- Nonlinear Solvers (SNES)

- Main Routine

- SPAI

- ILUDTP

- PC

- PVODE

- Timestepping Solvers (TS)

## Levels of Abstraction in Mathematical Software

- 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

## Focus On Implicit Methods

**Implicit**: Some or all variables are updated in a global linear or nonlinear solve
- Explicit and semi-explicit are easier cases
- No direct PETSc support for
- ADI-type schemes
- spectral methods
- particle-type methods

## Numerical Methods Paradigm

- Encapsulate the latest numerical algorithms in a consistent, application-friendly manner
- Use mathematical and algorithmic objects, not low-level programming language objects
- Application code focuses on mathematics of the global problem, not parallel programming details

## PETSc Programming Aids

- Correctness Debugging
- Automatic generation of tracebacks
- Detecting memory corruption and leaks
- Optional user-defined error handlers
- Performance Debugging
- Integrated profiling using -log_summary
- Profiling by stages of an application
- User-defined events

## 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 (processes involved in a computation)
- All PETSc routines for creating solver and data objects are collective with respect to a communicator, e.g.,
- VecCreate(MPI_Comm comm, int m, int M, Vec *x)
- Use **PETSC_COMM_WORLD** for all processes (like MPI_COMM_WORLD, but allows the same code to work when PETSc is started with a smaller set of processes)
- 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 by each process.

## 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;}

## Hello World (Fortran)

-     program main
-     integer ierr, rank
- #include "include/finclude/petsc.h"
-     call PetscInitialize( PETSC_NULL_CHARACTER, ierr )
-     call MPI_Comm_rank( PETSC_COMM_WORLD, rank, ierr )
-     if (rank .eq. 0) then
-         print *, ‘Hello World’
-     endif
-     call PetscFinalize(ierr)
-     end

## Fancier Hello World

- #include “petsc.h”int main( int argc, char *argv[] ){  int rank;
-   PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);
-   MPI_Comm_rank(PETSC_COMM_WORLD,&rank );  PetscSynchronizedPrintf(PETSC_COMM_WORLD,
-                                “Hello World from %d\n”,rank);
-   PetscSynchronizedFlush(PETSC_COMM_WORLD);  PetscFinalize();  return 0;}

## Data Objects

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

- 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

## Vectors

- 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(MPI_Comm,Vec *)
- MPI_Comm - processes that share the vector
- VecSetSizes( Vec, int, int )
- number of elements local to this process
- or total number of elements
- VecSetType(Vec,VecType)
- Where VecType is
- VEC_SEQ, VEC_MPI, or VEC_SHARED
- VecSetFromOptions(Vec) lets you set the type at *runtime*

- data objects: vectors

- beginner

- proc 3

- proc 2

- proc 0

- proc 4

- proc 1

## Creating a vector

- Vec x;int n;..._[PetscInitialize](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscInitialize.html#PetscInitialize)_(&argc,&argv,(char*)0,help);_[PetscOptionsGetInt](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscOptionsGetInt.html#PetscOptionsGetInt)_(PETSC_NULL,"-n",&n,PETSC_NULL);... _[VecCreate](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecCreate.html#VecCreate)_(PETSC_COMM_WORLD,&x);_[VecSetSizes](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetSizes.html#VecSetSizes)_(x,PETSC_DECIDE,n);_[VecSetType](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetType.html)_(x,VEC_MPI);_[VecSetFromOptions](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetFromOptions.html#VecSetFromOptions)_(x);

- Global size

- PETSc determines local size

- Use PETSc to get value from command line

## How Can We Use a PETSc Vector

- In order to support the distributed memory “shared nothing” model, as well as single processors and shared memory systems, a PETSc vector is a “handle” to the real vector
- Allows the vector to be distributed across many processes
- To access the *elements* of the vector, we cannot simply do  		for (i=0; i<n; i++) v[i] = i;
- We do not want to *require *that the programmer work only with the “local” part of the vector; we want to permit operations, such as setting an element of a vector, to be performed by any process.
- Recall how data is stored in the distributed memory programming model...

## Distributed Memory Model

- Integer A(10)...**print *, A**

- A(10)

- A(10)

- Integer A(10)do i=1,10  A(i) = ienddo...

- Process 0

- Process 1

*Different Variables!*

- Address Space

- This A is completely different from this one

## Vector Assembly

- A three step process
- Each process tells PETSc what values to set or add to a vector component.  Once *all* values provided,
- Begin communication between processes to ensure that values end up where needed
- (allow other operations, such as some computation, to proceed)
- Complete the communication
- 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

## Parallel Matrix and Vector Assembly

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

- data objects: vectors and matrices

- beginner

## One Way to Set the Elements of A Vector 

- _[VecGetSize](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetSize.html#VecGetSize)_(x,&N);  /* Global size */MPI_Comm_rank( PETSC_COMM_WORLD, &rank );if (rank == 0) {    for (i=0; i<N; i++)        	 		_[VecSetValues](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetValues.html#VecSetValues)_(x,1,&i,&v,INSERT_VALUES);}/* These two routines ensure that the data is distributed to the other processes */_[VecAssemblyBegin](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecAssemblyBegin.html#VecAssemblyBegin)_(x);_[VecAssemblyEnd](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecAssemblyEnd.html#VecAssemblyEnd)_(x);

- Vector index

- Vector value

## A Parallel Way to Set the Elements of A Distributed Vector 

- _[VecGetOwnershipRange](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecGetOwnershipRange.html#VecGetOwnershipRange)_(x,&low,&high); for (i=low; i<high; i++)        	 		_[VecSetValues](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetValues.html#VecSetValues)_(x,1,&i,&i,INSERT_VALUES);/* These two routines must be called (in case some other process contributed a value owned by another process) */_[VecAssemblyBegin](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecAssemblyBegin.html#VecAssemblyBegin)_(x);_[VecAssemblyEnd](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecAssemblyEnd.html#VecAssemblyEnd)_(x);

## Selected Vector Operations

- beginner

- data objects: vectors

## Simple Example Programs

- ex2.c 			- synchronized printing

*And many more examples ...*

-  Location:  petsc/src/sys/examples/tutorials/

- 1

- beginner

- 1

**E**

- data objects: vectors

- ex1.c, ex1f.F, ex1f90.F 	- basic vector routines
- ex3.c, ex3f.F           	             - parallel vector layout

- Location:  petsc/src/vec/examples/tutorials/

**E**

**E**

- 1

## A Complete PETSc Program

- #include _[petscvec](http://www-unix.mcs.anl.gov/petsc/petsc-current/include/petscvec.h.html)__[.h](http://www-unix.mcs.anl.gov/petsc/petsc-current/include/petscvec.h.html)_**int main(int argc,char **argv)**{  _[Vec](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/Vec.html#Vec)_ x;  int n = 20,ierr;  _[PetscTruth](http://www-unix.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscTruth.html#PetscTruth)_ flg;  PetscScalar one = 1.0, dot;
-   PetscInitialize(&argc,&argv,0,0);  PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);  VecCreate(PETSC_COMM_WORLD,&x);  VecSetSizes(x,PETSC_DECIDE,n);  VecSetFromOptions(x);  VecSet(&one,x);  VecDot(x,x,&dot);   PetscPrintf(PETSC_COMM_WORLD,"Vector length %dn",(int)dot);   VecDestroy(x);   PetscFinalize();  return 0;}

## Working With Local Vectors

- It is sometimes more efficient to directly access the storage for the local part of a PETSc Vec.
- E.g., for finite difference computations involving elements of the vector
- PETSc allows you to access the local storage with
- VecGetArray( Vec, double *[] );
- You must return the array to PETSc when you are done with it
- VecRestoreArray( Vec, double *[] )
- This allows PETSc to handle and data structure conversions
- For most common uses, these routines are inexpensive and do *not *involve a copy of the vector.

## Example of VecGetArray

- Vec      vec;
- Double *avec;...
- VecCreate(PETSC_COMM_SELF,&vec);VecSetSizes(vec,PETSC_DECIDE,n);VecSetFromOptions(vec);VecGetArray(vec,&avec);/* compute with avec directly, e.g., */
- PetscPrintf(PETSC_COMM_WORLD, “First element of local array of vec in each process is %f\n”, avec[0] );VecRestoreArray(vec,&avec);

## Matrices

- What are PETSc matrices?
- Fundamental objects for storing linear operators (e.g., Jacobians)
- Create matrices via
- MatCreate(...,Mat *)
- MPI_Comm - processes 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.
- MatSetFromOptions(Mat) lets you set the MatType at *runtime*.

- data objects: matrices

- beginner

## Matrices and Polymorphism

- Single user interface, e.g.,
- Matrix assembly
- MatSetValues()
- Matrix-vector multiplication
- MatMult()
- Matrix viewing
- MatView()
- Multiple underlying implementations
- AIJ, block AIJ, symmetric block AIJ, block diagonal, dense, matrix-free, etc.
- A matrix is defined by its *properties*, the operations that you can perform with it.
- Not by its data structures

- beginner

- data objects: matrices

## Matrix Assembly

- Same form as for PETSc Vectors:
- 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

## Matrix Assembly Example

- Mat      A;
- int        column[3], i;
- double value[3];
- ...
- MatCreate(PETSC_COMM_WORLD,                  PETSC_DECIDE,PETSC_DECIDE,n,n,&A);MatSetFromOptions(A);/* mesh interior */value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
- if (rank == 0) {  /* Only one process creates matrix */
-     for (i=1; i<n-2; i++) {
-         column[0] = i-1; column[1] = i; column[2] = i+1;
-         MatSetValues(A,1,&i,3,column,value,INSERT_VALUES);    }
- }
- /* also must set boundary points  (code for global row 0 and n-1 omitted) */
- MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

- simple 3-point stencil for 1D discretization

- beginner

- data objects: matrices

- Choose the global
- Size of the matrix

- Let PETSc decide how to allocate matrix across processes

## Parallel Matrix Distribution

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

- 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

## Matrix Assembly Example With Parallel Assembly

- Mat      A;
- int        column[3], i, start, end,istart,iend;
- double value[3];
- ...
- MatCreate(PETSC_COMM_WORLD,                  PETSC_DECIDE,PETSC_DECIDE,n,n,&A);MatSetFromOptions(A);MatGetOwnershipRange(A,&start,&end);  /* mesh interior */istart = start; if (start == 0) istart = 1;
- iend = end; if (iend == n-1) iend = n-2;
- value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
- for (i=istart; i<iend; i++) {
-     column[0] = i-1; column[1] = i; column[2] = i+1;
-     MatSetValues(A,1,&i,3,column,value,INSERT_VALUES);
- }
- /* also must set boundary points  (code for global row 0 and n-1 omitted) */
- MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

- simple 3-point stencil for 1D discretization

- beginner

- data objects: matrices

- Choose the global
- Size of the matrix

- Let PETSc decide how to allocate matrix across processes

## Why Are PETSc Matrices The Way They Are?

- No one data structure is appropriate for all problems
- Blocked and diagonal formats provide significant performance benefits
- PETSc provides a large selection of formats and makes it (relatively) easy to extend PETSc by adding new data structures
- Matrix assembly is difficult enough without being forced to worry about data partitioning
- PETSc provides parallel assembly routines
- Achieving high performance still requires making most operations local to a process but programs can be incrementally developed.
- Matrix decomposition by consecutive rows across processes is simple and makes it easier to work with other codes.
- For applications with other ordering needs, PETSc provides “Application Orderings” (AO), described later.

## Blocking: Performance Benefits

- 3D compressible Euler code
- Block size 5
- IBM Power2

- data objects: matrices

- beginner

- More issues discussed in full tutorials available via PETSc web site.

## Viewers

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

- tutorial outline: viewers

- beginner

- beginner

- intermediate

## 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

## Viewing Vector Fields

- VecView(Vec x,PetscViewer v);
- Default viewers
- ASCII (sequential):	PETSC_VIEWER_STDOUT_SELF
- ASCII (parallel):	PETSC_VIEWER_STDOUT_WORLD
- X-windows:	PETSC_VIEWER_DRAW_WORLD
- Default ASCII formats
- PETSC_VIEWER_ASCII_DEFAULT
- PETSC_VIEWER_ASCII_MATLAB
- PETSC_VIEWER_ASCII_COMMON
- PETSC_VIEWER_ASCII_INFO
- etc.

- viewers

- beginner

- Solution components,
- using runtime option
-        -snes_vecmonitor

- velocity: *u*

- velocity: *v*

- temperature: *T*

- vorticity:

## Viewing Matrix Data

- MatView(Mat A, PetscViewer 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

## More Viewer Info

- Viewer creation
- ViewerASCIIOpen()
- ViewerDrawOpen()
- ViewerSocketOpen()
- Binary I/O of vectors and matrices
- output: ViewerBinaryOpen()
- input: VecLoad(), MatLoad()

- viewers

- intermediate

## Solvers:  Usage Concepts

- Linear (SLES)
- Nonlinear (SNES)
- Timestepping (TS)

- Context variables
- Solver options
- Callback routines
- Customization

_**Solver Classes**_

_**Usage Concepts**_

- tutorial outline: solvers

## Linear PDE Solution

- PETSc

- Application
- Initialization

- Evaluation of* A* and *b*

- Post-
- Processing

- Solve
*Ax = b*

- PC

- KSP

- Linear Solvers (SLES)

- PETSc code

- User code

- Main Routine

- solvers:linear

- beginner

## Linear Solvers

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

- solvers:linear

- beginner

## Sample Linear Application:Exterior Helmholtz Problem 

- Imaginary

- Real

- Solution Components

*Collaborators: H. M. Atassi, D. E. Keyes, *
*                         **L. C. McInnes, R. Susan-Resiga*

- solvers:linear

- beginner

## Helmholtz: The Linear System

- Logically regular grid, parallelized with DAs
- Finite element discretization (bilinear quads)
- Nonreflecting exterior BC (via DtN map)
- Matrix sparsity structure (option: -mat_view_draw)

- solvers:linear

- beginner

## Linear Solvers (SLES)

- 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

*SLES: Scalable Linear Equations Solvers*

- tutorial outline: solvers:               linear

- beginner

- beginner

- intermediate

- intermediate

- advanced

- advanced

- beginner

- beginner

## Objects In PETSc

- How should a matrix be described in a program?
- Old way:
- Dense matrix:
- double precision A(10,10)
- Sparse matrix:
- integer ia(11), ja(max_nz)
- double precision a(max_nx)
- New way:
- Mat M
- Hides the choice of data structure
- Of course, the library still needs to represent the matrix with some choice of data structure, but this is an *implementation detail*
- Benefit
- Programs become independent of particular choices of data structure, making it easier to modify and adapt

## Operations In PETSc

- How should operations like “solve linear system” be described in a program?
- Old way:
- mpiaijgmres( ia, ja, a, comm, x, b, nlocal, nglobal, ndir, orthmethod, convtol, &its )
- New way:
- SLESSolve( sles, b, x, &its )
- Hides the choice of algorithm
- Algorithms are to operations as data structures are to objects
- Benefit
- Programs become independent of  particular choices of algorithm, making it easier to explore algorithmic choices and to adapt to new methods.
- In PETSc, operations have their own handle, called a *“context variable”*

## Context Variables

- 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)

- solvers:linear

- beginner

## Creating the SLES Context

- C/C++ version
- ierr = SLESCreate(PETSC_COMM_WORLD,&sles);
- Fortran version
- call SLESCreate(PETSC_COMM_WORLD,sles,ierr)
- Provides an** identical** user interface for all linear solvers
- uniprocess and parallel
- real and complex numbers

- solvers:linear

- beginner

## SLES Structure

- Each SLES object actually contains two other objects:
- KSP — Krylov Space Method
- The iterative method
- The context contains information on method parameters (e.g., GMRES search directions), work spaces, etc
- PC — Preconditioners
- Knows how to apply a preconditioner
- The context contains information on the preconditioner, such as what routine to call to apply it

## Linear Solvers in PETSc 2.0

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

- Block Jacobi
- Overlapping Additive Schwarz
- ICC, ILU via BlockSolve95
- ILU(k), LU (direct solve, sequential only)
- Arbitrary matrix
- etc.

_**Krylov**__*** ***__**Methods (KSP)**_

_**Preconditioners (PC)**_

- solvers:linear

- beginner

## 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(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&A);  MatSetFromOptions(A); 			/* (code to assemble matrix not shown) */
- VecCreate(PETSC_COMM_WORLD,&x);VecSetSizes(x,PETSC_DECIDE, n);VecSetFromOptions(x);
- VecDuplicate(x,&b);                                  /* (code to assemble RHS vector not shown)*/
- SLESCreate(PETSC_COMM_WORLD,&sles);
- SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN);
- SLESSetFromOptions(sles);
- SLESSolve(sles,b,x,&its);
- SLESDestroy(sles);

- solvers:linear

- beginner

- Indicate whether the preconditioner has the same nonzero pattern as the matrix *each time a system is solved.* This default works with *all* preconditioners.  Other values (e.g., SAME_NONZERO_PATTERN) can be used for particular preconditioners.  Ignored when solving only one system

## Basic Linear Solver Code (Fortran)

- SLES     sles
- Mat        A
- Vec        x, b
- integer  n, its, ierr
- call  MatCreate( PETSC_COMM_WORLD,PETSC_DECIDE,n,n,A,ierr )
- call  MatSetFromOptions( A, ierr )
- call  VecCreate( PETSC_COMM_WORLD,x,ierr )
- call  VecSetSizes( x, PETSC_DECIDE, n, ierr )
- call  VecSetFromOptions( x, ierr )
- call  VecDuplicate( x,b,ierr )
- call  SLESCreate(PETSC_COMM_WORLD,sles,ierr)
- call  SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN,ierr)
- call  SLESSetFromOptions(sles,ierr)
- call  SLESSolve(sles,b,x,its,ierr)
- call  SLESDestroy(sles,ierr)

- C     then assemble matrix and right-hand-side vector

- solvers:linear

- beginner

## Customization Options

- 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
- Procedural Interface
- Provides a great deal of control on a usage-by-usage basis inside a single code
- Gives full flexibility inside an application

- solvers:linear

- beginner

## -ksp_type  [cg,gmres,bcgs,tfqmr,...]
-pc_type  [lu,ilu,jacobi,sor,asm,...]

-ksp_max_it  <max_iters>
-ksp_gmres_restart  <restart>
-pc_asm_overlap  <overlap>
-pc_asm_type  [basic,restrict,interpolate,none]
etc ...

- -ksp_type  [cg,gmres,bcgs,tfqmr,...]
- -pc_type  [lu,ilu,jacobi,sor,asm,...]
- -ksp_max_it  <max_iters>
- -ksp_gmres_restart  <restart>
- -pc_asm_overlap  <overlap>
- -pc_asm_type  [basic,restrict,interpolate,none]
- etc ...

- Setting Solver Options at Runtime

- solvers:linear

- 1

- 2

## Linear Solvers: Monitoring Convergence

- -ksp_monitor          - Prints preconditioned residual norm
- -ksp_xmonitor         - Plots preconditioned residual norm
- -ksp_truemonitor     - Prints true residual norm || b-Ax ||
- -ksp_xtruemonitor   - Plots true residual norm || b-Ax ||
- User-defined monitors, using callbacks

- solvers:linear

- advanced

- 3

- 1

- 2

- 3

## Setting Solver Options within Code

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

- solvers:linear

- beginner

## Recursion: Specifying Solvers for Schwarz  Preconditioner Blocks

- Specify SLES solvers and options with “-sub” prefix, e.g.,
- Full or incomplete factorization
- -sub_pc_type lu
- -sub_pc_type ilu  -sub_pc_ilu_levels <levels>
- Can also use inner Krylov iterations, e.g.,
- -sub_ksp_type  gmres  -sub_ksp_rtol <rtol>
- -sub_ksp_max_it <maxit>

- solvers: linear: preconditioners

- beginner

## Helmholtz: Scalability

- 128x512 grid,  wave number = 13,  IBM SP
- GMRES(30)/Restricted Additive Schwarz
- 1 block per proc, 1-cell overlap, ILU(1) subdomain solver

- beginner

- solvers:linear

## 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

- beginner

- solvers:linear

## SLES: Review of Selected Preconditioner Options

*And many more options...*

- solvers: linear: preconditioners

- 1

- 2

## SLES: Review of Selected Krylov Method Options

- solvers: linear: Krylov methods

*And many more options...*

- 1

- 2

## Why Polymorphism?

- Programs become independent of the choice of algorithm
- Consider the question:
- What is the best combination of iterative method and preconditioner for my problem?
- How can you answer this experimentally?
- Old way:
- Edit code.  Make. Run.  Edit code. Make. Run. Debug. Edit. ...
- New way:...

## SLES: Runtime Script Example

- solvers:linear

- intermediate

## Viewing SLES Runtime Options

- solvers:linear

- intermediate

## Providing Different Matrices to Define Linear System and Preconditioner

- Krylov method:  Use *A *for matrix-vector products
- Build preconditioner using either
*A *- matrix that defines linear system
- or *P *- a different matrix (cheaper to assemble)
- SLESSetOperators(SLES sles,
- Mat A,
- Mat P,
- MatStructure flag)

- solvers:linear

- advanced

## Matrix-Free Solvers

- Use “shell” matrix data structure
- MatCreateShell(..., Mat *mfctx)
- Define operations for use by Krylov methods
- MatShellSetOperation(Mat mfctx,
- MatOperation MATOP_MULT,
- (void *) int (UserMult)(Mat,Vec,Vec))
- Names of matrix operations defined in petsc/include/mat.h
- Some defaults provided for nonlinear solver usage

- advanced

- solvers:linear

## User-defined Customizations

- Restricting the available solvers
- Customize PCRegisterAll( ), KSPRegisterAll( )
- Adding user-defined preconditioners via
- PCShell preconditioner type
- Adding preconditioner and Krylov methods in library style
- Method registration via PCRegister( ), KSPRegister( )
- Heavily commented example implementations
- Jacobi preconditioner: petsc/src/sles/pc/impls/jacobi.c
- Conjugate gradient: petsc/src/sles/ksp/impls/cg/cg.c

- solvers:linear

- 3

- 4

## SLES: Example Programs

- ex1.c, ex1f.F - basic uniprocess codes
- ex23.c           - basic parallel code
- ex11.c          - using complex numbers
- ex4.c              - using different linear system and
-                                preconditioner matrices
- ex9.c              - repeatedly solving different linear systems
- ex22.c	   - 3D Laplacian using multigrid
- ex15.c            - setting a user-defined preconditioner

*And many more examples ...*

- Location:  petsc/src/sles/examples/tutorials/

- solvers:linear

- 1

- 2

- 3

**E**

**E**

## Now What?

- If you have a running program, are you done?
- Yes, if your program answers your question.
- No, if you now need to run much larger problems or many more problems.
- How can you tune an application for performance?
- PETSc approach:
- Debug the application by focusing on the mathematically operations.  The parallel matrix assembly operations are an example.
- Once the code *works**, *
- Understand the performance
- Identify performance problems
- Use special PETSc routines and other techniques to optimize only the code that is underperforming

## Profiling and Performance Tuning

- Integrated profiling using -log_summary
- User-defined events
- Profiling by stages of an application

**Profiling:**

**Performance Tuning:**

- tutorial outline: profiling and performance tuning

- beginner

- intermediate

- intermediate

- Matrix optimizations
- Application optimizations
- Algorithmic tuning

- intermediate

- advanced

- intermediate

## Profiling

- Integrated monitoring of
- time
- floating-point performance
- memory usage
- communication
- All PETSc events are logged if PETSc was compiled with  -DPETSC_LOG (default); can also profile application code segments
- Print summary data with option:  -log_summary
- Print redundant information from PETSc routines: -log_info
- Print the trace of the functions called: -log_trace

- profiling and performance tuning

- beginner

## User-defined Events

- int USER_EVENT;
- int user_event_flops
- PetscLogEventRegister(&USER_EVENT,”User event name”, “eventColor”);
- PetscLogEventBegin(USER_EVENT,0,0,0,0);
- [ code to monitor]
- PetscLogFlops(user_evnet_flops);
- PetscLogEvent End(USER_EVENT,0,0,0,0);

- profiling and performance tuning

- intermediate

## Nonlinear Solvers (SNES)

*SNES: Scalable Nonlinear Equations Solvers*

- beginner

- beginner

- intermediate

- advanced

- advanced

- Application code interface
- Choosing the solver
- Setting algorithmic options
- Viewing the solver
- Determining and monitoring convergence
- Matrix-free solvers
- User-defined customizations

- tutorial outline: solvers:               nonlinear

- beginner

- beginner

## Nonlinear 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)

- Solve
*F(u) = 0*

- solvers: nonlinear

- beginner

## Nonlinear Solvers

**Goal**:  For problems arising from PDEs,
- support the general solution of  *F(u) = 0*
- 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
- AD support via collaboration with P. Hovland and B. Norris
- Coming in next PETSc release via automated interface to ADIFOR and ADIC (see _[http://www.mcs.anl.gov/](http://www.mcs.anl.gov/autodiff)__[autodiff](http://www.mcs.anl.gov/autodiff)_)

- solvers: nonlinear

- beginner

## Nonlinear Solvers (SNES)

- Newton-based methods, including
- Line search strategies
- Trust region approaches
- Pseudo-transient continuation
- Matrix-free variants
- User can customize all phases of the solution process

- solvers: nonlinear

- beginner

## Sample Nonlinear Application:Driven Cavity Problem

*Application code author: D. E. Keyes**                          *

- Velocity-vorticity formulation
- Flow driven by lid and/or bouyancy
- Logically regular grid, parallelized with DAs
- Finite difference discretization
- source code:

- solvers: nonlinear

- beginner

- petsc/src/snes/examples/tutorials/ex19.c

## Basic Nonlinear Solver Code (C/C++)

- 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  */
- ...
- MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,&J);MatSetFromOptions(J);
- VecCreate(PETSC_COMM_WORLD,&x);VecSetSizes(x,PETSC_DECIDE,n);
- VecSetFromOptions(x);
- VecDuplicate(x,&F);
- SNESCreate(PETSC_COMM_WORLD,&snes);
- SNESSetFunction(snes,F,EvaluateFunction,usercontext);
- SNESSetJacobian(snes,J,J,EvaluateJacobian,usercontext);
- SNESSetFromOptions(snes);
- SNESSolve(snes,x,&its);
- SNESDestroy(snes);

- solvers: nonlinear

- beginner

## Basic Nonlinear Solver Code (Fortran)

-       SNES   snes
-       Mat      J
-       Vec      x, F
-       int        n, its
-          **...**
-       call  MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE, PETSC_DECIDE,n,n,J,ierr)
-       call  MatSetFromOptions( J, ierr )
-       call  VecCreate(PETSC_COMM_WORLD,x,ierr)
-       call  VecSetSizes(x,PETSC_DECIDE,n,ierr)
-       call  VecSetFromOptions(x,ierr)
-       call  VecDuplicate(x,F,ierr)
-       call  SNESCreate(PETSC_COMM_WORLD, snes,ierr)
-       call  SNESSetFunction(snes,F,EvaluateFunction,PETSC_NULL,ierr)
-       call  SNESSetJacobian(snes,J,J,EvaluateJacobian,PETSC_NULL,ierr)
-       call  SNESSetFromOptions(snes,ierr)
-       call  SNESSolve(snes,x,its,ierr)
-       call  SNESDestroy(snes,ierr)

- solvers: nonlinear

- beginner

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

- solvers: nonlinear

- beginner

## Sample Application Context: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;

- solvers: nonlinear

- beginner

## Sample Function Evaluation Code: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 */
-      ....
-      /* (Code to communicate nonlocal ghost point data not shown) */
-      VecGetArray( F, &f );
-      /* (Code to compute local function components; insert into f[ ] shown on next slide) */
-      VecRestoreArray( F, &f );
-      ....
- return 0;
- }

- solvers: nonlinear

- beginner

## Sample Local Computational Loops: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;
-  } }
-  ....

- solvers: nonlinear

-   The PDE’s 4 components (U,V,Omega,Temp) are interleaved in the unknown vector.
-    #define statements provide easy access to each component.

- beginner

## Finite Difference Jacobian Computations

- Compute and explicitly store Jacobian via 1st-order FD
- Sparse via colorings: MatFDColoringCreate(), SNESDefaultComputeJacobianColor()
- Matrix-free Newton-Krylov via 1st-order FD, no preconditioning unless specifically set by user
- -snes_mf
- Matrix-free Newton-Krylov via 1st-order FD, user-defined preconditioning matrix
- -snes_mf_operator

- beginner

- solvers: nonlinear

## Uniform access to all linear and nonlinear solvers

- -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_convergence <tolerance>
- etc...

- solvers: nonlinear

- 1

- 2

## Customization via Callbacks: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( )

- solvers: nonlinear

- 2

- 3

## SNES:  Review of Basic Usage

- 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

- beginner

- solvers: nonlinear

## SNES: Review of Selected Options

*And many more options...*

- 1

- 2

- solvers: nonlinear

## SNES: Example Programs

- ex1.c, ex1f.F		     - basic uniprocess codes
- ex2.c 		         	     - uniprocess nonlinear PDE 					        (1 DoF per node)
- ex5.c, ex5f.F, ex5f90.F  	     - parallel nonlinear PDE (1 DoF per node)
- ex18.c		     - parallel radiative transport problem with  			       multigrid
- ex19.c          		     - parallel driven cavity problem with 				       multigrid

*And many more examples ...*

- Location:  petsc/src/snes/examples/tutorials/

- 1

- 2

- beginner

- 1

- intermediate

- 2

- solvers: nonlinear

**E**

**E**

**E**

## Timestepping Solvers (TS)

- (and ODE Integrators)

- tutorial outline: solvers: timestepping

- beginner

- beginner

- intermediate

- advanced

- Application code interface
- Choosing the solver
- Setting algorithmic options
- Viewing the solver
- Determining and monitoring convergence
- User-defined customizations

- beginner

- beginner

## Time-Dependent 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)

- Solve
* **U** t **= F(U,U**x,**U**xx**)*

- solvers: timestepping

- beginner

## Timestepping Solvers

**Goal**:  Support the (real and pseudo) time
- evolution of PDE systems
*                      **U**t**  = F(U,U**x**,U**xx**,t)*
- User provides:
- Code to evaluate *F(U,U**x**,U**xx**,t)*
- Code to evaluate Jacobian of* F(U,U**x**,U**xx**,t)*
- or use sparse finite difference approximation
- or use automatic differentiation (coming soon!)

- solvers: timestepping

- beginner

## Sample Timestepping Application:Burger’s Equation

- Sample Timestepping Application:Burger’s Equation

- solvers: timestepping

- beginner

## Ut = F(t,U) = Ui (Ui+1 - U i-1)/(2h) +
			    (Ui+1 - 2Ui + U i-1)/(h*h)

Do 10, i=1,localsize
F(i) = (.5/h)*u(i)*(u(i+1)-u(i-1)) +
	  (e/(h*h))*(u(i+1) - 2.0*u(i) + u(i-1))
10 continue  

*	**	**	**    ****(U**i+1** - 2U**i** + U **i-1**)/(h*h)*
- Do 10, i=1,localsize
- F(i) = (.5/h)*u(i)*(u(i+1)-u(i-1)) +
- 	  (e/(h*h))*(u(i+1) - 2.0*u(i) + u(i-1))
- 10 continue

- Actual Local Function Code

- solvers: timestepping

- beginner

## Timestepping Solvers

- Euler
- Backward Euler
- Pseudo-transient continuation
- Interface to PVODE, a sophisticated parallel ODE solver package by Hindmarsh et al. of LLNL
- Adams
- BDF

- solvers: timestepping

- beginner

## Timestepping Solvers

- Allow full access to all of the PETSc
- nonlinear solvers
- linear solvers
- distributed arrays, matrix assembly tools, etc.
- User can customize all phases of the solution process

- solvers: timestepping

- beginner

## TS: Review of Basic Usage

- TSCreate( )		- Create TS context
- TSSetRHSFunction( ) 	- Set function eval. routine
- TSSetRHSJacobian( ) 	- Set Jacobian eval. routine
- TSSetFromOptions( ) 	- Set runtime solver options   				   	   for [TS,SNES,SLES,KSP,PC]
- TSSolve( ) 		- Run timestepping solver
- TSView( ) 			- View solver options 					   	   actually used at runtime 				  	   (alternative: -ts_view)
- TSDestroy( ) 		- Destroy solver

- beginner

- solvers: timestepping

## TS: Review of Selected Options

- solvers: timestepping

- 1

- 2

*And many more options...*

## TS: Example Programs

- ex1.c, ex1f.F 	- basic uniprocess codes (time-				  dependent nonlinear PDE)
- ex2.c, ex2f.F 	- basic parallel codes (time-dependent 			   nonlinear PDE)
- ex3.c		- uniprocess heat equation
- ex4.c		- parallel heat equation

*And many more examples ...*

- Location:  petsc/src/ts/examples/tutorials/

- 1

- 2

- beginner

- 1

- intermediate

- 2

- solvers: timestepping

**E**

## Mesh Definitions:  For Our Purposes

**Structured**:  Determine neighbor relationships purely from logical I, J, K coordinates
**Semi-Structured**:  In well-defined regions, determine neighbor relationships purely from logical I, J, K coordinates
**Unstructured**:** ** Do not explicitly use logical I, J, K coordinates

- data layout

- beginner

## Structured Meshes

- PETSc support provided via DA objects

- data layout

- beginner

## Unstructured Meshes

- One is always free to manage the mesh data as if unstructured
- PETSc does not currently have *high-level* tools for managing such meshes
- PETSc can manage unstructured meshes though lower-level VecScatter utilities

- data layout

- beginner

## Semi-Structured Meshes

- No explicit PETSc support
- OVERTURE-PETSc for composite meshes
- SAMRAI-PETSc for AMR

- data layout

- beginner

## Parallel Data Layout and Ghost Values:  Usage Concepts

- Structured
- DA objects
- Unstructured
- VecScatter objects

- Geometric data
- Data structure creation
- Ghost point updates
- Local numerical computation

_**Mesh Types**_

_**Usage Concepts**_

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

- tutorial outline: data layout

## Ghost Values 

- data layout

- beginner

**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

- data layout

**Communication**

- Data Structure
- Creation

- Ghost Point
- Data Structures

- Ghost Point
- Updates

**Local**
**Numerical**
**Computation**

- Geometric
- Data

- DA
- AO

- DACreate(  )

- DAGlobalToLocal(  )

- Loops over
- I,J,K
- indices

- stencil
- [implicit]

- VecScatter
- AO

- VecScatterCreate(  )

- VecScatter(  )

- Loops over
- entities

- elements
- edges
- vertices

***unstructured meshes***

***structured meshes***

- 1

- 2

## DA: Parallel Data Layout and Ghost Values for Structured Meshes

- Local and global indices
- Local and global vectors
- DA creation
- Ghost point updates
- Viewing

- tutorial outline: data layout: distributed arrays

- beginner

- beginner

- intermediate

- intermediate

- beginner

## Communication and Physical Discretization:Structured Meshes

- data layout: distributed arrays

- beginner

**Communication**

- Data Structure
- Creation

- Ghost Point
- Data Structures

- Ghost Point
- Updates

**Local**
**Numerical**
**Computation**

- Geometric
- Data

- DA
- AO

- DACreate(  )

- DAGlobalToLocal(  )

- Loops over
- I,J,K
- indices

- stencil
- [implicit]

***structured meshes***

## Global and Local Representations 

- data layout: distributed arrays

- beginner

***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

## Global and Local Representations (cont.)

**Global Representation:**

**Local Representations:**

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 9

- 10

- 11

- Proc 1

- Proc 0

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 1

- 2

- 3

- 4

- 5

- 6

- 0

- 7

- 8

- Proc 0

**  ****6  7   8 **
**  ****3  4   5**
**  ****0   1   2 **
**  ****6   7   8**
**  ****3   4   5**
**  ****0   1   2**

**  ****9  10 11**
**  ****6   7   8**
**  ****3   4   5**
**  ****0   1   2**

- Proc 0

- Proc 1

- Proc 1

- Proc 0

- beginner

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 3

- 4

- 5

- 6

- 7

- 8

- 9

- 10

- 11

- Proc 1

- data layout: distributed arrays

## Logically Regular Meshes

- beginner

- DA - Distributed Array: object containing information about vector layout across the processes and communication of ghost values
- Form a DA
- DACreate1d(....,DA *)
- DACreate2d(....,DA *)
- DACreate3d(....,DA *)
- Create the corresponding PETSc vectors
- DACreateGlobalVector( DA, Vec *) or
- DACreateLocalVector( DA, Vec *)
- Update ghostpoints (scatter global vector into local parts, including ghost points)
- DAGlobalToLocalBegin(DA, ...)
- DAGlobalToLocalEnd(DA,...)

## Distributed Arrays

- Proc 10

- Proc 0

- Proc 1

- Proc 10

- Proc 0

- Proc 1

*Box-type *
*stencil*

*Star-type *
*stencil*

- data layout: distributed arrays

- beginner

- Data layout and ghost values

## Vectors and DAs

- 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 process stores its local portion plus ghost values
- DACreateLocalVector(DA da,Vec *lvec);
- uses “natural” local numbering of indices (0,1,...*nlocal*-1)

- data layout: distributed arrays

- beginner

## DACreate1d(...,*DA)

- DACreate1d(MPI_Comm comm,DAPeriodicType wrap,int M,int dof,int s,int *lc,DA *inra)
- MPI_Comm — processes containing array
- DA_[NONPERIODIC,XPERIODIC]
- number of grid points in x-direction
- degrees of freedom per node
- stencil width
- Number of nodes for each cell (use PETSC_NULL for the default)

- data layout: distributed arrays

- beginner

## DACreate2d(...,*DA)

- DACreate2d(MPI_Comm comm,DAPeriodicType wrap,DAStencilType stencil_type, int M,int N,int m,int n,int dof,int s,int *lx,int *ly,DA *inra)
- DA_[NON,X,Y,XY]PERIODIC
- DA_STENCIL_[STAR,BOX]
- number of grid points in x- and y-directions
- processes in x- and y-directions
- degrees of freedom per node
- stencil width
- Number of nodes for each cell (use PETSC_NULL for the default) as tensor-product

- data layout: distributed arrays

- beginner

- And similarly for DACreate3d()

## Updating the Local Representation

- DAGlobalToLocalBegin(DA, Vec global_vec, insert,                                       Vec local_vec )
- global_vec provides data
- Insert is either INSERT_VALUES or ADD_VALUES and specifies how to update values in the local vector, local_vec (a pre-existing local vector)
- DAGlobalToLocal End(DA,...)
- Takes same arguments

- Two-step process that enables overlapping computation and communication

- data layout: distributed arrays

- beginner

## Ghost Point Scatters: Burger’s Equation Example

-      call DAGlobalToLocalBegin(da,u_global,INSERT_VALUES,ulocal,ierr)
-       call DAGlobalToLocalEnd(da,u_global,INSERT_VALUES,ulocal,ierr)
-       call VecGetArray( ulocal, uv, ui, ierr )
- #define u(i) uv(ui+i)
- C  Do local computations (here u and f are local vectors)
-       do 10, i=1,localsize
-      f(i) = (.5/h)*u(i)*(u(i+1)-u(i-1)) +
-      	  (e/(h*h))*(u(i+1) - 2.0*u(i) + u(i-1))
- 10  continue      call VecRestoreArray( ulocal, uv, ui, ierr )
-       call DALocalToGlobal(da,f,INSERT_VALUES,f_global,ierr)

- beginner

- data layout: distributed arrays

## Global Numbering used by DAs

- Natural numbering, corresponding to the entire problem domain

- 26  27  28      29  30
- 21  22  23      24  25
- 16  17  18      19  20
- 11  12  13      14  15
-   6    7    8        9  10
-   1    2    3        4    5

- Proc 0

- Proc 3

- Proc 2

- Proc1

- 22  23  24      29  30
- 19  20  21      27  28
- 16  17  18      25  26
-   7    8    9      14  15
-   4    5    6      12  13
-   1    2    3      10  11

- Proc 0

- Proc 3

- Proc 2

- Proc1

- PETSc numbering used by DAs

- data layout: distributed arrays

- intermediate

## Mapping Between Global Numberings

- Natural global numbering
- convenient for visualization of global problem, specification of certain boundary conditions, etc.
- Can convert between various global numbering schemes using AO (Application Orderings)
- DAGetAO(DA da, AO *ao);
- AO usage explained in next section
- Some utilities (e.g., VecView()) automatically handle this mapping for global vectors attained from DAs

- data layout: distributed arrays

- intermediate

## Distributed Array Example

- src/snes/examples/tutorial/ex5.c, ex5f.F
- Uses functions that help construct vectors and matrices naturally associated a DA
- DAGetMatrix
- DASetLocalFunction
- DASetLocalJacobian

## Unstructured Meshes

- beginner

- Setting up communication patterns is much more complicated than the structured case due to
- mesh dependence
- discretization dependence
- cell-centered
- vertex-centered
- cell and vertex centered (e.g., staggered grids)
- mixed triangles and quadrilaterals
- Can use VecScatter
- Introduction in this tutorial
- See additional tutorial material available via PETSc web site

## Communication and Physical Discretization

- data layout

## Unstructured Mesh Concepts

- data layout: vector scatters

- AO: Application Orderings
- map between various global numbering schemes
- IS: Index Sets
- indicate collections of nodes, cells, etc.
- VecScatter:
- update ghost points using vector scatters/gathers
- ISLocalToGlobalMapping
- map from a local (calling process) numbering of indices to a global (across-processes) numbering

## Setting Up the Communication Pattern

- Renumber objects so that contiguous entries are adjacent
- Use AO to map from application to PETSc ordering
- Determine needed neighbor values
- Generate local numbering
- Generate local and global vectors
- Create communication object (VecScatter)

- data layout: vector scatters

- intermediate

- (the steps in creating VecScatter)

## Cell-based Finite VolumeApplication Numbering

- data layout: vector scatters

- intermediate

- global indices defined by application

- Proc 0

- Proc1

## PETSc Parallel Numbering

- data layout: vector scatters

- intermediate

- global indices numbered contiguously on each process

- Proc 0

- Proc 1

## Remapping Global Numbers: An Example

- Process 0
- nlocal: 6
- app_numbers: {1,2,0,5,11,10}
- petsc_numbers: {0,1,2,3,4,5}
- Process 1
- n_local: 6
- app_numbers: {3,4,6,7,9,8}
- petsc_numbers: {6,7,8,9,10,11}

- data layout: vector scatters

- intermediate

- Proc 0

- Proc 1

- PETSc
- numbers

- application
- numbers

## Remapping Numbers (1)

- Generate a parallel object (AO) to use in mapping between numbering schemes
- AOCreateBasic( MPI_Comm comm,
- int nlocal,
- int app_numbers[ ],
- int petsc_numbers[ ],
- AO *ao );

- data layout: vector scatters

- intermediate

## Remapping Numbers (2)

- AOApplicationToPetsc(AO *ao,
- int number_to_convert,
- int indices[ ]);
- For example, if  indices[ ] contains the cell neighbor lists in an application numbering, then apply AO to convert to new numbering

- data layout: vector scatters

- intermediate

## Neighbors

- data layout: vector scatters

- intermediate

- Process 0

- Process1

- 0

- 1

- 6

- 7

- 2

- 3

- 8

- 9

- 11

- 10

- 5

- 4

- 3

- 1

- 5

- 6

- 8

- 10

- ghost cells

- using global PETSc numbers (ordering of the global PETSc vector)

## Local Numbering

- data layout: vector scatters

- intermediate

- Proc 0

- Proc1

- 0

- 1

- 0

- 1

- 2

- 3

- 2

- 3

- 5

- 4

- 5

- 4

- 7

- 6

- 8

- 6

- 7

- 8

- ghost cells

- At *end* of local vector

## Summary of the Different Orderings

- Local points at the *end*

## Global and Local Representations

- Global representation
- parallel vector with no ghost locations
- suitable for use by PETSc parallel solvers (SLES, SNES, and TS)
- Local representation
- sequential vectors with room for ghost points
- used to evaluate functions, Jacobians, etc.

- data layout: vector scatters

- intermediate

## Global and Local Representations

- Local Representations:

- 0

- 1

- 2

- 3

- 4

- 5

- 6

- 7

- 8

- 6

- 7

- 8

- 9

- 10

- 11

- 1

- 3

- 5

- data layout: vector scatters

- intermediate

- Proc1

- Proc 0

- Proc 0

- Proc1

## Creating Vectors

- Sequential
-     VecCreateSeq(PETSC_COMM_SELF, 9, Vec *lvec);
- Parallel
-     VecCreateMPI(PETSC_COMM_WORLD,
-                                6,PETSC_DETERMINE,Vec *gvec)

- data layout: vector scatters

- intermediate

## Local and Global Indices

- Process 0
- int indices[] = {0,1,2,3,4,5,6,8,10};ISCreateGeneral(PETSC_COMM_WORLD,                              9, indices, IS *isg);
- ISCreateStride(PETSC_COMM_SELF, 9,0,1,IS* isl);
- Process 1
- int indices = {6,7,8,9,10,11,1,3,5};ISCreateGeneral(PETSC_COMM_WORLD,                              9, indices, IS *isg);
- ISCreateStride(PETSC_COMM_SELF, 9,0,1,IS* isl);

- data layout: vector scatters

- intermediate

- ISCreateGeneral() - Specify *global* numbers of locally owned cells, including ghost cells

- ISCreateStride() - Specify *local *numbers of locally owned cells, including ghost cells

## Creating Communication Objects

- VecScatterCreate(Vec gvec,
- IS gis,
- Vec lvec,
- IS lis
- VecScatter gtol);
- Determines all required messages for mapping data from a global vector to local (ghosted) vectors

- data layout: vector scatters

- intermediate

## Performing a Global-to-Local Scatter

- VecScatterBegin(VecScatter gtol,
- Vec gvec,
- Vec lvec,
- INSERT_VALUES
- SCATTER_FORWARD);
- VecScatterEnd(...);

- Two-step process that enables overlapping computation and communication

- data layout: vector scatters

- intermediate

## Performing a Local-to-Global Scatter

- VecScatterBegin(VecScatter gtol,
- Vec lvec,
- Vec gvec,
- ADD_VALUES,
- SCATTER_REVERSE);
- VecScatterEnd(...);

- data layout: vector scatters

- intermediate

## Setting Values in Global Vectors and Matrices using a Local Numbering

- Create mapping
- ISLocalToGlobalMappingCreateIS(IS gis, ISLocalToGlobalMapping *lgmap);
- Set mapping
- VecSetLocalToGlobalMapping(Vec gvec, ISLocalToGlobalMapping lgmap);
- MatSetLocalToGlobalMapping(Mat gmat, ISLocalToGlobalMapping lgmap);
- Set values with local numbering
- VecSetValuesLocal(Vec gvec,int ncols, int localcolumns, Scalar *values, ...);
- MatSetValuesLocal(Mat gmat, ...);
- MatSetValuesLocalBlocked(Mat gmat, ...);

- data layout: vector scatters

- intermediate

## Sample Function Evaluation

- data layout: vector scatters

- intermediate

- int FormFunction(SNES snes, Vec Xglobal, Vec Fglobal, void *ptr)
- {
-     AppCtx  *user = (AppCtx *) ptr;
-     double   x1, x2, f1, f2, *x, *f ;
-     int          *edges = user->edges;
-     VecScatterBegin(user->scatter, Xglobal, user->Xlocal, SCATTER_FORWARD, INSERT_VALUES);
-     VecScatterEnd(user->scatter, Xglobal, user->Xlocal, SCATTER_FORWARD, INSERT_VALUES);
-     VecGetArray(Xlocal,&X);   VecGetArray(Flocal,&F);
-     for (i=0; i < user->nlocal; i++)  {
-          x1 = X[edges[2*i]];  x2 = X[edges[2*i+1]];  /* then compute f1, f2 */
-          F[edges[2*i]] += f1;  F[edges[2*i+1]] += f2;
-     }
-     VecRestoreArray(Xlocal,&X);   VecRestoreArray(Flocal,&F);
-     VecScatterBegin(user->scatter, user->Flocal, Fglobal, SCATTER_REVERSE, INSERT_VALUES);
-     VecScatterEnd(user->scatter, user->Flocal, Fglobal, SCATTER_REVERSE, INSERT_VALUES);
-     return 0;
- }

## Communication and Physical Discretization

- data layout

## Unstructured Meshes: Example Programs

- Location:
- petsc/src/snes/examples/tutorials/ex10d/ex10.c
- 2 process example (PETSc requires that a separate tool be used to divide the unstructured mesh among the processes).

## Matrix Partitioning and Coloring

- Partitioner object creation and use
- Sparse finite differences for Jacobian computation (using colorings)

- tutorial outline:     data objects: matrix partitioning and coloring

- Intermediate

- intermediate

## Partitioning for Load Balancing

- MatPartitioning - object for managing the partitioning of meshes, vectors, matrices, etc.
- Interface to parallel partitioners, i.e., the adjacency matrix, is stored in parallel
- MatPartitioningCreate(
- MPI_Comm comm
- MatPartitioning *matpart);

- data objects: matrices: partitioning

- intermediate

## Setting Partitioning Software

- MatPartitioningSetType(
- MatPartitioning matpart,
- MatPartitioningType MATPARTITIONING_PARMETIS);
- PETSc interface could support other packages (e.g., Pjostle, just no one has written wrapper code yet)

- data objects: matrices: partitioning

- intermediate

## Setting Partitioning Information

- MatPartitioningSetAdjacency(
- MatPartitioning matpart,
- Mat adjacency_matrix);
- MatPartitioningSetVertexWeights(
- MatPartitioning matpart,
- int *weights);
- MatPartitioningSetFromOptions(
- MatPartitioning matpart);

- intermediate

- data objects: matrices: partitioning

## Partitioning

- MatPartitioningApply(
- MatPartitioning matpart,
- IS* processforeachlocalnode);
- ...
- MatPartitioningDestroy(
- MatPartitioning matpart);

- intermediate

- data objects: matrices: partitioning

## Constructing Adjacency Matrix

- MatCreateMPIAdj( (Uses CSR input format)
- MPI_Comm comm,
- int numberlocalrows,
- int numbercolums,
- int rowindices[ ],
- int columnindices[ ],
- int values[ ],
- Mat *adj)
- Other input formats possible, e.g., via MatSetValues( )

- intermediate

- data objects: matrices: partitioning

## Finite Difference Approximation of Sparse Jacobians

- Compute a coloring of the Jacobian (e.g., determine columns of the Jacobian that may be computed using a single function evaluation)
- Construct a MatFDColoring object that uses the coloring information to compute the Jacobian
- example: driven cavity model: petsc/src/snes/examples/tutorials/ex8.c

- data objects: matrices: coloring

- Two-stage process:

- intermediate

## Coloring for Structured Meshes

- DAGetColoring (
- DA da,
- ISColorType ctype,
- ISColoring *coloringinfo );
- DAGetMatrix(
- DA da,
- MatType mtype,
- Mat *sparsematrix );
- For structured meshes (using DAs) we provide parallel colorings

- intermediate

- data objects: matrices: coloring

## Coloring for Unstructured Meshes

- MatGetColoring (
- Mat A,
- MatColoringType
- MATCOLORING_NATURAL
- MATCOLORING_SL
- MATCOLORING_LF
- MATCOLORING_ID
- ISColoring *coloringinfo)
- Automatic coloring of sparse matrices currently implemented only for sequential matrices
- If using a “local” function evaluation, the sequential coloring is enough

- intermediate

- data objects: matrices: coloring

## Actual Computation of Jacobians

- MatFDColoringCreate (
- Mat J,
- ISColoring coloringinfo,
- MatFDColoring *matfd)
- MatFDColoringSetFunction(
- MatFDColoring matfd,
- int (*f)(void),void *fctx)
- MatFDColoringSetFromOptions(
- MatFDColoring)

- intermediate

- data objects: matrices: coloring

## Computation of Jacobians within SNES

- User calls:  SNESSetJacobian(
- SNES snes,
- Mat A ,
- Mat J,
- SNESDefaultComputeJacobianColor( ),
- fdcoloring);
- where SNESSetJacobian() actually uses .... MatFDColoringApply(
- Mat J,
- MatFDColoring coloring,
- Vec x1,
- MatStructure *flag,
- void *sctx)
- to form the Jacobian

- data objects: matrices: coloring

- intermediate

## Driven Cavity Model

- Velocity-vorticity formulation, with flow driven by lid and/or bouyancy
- Finite difference discretization with 4 DoF per mesh point

- solvers: nonlinear

- Example code: petsc/src/snes/examples/tutorials/ex19.c

*[u,v*,,*T]*

- beginner

- 1

- intermediate

- 2

## Driven Cavity Program

**Part A**: Parallel data layout** **
**Part***** *****B**: Nonlinear solver creation, setup, and usage
**Part C**: Nonlinear function evaluation
- ghost point updates
- local function computation
**Part D**: Jacobian evaluation
- default colored finite differencing approximation
- Experimentation

- beginner

- 1

- intermediate

- 2

- solvers: nonlinear

## Driven Cavity Solution Approach

- PETSc code

- User code

- Application
- Initialization

- Function
- Evaluation

- Jacobian
- Evaluation

- Post-
- Processing

- PC

- KSP

- PETSc

- Main Routine

- Linear Solvers (SLES)

- Nonlinear Solvers (SNES)

- Solve
*F(u) = 0*

- solvers: nonlinear

**A**

**C**

**D**

**B**

## Driven Cavity: Running the program (1) 

- 1 process: (thermally-driven flow)
- mpirun -np 1 ex19 -snes_mf -snes_monitor -grashof 1000.0 -lidvelocity 0.0
- 2 processes, view DA (and pausing for mouse input):
- mpirun -np 2 ex19 -snes_mf -snes_monitor                                                         -da_view_draw  -draw_pause  -1
- View contour plots of converging iterates
- mpirun ex19 -snes_mf -snes_monitor  -snes_vecmonitor

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

- solvers: nonlinear

- beginner

## Driven Cavity: 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

- solvers: nonlinear

- intermediate

## Debugging and Error Handling

- Automatic generation of tracebacks
- Detecting memory corruption and leaks
- Optional user-defined error handlers

- tutorial outline: debugging and errors

- beginner

- beginner

- developer

## Debugging

- -start_in_debugger  [gdb,dbx,noxterm]
- -on_error_attach_debugger [gb,dbx,noxterm]
- -on_error_abort
- -debugger_nodes 0,1
- -display machinename:0.0

- beginner

- Support for parallel debugging

- When debugging, it is often useful to place a breakpoint in the function PetscError( ).

- debugging and errors

## Sample Error Traceback

- Breakdown in ILU factorization due to a zero pivot

- debugging and errors

- beginner

## Sample Memory Corruption Error

- beginner

- debugging and errors

## Sample Out-of-Memory Error

- beginner

- debugging and errors

## Sample Floating Point Error

- beginner

- debugging and errors

## Matrix Memory Pre-allocation

- PETSc sparse matrices are dynamic data structures. Can add additional nonzeros freely
- Dynamically adding many nonzeros
- requires additional memory allocations
- requires copies
- can kill performance
- Memory pre-allocation provides the freedom of dynamic data structures plus good performance

- intermediate

- profiling and performance tuning

## Indicating Expected NonzerosSequential Sparse Matrices

- Indicating Expected NonzerosSequential Sparse Matrices

- MatCreateSeqAIJ(...., int *nnz,Mat *A)
- nnz[0] - expected number of nonzeros in row 0
- nnz[1] - expected number of nonzeros in row 1

- row 0
- row 1
- row 2
- row 3

- another sample
- nonzero pattern

- sample nonzero pattern

- intermediate

- profiling and performance tuning

## Symbolic Computation of Matrix Nonzero Structure

- Write code that forms the nonzero structure of the matrix
- loop over the grid for finite differences
- loop over the elements for finite elements
- etc.
- Then create matrix via MatCreateSeqAIJ(), MatCreateMPIAIJ(), ...

- intermediate

- profiling and performance tuning

## Parallel Sparse Matrices

- proc 3: locally owned rows

- proc 0

- proc 3

- proc 2

- proc 1

- proc 4

- diagonal portions

- off-diagonal portions

- Each process locally owns a submatrix of contiguously numbered global rows.
- Each submatrix consists of **diagonal **and **off-diagonal** parts.

- intermediate

- profiling and performance tuning

## Indicating Expected NonzerosParallel Sparse Matrices

- MatCreateMPIAIJ(...., int d_nz, int *d_nnz,int o_nz,
-                               int *o_nnz,Mat *A)
- d_nnz[ ] - expected number of nonzeros per row in diagonal portion of local submatrix
- o_nnz[ ] - expected number of nonzeros per row in off-diagonal portion of local submatrix

- intermediate

- profiling and performance tuning

## Verifying Predictions

- Use runtime option:  -log_info
- Output:
-        [proc #] Matrix size: % d X % d; storage space: % d unneeded, % d used
-       [proc #] Number of mallocs during MatSetValues( )  is % d

- intermediate

- profiling and performance tuning

## Other Features of PETSc

- Summary
- New features
- Interfacing with other packages
- Extensibility issues
- References

## Summary

- Creating data objects
- Setting algorithmic options for linear, nonlinear and ODE solvers
- Using callbacks to set up the problems for nonlinear and ODE solvers
- Managing data layout and ghost point communication
- Evaluating parallel functions and Jacobians
- Consistent profiling and error handling

## General multigrid support
PC framework wraps MG for use as preconditioner 
See MGSetXXX(), MGGetXXX()
Can access via -pc_type mg
User provides coarse grid solver, smoothers, and interpolation/restriction operators
DMMG - simple MG interface for structured meshes
User provides
“Local” function evaluation 
[Optional] local Jacobian evaluation

- PC framework wraps MG for use as preconditioner
- See MGSetXXX(), MGGetXXX()
- Can access via -pc_type mg
- User provides coarse grid solver, smoothers, and interpolation/restriction operators
- DMMG - simple MG interface for structured meshes
- User provides
- “Local” function evaluation
- [Optional] local Jacobian evaluation

- Multigrid Structured Mesh Support:
- DMMG: New Simple Interface

## int Function(DALocalInfo *info,double **u,double **f,AppCtx *user)
  ...
  lambda = user->param;
  hx     = 1.0/(info->mx-1);
  hy     = 1.0/(info->my-1);
  for (j=info->ys; j<info->ys+info->ym; j++) {
    for (i=info->xs; i<info->xs+info->xm; i++) {
       f[j][i] = ... u[j][i] ... u[j-1][i] ....
    }
  }

-   ...
-   lambda = user->param;
-   hx     = 1.0/(info->mx-1);
-   hy     = 1.0/(info->my-1);
-   for (j=info->ys; j<info->ys+info->ym; j++) {
-     for (i=info->xs; i<info->xs+info->xm; i++) {
-        f[j][i] = ... u[j][i] ... u[j-1][i] ....
-     }
-   }

- Multigrid Structured Mesh Support:
- Sample Function Computation

## int Jacobian (DALocalInfo *info,double **u,Mat J,AppCtx *user)
  MatStencil mrow,mcols[5];
  double       v[5];
  ...
  for (j=info->ys; j<info->ys+info->ym; j++) {
    row.j = j;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        v[0] = ...;                                           col[0].j = j - 1;    col[0].i = i;
        v[1] = ...;                                           col[1].j = j;         col[1].i = i-1;
        v[2] = ...;                                           col[2].j = j;         col[2].i = i;
        v[3] = ...;                                           col[3].j = j;         col[3].i = i+1;
        v[4] = ...;                                           col[4].j = j + 1;   col[4].i = i;                         
        MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
     }
   }

-   MatStencil mrow,mcols[5];
-   double       v[5];
-   ...
-   for (j=info->ys; j<info->ys+info->ym; j++) {
-     row.j = j;
-     for (i=info->xs; i<info->xs+info->xm; i++) {
-         v[0] = ...;                                           col[0].j = j - 1;    col[0].i = i;
-         v[1] = ...;                                           col[1].j = j;         col[1].i = i-1;
-         v[2] = ...;                                           col[2].j = j;         col[2].i = i;
-         v[3] = ...;                                           col[3].j = j;         col[3].i = i+1;
-         v[4] = ...;                                           col[4].j = j + 1;   col[4].i = i;
-         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
-      }
-    }

- Multigrid Structured Mesh Support:
- Sample Jacobian Computation

## DMMG   dmmg;
DMMGCreate(comm,nlevels,user,&dmmg)
DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,4,                 	        4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da)
DMMGSetDM(dmmg,da)  
DMMGSetSNESLocal(dmmg,Function,Jacobian,0,0)
DMMGSolve(dmmg)
solution = DMMGGetx(damg)

- DMMGCreate(comm,nlevels,user,&dmmg)
- DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,4,                 	        4,PETSC_DECIDE,PETSC_DECIDE,4,1,0,0,&da)
- DMMGSetDM(dmmg,da)
- DMMGSetSNESLocal(dmmg,Function,Jacobian,0,0)
- DMMGSolve(dmmg)
- solution = DMMGGetx(damg)

-  2-dim nonlinear problem with 4 degrees of freedom per mesh point
-  Function() and Jacobian() are user-provided functions

- All standard SNES, SLES, PC and MG options apply.

- Multigrid Structured Mesh Support:
- Nonlinear Example

## Collaboration with P. Hovland and B. Norris (see http://www.mcs.anl.gov/autodiff)
Additional alternatives 
Compute sparse Jacobian explicitly using AD
DMMGSetSNESLocal(dmmg,Function,0,ad_Function,0)
PETSc + ADIC automatically generate ad_Function

Provide a “matrix-free” application of the Jacobian using AD
DMMGSetSNESLocal(dmmg,Function, 0,0,admf_Function)
PETSc + ADIC automatically generate admf_Function
Similar situation for Fortran and ADIFOR

- Additional alternatives
- Compute sparse Jacobian explicitly using AD
- DMMGSetSNESLocal(dmmg,Function,0,ad_Function,0)
- PETSc + ADIC automatically generate ad_Function
- Provide a “matrix-free” application of the Jacobian using AD
- DMMGSetSNESLocal(dmmg,Function, 0,0,admf_Function)
- PETSc + ADIC automatically generate admf_Function
- Similar situation for Fortran and ADIFOR

- Multigrid Structured Mesh Support:
- Jacobian via Automatic Differentiation

## Interface Approach
Based on interfacing at the matrix level, where external linear solvers typically use a variant of compressed sparse row matrix storage
Usage
Install PETSc indicating presence of any optional external packages in the file petsc/bmake/$PETSC_ARCH/base.site, e.g.,
PETSC_HAVE_SPAI = -DPETSC_HAVE_SPAI
SPAI_INCLUDE = -I/home/username/software/spai_3.0/include
SPAI_LIB = /home/username/software/spai_3.0/lib/${PETSC_ARCH}/libspai.a
Set preconditioners via the usual approach
Procedural interface: PCSetType(pc,”spai”)
Runtime option: -pc_type spai
Set preconditioner-specific options via the usual approach, e.g., 
PCSPAISetEpsilon(), PCSPAISetVerbose(), etc.
-pc_spai_epsilon <eps>   -pc_spai_verbose etc.

- Based on interfacing at the matrix level, where external linear solvers typically use a variant of compressed sparse row matrix storage
- Usage
- Install PETSc indicating presence of any optional external packages in the file petsc/bmake/$PETSC_ARCH/base.site, e.g.,
- PETSC_HAVE_SPAI = -DPETSC_HAVE_SPAI
- SPAI_INCLUDE = -I/home/username/software/spai_3.0/include
- SPAI_LIB = /home/username/software/spai_3.0/lib/${PETSC_ARCH}/libspai.a
- Set preconditioners via the usual approach
- Procedural interface: PCSetType(pc,”spai”)
- Runtime option: -pc_type spai
- Set preconditioner-specific options via the usual approach, e.g.,
- PCSPAISetEpsilon(), PCSPAISetVerbose(), etc.
- -pc_spai_epsilon <eps>   -pc_spai_verbose etc.

- Using PETSc with Other Packages:Linear Solvers

## AMG 
Algebraic multigrid code by J. Ruge, K. Steuben, and R. Hempel (GMD)
http://www.mgnet.org/mgnet-codes-gmd.html
PETSc interface by D. Lahaye (K.U.Leuven), uses MatSeqAIJ
BlockSolve95
Parallel, sparse ILU(0) for symmetric nonzero structure and ICC(0)
M. Jones (Virginia Tech.) and P. Plassmann (Penn State Univ.)
http://www.mcs.anl.gov/BlockSolve95
PETSc interface uses MatMPIRowbs
ILUTP
Drop tolerance ILU by Y. Saad (Univ. of Minnesota), in SPARSKIT
http://www.cs.umn.edu/~saad/
PETSc interface uses MatSeqAIJ 

- Algebraic multigrid code by J. Ruge, K. Steuben, and R. Hempel (GMD)
- _[http://www.mgnet.org/mgnet-codes-gmd.html](http://www.mgnet.org/mgnet-codes-gmd.html)_
- PETSc interface by D. Lahaye (K.U.Leuven), uses MatSeqAIJ
- BlockSolve95
- Parallel, sparse ILU(0) for symmetric nonzero structure and ICC(0)
- M. Jones (Virginia Tech.) and P. Plassmann (Penn State Univ.)
- _[http://www.mcs.anl.gov/BlockSolve95](http://www.mcs.anl.gov/BlockSolve95)_
- PETSc interface uses MatMPIRowbs
- ILUTP
- Drop tolerance ILU by Y. Saad (Univ. of Minnesota), in SPARSKIT
- _[http://www.](http://www.cs.umn.edu/~saad/)__[cs](http://www.cs.umn.edu/~saad/)__[.](http://www.cs.umn.edu/~saad/)__[umn](http://www.cs.umn.edu/~saad/)__[.](http://www.cs.umn.edu/~saad/)__[edu](http://www.cs.umn.edu/~saad/)__[/~](http://www.cs.umn.edu/~saad/)__[saad](http://www.cs.umn.edu/~saad/)__[/](http://www.cs.umn.edu/~saad/)_
- PETSc interface uses MatSeqAIJ

- Using PETSc with Other Packages:Linear Solvers

## LUSOL
Sparse LU, part of MINOS
M. Saunders (Stanford Univ)
http://www.sbsi-sol-optimize.com
PETSc interface by T. Munson (ANL), uses MatSeqAIJ
SPAI
Sparse approximate inverse code by S. Barnhard (NASA Ames) and M. Grote (ETH Zurich)
http://www.sam.math.ethz.ch/~grote/spai
PETSc interface converts from any matrix format to SPAI matrix
SuperLU
Parallel, sparse LU
J. Demmel, J. Gilbert, (U.C. Berkeley) and X. Li (NERSC)
http://www.nersc.gov/~xiaoye/SuperLU
PETSc interface uses MatSeqAIJ
Currently only sequential interface supported; parallel interface under development

- Sparse LU, part of MINOS
- M. Saunders (Stanford Univ)
- _[http://www.](http://www.sbsi-sol-optimize.com/)__[sbsi](http://www.sbsi-sol-optimize.com/)__[-sol-optimize.com](http://www.sbsi-sol-optimize.com/)_
- PETSc interface by T. Munson (ANL), uses MatSeqAIJ
- SPAI
- Sparse approximate inverse code by S. Barnhard (NASA Ames) and M. Grote (ETH Zurich)
- _[http://www.sam.math.ethz.ch/~grote/spai](http://www.sam.math.ethz.ch/~grote/spai)_
- PETSc interface converts from any matrix format to SPAI matrix
- SuperLU
- Parallel, sparse LU
- J. Demmel, J. Gilbert, (U.C. Berkeley) and X. Li (NERSC)
- _[http://www.nersc.gov/~xiaoye/](http://www.nersc.gov/~xiaoye/SuperLU)__[SuperLU](http://www.nersc.gov/~xiaoye/SuperLU)_
- PETSc interface uses MatSeqAIJ
- Currently only sequential interface supported; parallel interface under development

- Using PETSc with Other Packages:Linear Solvers (cont.)

## TAO - Toolkit for Advanced Optimization 
Software for large-scale optimization problems 
S. Benson, L. McInnes, and J. Moré
http://www.mcs.anl.gov/tao
Initial TAO design uses PETSc for 
Low-level system infrastructure - managing portability
Parallel linear algebra tools (SLES)
Veltisto (library for PDE-constrained optimization by G. Biros, see http://www.cs.nyu.edu/~biros/veltisto) – uses a similar interface approach
TAO is evolving toward
CCA-compliant component-based design (see http://www.cca-forum.org)
Support for ESI interfaces to various linear algebra libraries (see http://z.ca.sandia.gov/esi)

- Software for large-scale optimization problems
- S. Benson, L. McInnes, and J. Moré
- _[http://www.mcs.anl.gov/](http://www.mcs.anl.gov/tao)__[tao](http://www.mcs.anl.gov/tao)_
- Initial TAO design uses PETSc for
- Low-level system infrastructure - managing portability
- Parallel linear algebra tools (SLES)
- Veltisto (library for PDE-constrained optimization by G. Biros, see _[http://www.cs.nyu.edu/~biros/](http://www.cs.nyu.edu/~biros/veltisto)__[veltisto](http://www.cs.nyu.edu/~biros/veltisto)_) – uses a similar interface approach
- TAO is evolving toward
- CCA-compliant component-based design (see _[http://www.cca-forum.org](http://www.cca-forum.org/)_)
- Support for ESI interfaces to various linear algebra libraries (see _[http://z.ca.sandia.gov/](http://z.ca.sandia.gov/esi)__[esi](http://z.ca.sandia.gov/esi)_)

- Using PETSc with Other Packages:TAO – Optimization Software

## TAO Interface

**TAO_SOLVER**     tao;                     /*  optimization solver  */
- Vec                      x, g;                      /*  solution and gradient vectors  */
- ApplicationCtx     usercontext;         /*  user-defined context  */
**TaoInitialize();**
- /* Initialize Application -- Create variable and gradient vectors x and g */   ...
**TaoCreate**(MPI_COMM_WORLD,”tao_lmvm”,&tao);
**TaoSetFunctionGradient**(tao,x,g, FctGrad,(void*)&usercontext);
**TaoSolve**(tao);
- /* Finalize application -- Destroy vectors x and g */   ...
**TaoDestroy**(tao);
**TaoFinalize**();

- Similar Fortran interface, e.g., call TaoCreate(...)

- software interfacing: TAO

## Using PETSc with Other Packages:PVODE – ODE Integrators

- PVODE
- Parallel, robust, variable-order stiff and non-stiff ODE integrators
- A. Hindmarsh et al. (LLNL)
- _[http://www.](http://www.llnl.gov/CASC/PVODE)__[llnl](http://www.llnl.gov/CASC/PVODE)__[.gov/CASC/PVODE](http://www.llnl.gov/CASC/PVODE)_
- L. Xu developed PVODE/PETSc interface
- Interface Approach
- PVODE
- ODE integrator – evolves field variables in time
- vector – holds field variables
- preconditioner placeholder
- Usage
- TSCreate(MPI_Comm,TS_NONLINEAR,&ts)
- TSSetType(ts,TS_PVODE)
- ..... regular TS functions
- TSPVODESetType(ts,PVODE_ADAMS)
- .... other PVODE options
- TSSetFromOptions(ts) – accepts PVODE options

-      PETSc
-     ODE integrator placeholder
-      vector
-      sparse matrix and preconditioner

## SUMAA3d
Scalable Unstructured Mesh Algorithms and Applications
L. Freitag (ANL), M. Jones (VA Tech), P. Plassmann (Penn State)
http://www.mcs.anl.gov/sumaa3d
L. Freitag and M. Jones developed SUMAA3d/PETSc interface
SAMRAI
Structured adaptive mesh refinement
R. Hornung, S. Kohn (LLNL) 
http://www.llnl.gov/CASC/SAMRAI
SAMRAI team developed SAMRAI/PETSc interface
Overture
Structured composite meshes and discretizations
D. Brown, W. Henshaw, D. Quinlan (LLNL)
http://www.llnl.gov/CASC/Overture
K. Buschelman and Overture team developed Overture/PETSc interfaces

- Scalable Unstructured Mesh Algorithms and Applications
- L. Freitag (ANL), M. Jones (VA Tech), P. Plassmann (Penn State)
- _[http://www.mcs.anl.gov/sumaa3d](http://www.mcs.anl.gov/sumaa3d)_
- L. Freitag and M. Jones developed SUMAA3d/PETSc interface
- SAMRAI
- Structured adaptive mesh refinement
- R. Hornung, S. Kohn (LLNL)
- _[http://www.llnl.gov/CASC/SAMRAI](http://www.llnl.gov/CASC/SAMRAI)_
- SAMRAI team developed SAMRAI/PETSc interface
- Overture
- Structured composite meshes and discretizations
- D. Brown, W. Henshaw, D. Quinlan (LLNL)
- _[http://www.](http://www.llnl.gov/CASC/Overture)__[llnl](http://www.llnl.gov/CASC/Overture)__[.gov/CASC/Overture](http://www.llnl.gov/CASC/Overture)_
- K. Buschelman and Overture team developed Overture/PETSc interfaces

- Using PETSc with Other Packages:Mesh Management and Discretization

## Using PETSc with Other Packages:Matlab

- Matlab
- _[http://www.mathworks.com](http://www.mathworks.com/)_
- Interface Approach
- PETSc socket interface to Matlab
- Sends matrices and vectors to interactive Matlab session
- PETSc interface to MatlabEngine
- MatlabEngine – Matlab library that allows C/Fortran programmers to use Matlab functions in programs
- PetscMatlabEngine – unwraps PETSc vectors and matrices so that MatlabEngine can understand them
- Usage
- PetscMatlabEngineCreate(MPI_Comm,machinename,  	PetscMatlabEngine eng)
- PetscMatlabEnginePut(eng,PetscObject obj)
- Vector
- Matrix
- PetscMatlabEngineEvaluate(eng,”R = QR(A);”)
- PetscMatlabEngineGet(eng,PetscObject obj)

## Using PETSc with Other Packages:ParMETIS – Graph Partitioning

- ParMETIS
- Parallel graph partitioning
- G. Karypis (Univ. of Minnesota)
- _[http://www.](http://www.cs.umn.edu/~karypis/metis/parmetis)__[cs](http://www.cs.umn.edu/~karypis/metis/parmetis)__[.](http://www.cs.umn.edu/~karypis/metis/parmetis)__[umn](http://www.cs.umn.edu/~karypis/metis/parmetis)__[.](http://www.cs.umn.edu/~karypis/metis/parmetis)__[edu](http://www.cs.umn.edu/~karypis/metis/parmetis)__[/~](http://www.cs.umn.edu/~karypis/metis/parmetis)__[karypis](http://www.cs.umn.edu/~karypis/metis/parmetis)__[/](http://www.cs.umn.edu/~karypis/metis/parmetis)__[metis](http://www.cs.umn.edu/~karypis/metis/parmetis)__[/](http://www.cs.umn.edu/~karypis/metis/parmetis)__[parmetis](http://www.cs.umn.edu/~karypis/metis/parmetis)_
- Interface Approach
- Use PETSc MatPartitioning() interface and MPIAIJ or MPIAdj matrix formats
- Usage
- MatPartitioningCreate(MPI_Comm,MatPartitioning ctx)
- MatPartitioningSetAdjacency(ctx,matrix)
- Optional – MatPartitioningSetVertexWeights(ctx,weights)
- MatPartitioningSetFromOptions(ctx)
- MatPartitioningApply(ctx,IS *partitioning)

## Recent Additions

- Hypre (www.llnl.gov/casc/hypre) via PCHYPRE, includes
- EUCLID (parallel ILU(k) by David Hysom)
- BoomerAMG
- PILUT
- DSCPACK, a Domain-Separator Cholesky Package for solving sparse spd problems, by Padma Raghavan
- SPOOLES (SParse Object Oriented Linear Equations Solver), by Cleve Ashcraft
- Both SuperLU and SuperLU_Dist sequential and parallel direct sparse solvers, by Xiaoye Li et al.

## Extensibility Issues

- 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

- 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 and performance issues by writing to petsc-maint@mcs.anl.gov.

## References

- Documentation: _[http://www.mcs.anl.gov/petsc/docs](http://www.mcs.anl.gov/petsc/docs)_
- PETSc Users manual
- Manual pages
- Many hyperlinked examples
- FAQ, Troubleshooting info, installation info, etc.
- Publications: _[http://www.mcs.anl.gov/petsc/publications](http://www.mcs.anl.gov/petsc/publications)_
- Research and publications that make use PETSc
- MPI Information: _[http://www.mpi-forum.org](http://www.mpi-forum.org/)_
***Using MPI***** **(2nd Edition), by Gropp, Lusk, and Skjellum
***Domain Decomposition***, by Smith, Bjorstad, and Gropp