PETSc (Portable Extensible Toolkit for Scientific computation) est une bibliothèque de fonctions mathématiques pour la résolution des équations aux dérivées partielles. Très utilisée en calcul parallèle avec MPI (Message Passing Interface). Elle est principalement réservée aux langages C, C++ et Fortran.

PETSc is intended for use in large-scale application projects, many ongoing computational science projects are built around the PETSc libraries. PETSc is easy to use for beginners. Moreover, its careful design allows advanced users to have detailed control over the solution process. PETSc includes a large suite of parallel linear, nonlinear equation solvers and ODE integrators that are easily used in application codes written in C, C++, Fortran and now Python. PETSc provides many of the mechanisms needed within parallel application codes, such as simple parallel matrix and vector assembly routines that allow the overlap of communication and computation. In addition, PETSc includes support for parallel distributed arrays useful for finite difference methods.

Features include:

• Parallel vectors
• includes code for communicating ghost points
• Parallel matrices
• several sparse storage formats
• easy, efficient assembly
• Scalable parallel preconditioners
• Krylov subspace methods
• Parallel Newton-based nonlinear solvers
• Parallel timestepping (ODE) solvers
• Support for Nvidia GPU cards
• Complete documentation
• Automatic profiling of floating point and memory usage
• Consistent user interface
• Intensive error checking
• Portable to UNIX and Windows
• Over one hundred examples
• PETSc is supported and will be actively enhanced for many years.

PETSc dépend de MPI (Open MPI), pour compiler des programmes il faut charger le module: numlib/petsc

$ module load  numlib/petsc

Exemple : HelloWolrd en PETSc

hello_petsc.c

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

 
mpicc hello_petsc.c -o hello_petsc \
 -I$PETSC_HOME/include \
 -L$PETSC_HOME/lib \
 -lpetsc \
 

Makefile

ALL: hello
 
include ${PETSC_DIR}/conf/variables
 
include ${PETSC_DIR}/conf/rules
 
 
hello: hello.o 
 
     ${CLINKER} -o hello hello.o ${PETSC_LIB}

$ module load  numlib/petsc && make

Exemple : exécution sur sur 4 processeurs

$ mpirun -np 4 ./hello_petsc
 
Hello World from rank 0
Hello World from rank 1
Hello World from rank 2
Hello World from rank 3

Généralement, ce type de programme ressemble à :

...
static char help []  = "blabla\n";
...
int main(int argc,char **args)
{
	PetscInitialize(&argc,&args,(char *)0,help);
	...
	PetscFinalize();
	return 0;
 }

Un objet de type PetscErrorCode est le plus souvent créé pour localiser directement les erreurs.

PetscErrorCode ierr;
...
ierr = PetscFinalize();CHKERRQ(ierr)

Ils s'utilisent avec petscvec.h, objet de type Vec. Trois commandes sont nécessaires à la création d'un vecteur, et une pour la destruction. Par exemple pour créer un vecteur b de dimension dim, mettre la valeur val dans sa première coordonnée, le visionner, et le supprimer, il faut taper :

...	
Vec b;
...
VecCreate(PETSC_COMM_WORLD,&b); 
VecSetSizes(b,PETSC_DECIDE,dim); 
VecSetFromOptions(b);
...
VecSetValue(b,0,val,INSERT_VALUES); //les tableaux commencent à l'indice 0
VecView(b,PETSC_VIEWER_STDOUT_SELF);
...
VecDestroy(b);
...

Vec.c

#include "petscvec.h"
 
int main(int argc,char **argv)
{
    Vec x;
    int n = 20,ierr;
 
    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(x,one);
 
    VecDot(x,x,&dot);
 
    PetscPrintf(PETSC_COMM_WORLD,"Vector length % d\n",(int)dot);
 
    VecDestroy(&x);
 
    PetscFinalize();
 
    return 0;
}

Ils sont accessible via petscmat.h, objet de type Mat. De même, trois commandes sont nécessaires à la création d'une matrice, et une pour la destruction. Par exemple pour créer une matrice A de dimension dim1*dim2, mettre la valeur val dans sa première ligne et deuxième colonne, la visionner, et la supprimer, il faut taper :

...	
Mat A;
...
MatCreate(PETSC_COMM_WORLD,&A); 
MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim1,dim2); 
MatSetFromOptions(A);
...
MatSetValue(A,0,1,val,INSERT_VALUES);
MatView(A,PETSC_VIEWER_STDOUT_SELF);
...
MatDestroy(A);
...

Type de Matrices

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.

Pour résoudre un système Ax=b avec A une matrice carrée, b et x deux vecteurs colonnes, il est recommandé d'utiliser les commandes KSP :

...	
KSP ksp; 	 
...
KSPCreate(PETSC_COMM_WORLD,&ksp); 
KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN); 
//-- A est stockée dans KSP
KSPSetTolerances(ksp,1.e-14,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 
//-- la tolérance souhaitée
KSPSetFromOptions(ksp); 
KSPSolve(ksp,b,x);  //-- le résultat est dans le vecteur x : Ax=b

Petsc4py est l'interface Python pour PETSc. Pour l'utiliser il faut charger python :

$ module load lang/python/2.7.11

Exemple d'utilisation

petsc.py

import petsc4py, sys
petsc4py.init(sys.argv)
 
from petsc4py import PETSc
 
from mpi4py import MPI
comm = MPI.COMM_WORLD
 
# grid size and spacing
m, n  = 101, 101
hx = 1.0/(m-1)
hy = 1.0/(n-1)
 
# create sparse matrix
A = PETSc.Mat()
A.create(PETSc.COMM_WORLD)
A.setSizes([m*n, m*n])
A.setType('aij') # sparse
A.setPreallocationNNZ(5)
 
# precompute values for setting
# diagonal and non-diagonal entries
diagv = 2.0/hx**2 + 2.0/hy**2
offdx = -1.0/hx**2
offdy = -1.0/hy**2
 
# loop over owned block of rows on this
# processor and insert entry values
Istart, Iend = A.getOwnershipRange()
for I in range(Istart, Iend) :
    A[I,I] = diagv
    i = I//n    # map row number to
    j = I - i*n # grid coordinates
    if i> 0  : J = I-n; A[I,J] = offdx
    if i< m-1: J = I+n; A[I,J] = offdx
    if j> 0  : J = I-1; A[I,J] = offdy
    if j< n-1: J = I+1; A[I,J] = offdy
# communicate off-processor values
# and setup internal data structures
# for performing parallel operations
A.assemble()
 
ksp = PETSc.KSP()
ksp.create(PETSc.COMM_WORLD)
# use conjugate gradients
ksp.setType('cg')
# and incomplete Cholesky
#ksp.getPC().setType('icc')
# obtain sol & rhs vectors
x, b = A.createVecs()
x.set(0)
b.set(1)
# and next solve
ksp.setOperators(A)
ksp.setFromOptions()
ksp.solve(b, x)
 
myviewer = PETSc.Viewer().createASCII("test.txt", format=PETSc.Viewer.Format.ASCII_VTK,comm= PETSc.COMM_WORLD)
 
x.view(myviewer)
comm.Barrier()

Exécution :

$ mpirun -np 4 python petsc.py

• Un très bon tutoriel sur l'utilisation PETSc pour la résolution des PDEs petscpde.pdf
• Tutoriel PETSc (en local) /Softs/libraries/petsc/3.1-p8/tutorials
• Le site officiel http://www.mcs.anl.gov/petsc/petsc-as/