PETSc tools

Experimental html version of downloadable textbook, see http://www.tacc.utexas.edu/~eijkhout/istc/istc.html
\[ \newcommand\inv{^{-1}}\newcommand\invt{^{-t}} \newcommand\bbP{\mathbb{P}} \newcommand\bbR{\mathbb{R}} \newcommand\defined{ \mathrel{\lower 5pt \hbox{${\equiv\atop\mathrm{\scriptstyle D}}$}}} \] 38.1 : Error checking and debugging
38.1.1 : Debug mode
38.1.2 : Error codes
38.1.3 : Memory corruption
38.1.3.1 : Valgrind
38.2 : Program output
38.2.1 : Screen I/O
38.2.1.1 : printf replacements
38.2.1.2 : scanf replacement
38.2.2 : Exporting internal data structures
38.2.2.1 : Viewer types
38.2.2.2 : Viewer formats
38.2.2.3 : Commandline option for viewers
38.2.2.4 : Naming objects
38.3 : Commandline options
38.3.1 : Adding your own options
38.3.2 : Options prefix
38.3.3 : Where to specify options
38.4 : Timing and profiling
38.4.1 : Logging
38.5 : Memory management
38.5.1 : GPU allocation
Back to Table of Contents

38 PETSc tools

38.1 Error checking and debugging

crumb trail: > petsc-tools > Error checking and debugging

38.1.1 Debug mode

crumb trail: > petsc-tools > Error checking and debugging > Debug mode

During installation (see section~ 31.3 ), there is an option of turning on debug mode. An installation with debug turned on:

  • Does more runtime checks on numerics, or array indices;
  • Does a memory analysis when you insert the CHKMEMQ macro (section~ 38.1.3 );
  • Has the macro PETSC_USE_DEBUG set to~1.

38.1.2 Error codes

crumb trail: > petsc-tools > Error checking and debugging > Error codes

PETSc performs a good amount of runtime error checking. Some of this is for internal consistency, but it can also detect certain mathematical errors. To facilitate error reporting, the following scheme is used.

  1. Every PETSc routine is a function returning a parameter of type PetscErrorCode .
  2. For a good traceback, surround the executable part of any subprogram with PetscFunctionBegin and PetscFunctionReturn , where the latter has the return value as parameter.
  3. Calling the macro CHKERRQ on the error code will cause an error to be printed and the current routine to be terminated. Recursively this gives a traceback of where the error occurred.

    PetscErrorCode ierr;
    ierr = AnyPetscRoutine( arguments ); CHKERRQ(ierr);
    

  4. You can effect your own error return by using
    #include 
    PetscErrorCode SETERRQ (MPI_Comm comm,PetscErrorCode ierr,char *message)
    PetscErrorCode SETERRQ1(MPI_Comm comm,PetscErrorCode ierr,char *formatmessage,arg1)
    PetscErrorCode SETERRQ2(MPI_Comm comm,PetscErrorCode ierr,char *formatmessage,arg1,arg2)
    PetscErrorCode SETERRQ3(MPI_Comm comm,PetscErrorCode ierr,char *formatmessage,arg1,arg2,arg3)
    
    Input Parameters:
    comm - A communicator, so that the error can be collective
    ierr - nonzero error code, see the list of standard error codes in include/petscerror.h
    message - error message in the printf format
    arg1,arg2,arg3 - argument (for example an integer, string or double)
    
    
    SETERRQ

Fortran note

In the main program, use CHKERRA and SETERRA . Also beware that these error `commands' are macros, and after expansion may interfere with Fortran line length , so they should only be used in .F90 files.

Example. We write a routine that sets an error:

// backtrace.c
PetscErrorCode this_function_bombs() {
  PetscFunctionBegin;
  SETERRQ(PETSC_COMM_SELF,1,"We cannot go on like this");
  PetscFunctionReturn(0);
}

Running this gives, in process zero, the output

[0]PETSC ERROR: We cannot go on like this
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.12.2, Nov, 22, 2019
[0]PETSC ERROR: backtrace on a [computer name]
[0]PETSC ERROR: Configure options [all options]
[0]PETSC ERROR: #1 this_function_bombs() line 20 in backtrace.c
[0]PETSC ERROR: #2 main() line 30 in backtrace.c

Fortran note

In Fortran the backtrace is not quite as elegant.

// backtrace.F90
Subroutine this_function_bombs(ierr)
  implicit none
  integer,intent(out) :: ierr

  SETERRQ(PETSC_COMM_SELF,1,"We cannot go on like this")
  ierr = -1

end Subroutine this_function_bombs

[0]PETSC ERROR: ----- Error Message ------------------------------
[0]PETSC ERROR: We cannot go on like this
[....]
[0]PETSC ERROR: #1 User provided function() line 0 in User file

Remark

In this example, the use of PETSC_COMM_SELF indicates that this error is individually generated on a process; use PETSC_COMM_WORLD only if the same error would be detected everywhere.

Exercise

Look up the definition of SETERRQ1 . Write a routine to compute square roots that is used as follows:

  x = 1.5; ierr = square_root(x,&rootx); CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"Root of %f is %f\n",x,rootx);
  x = -2.6; ierr = square_root(x,&rootx); CHKERRQ(ierr);
  PetscPrintf(PETSC_COMM_WORLD,"Root of %f is %f\n",x,rootx);

This should give as output:

Root of 1.500000 is 1.224745
[0]PETSC ERROR: ----- Error Message ----------------------------------------------
[0]PETSC ERROR: Cannot compute the root of -2.600000
[...]
[0]PETSC ERROR: #1 square_root() line 23 in root.c
[0]PETSC ERROR: #2 main() line 39 in root.c

38.1.3 Memory corruption

crumb trail: > petsc-tools > Error checking and debugging > Memory corruption

PETSc has its own memory management (section  38.5 ) and this facilitates finding memory corruption errors. The macro CHKMEMQ ( CHKMEMA in void functions) checks all memory that was allocated by PETSc, either internally or throug the allocation routines, for corruption. Sprinkling this macro through your code can detect memory problems before they lead to a segfault .

This testing is only done if the commandline argument -malloc_debug ( -malloc_test in debug mode) is supplied, so it carries no overhead for production runs.

38.1.3.1 Valgrind

crumb trail: > petsc-tools > Error checking and debugging > Memory corruption > Valgrind

Valgrind is rather verbose in its output. To limit the number of processs that run under valgrind:

mpiexec -n 3 valgrind --track-origins=yes ./app -args : -n 5 ./app -args

38.2 Program output

crumb trail: > petsc-tools > Program output

PETSc has as variety of mechanisms to export or visualize program data. We will consider a few possibilities here.

38.2.1 Screen I/O

crumb trail: > petsc-tools > Program output > Screen I/O

Printing screen output in parallel is tricky. If two processes execute a print statement at more or less the same time there is no guarantee as to in what order they may appear on screen. (Even attempts to have them print one after the other may not result in the right ordering.) Furthermore, lines from multi-line print actions on two processes may wind up on the screen interleaved.

38.2.1.1 printf replacements

crumb trail: > petsc-tools > Program output > Screen I/O > printf replacements

PETSc has two routines that fix this problem. First of all, often the information printed is the same on all processes, so it is enough if only one process, for instance process 0, prints it. This is done with

C:
PetscErrorCode  PetscPrintf(MPI_Comm comm,const char format[],...)

Fortran:
PetscPrintf(MPI_Comm, character(*), PetscErrorCode ierr)

Python:
PETSc.Sys.Print(type cls, *args, **kwargs)
kwargs:
comm : communicator object
PetscPrintf .

If all processes need to print, you can use

C:
PetscErrorCode  PetscSynchronizedPrintf(
    MPI_Comm comm,const char format[],...)

Fortran:
PetscSynchronizedPrintf(MPI_Comm, character(*), PetscErrorCode ierr)

python:
PETSc.Sys.syncPrint(type cls, *args, **kargs)
kwargs:
comm : communicator object
flush : if True, do synchronizedFlush
other keyword args as for python3 print function
PetscSynchronizedPrintf that forces the output to appear in process order.

To make sure that output is properly flushed from all system buffers use

C:
PetscErrorCode  PetscSynchronizedFlush(MPI_Comm comm,FILE *fd)
fd : output file pointer, needs to be valid on process zero

Fortran:
PetscSynchronizedFlush(comm,fd,err)
Integer :: comm
fd is usually PETSC_STDOUT
PetscErrorCode :: err

python:
PETSc.Sys.syncFlush(type cls, comm=None)
PetscSynchronizedFlush where for ordinary screen output you would use stdout for the file.

Fortran note

The Fortran calls are only wrappers around C routines, so you can use \n newline characters in the Fortran string argument to PetscPrintf .

The file to flush is typically PETSC_STDOUT .

Python note

Since the print routines use the python print call, they automatically include the trailing newline. You don't have to specify it as in the C calls.

38.2.1.2 scanf replacement

crumb trail: > petsc-tools > Program output > Screen I/O > scanf replacement

Using scanf in Petsc is tricky, since integers and real numbers can be of different sizes, depending on the installation. Instead, use

Synopsis

#include "petscviewer.h"
PetscErrorCode  PetscViewerRead(PetscViewer viewer, void *data, PetscInt num, PetscInt *count, PetscDataType dtype)

Collective

Input Parameters
	viewer 	- The viewer
	data 	- Location to write the data
	num 	- Number of items of data to read
	datatype 	- Type of data to read

Output Parameters
count -number of items of data actually read, or NULL
PetscViewerRead , which operates in terms of PetscDataType .

38.2.2 Exporting internal data structures

crumb trail: > petsc-tools > Program output > Exporting internal data structures

In order to export PETSc matrix or vector data structures there is a PetscViewer object type. This is a quite general concept of viewing: it encompasses ascii output to screen, binary dump to file, or communication to a running Matlab process. Calls such as MatView or KSPView accept a PetscViewer argument.

In cases where this makes sense, there is also an inverse `load' operation. See section 32.3.5 for vectors.

Some viewers are predefined, such as PETSC_VIEWER_STDOUT_WORLD for ascii rendering to standard out. (In C, specifying zero or NULL also uses this default viewer; for Fortran use

38.2.2.1 Viewer types

crumb trail: > petsc-tools > Program output > Exporting internal data structures > Viewer types

For activities such as dumping to file you first need create the viewer with PetscViewerCreate and set its type with PetscViewerSetType .

PetscViewerCreate(comm,&viewer);
PetscViewerSetType(viewer,PETSCVIEWERBINARY);

Popular types include PETSCVIEWERASCII , PETSCVIEWERBINARY , PETSCVIEWERSTRING , PETSCVIEWERDRAW , PETSCVIEWERSOCKET , PETSCVIEWERHDF5 , PETSCVIEWERVTK ; the full list can be found in include/petscviewer.h .

38.2.2.2 Viewer formats

crumb trail: > petsc-tools > Program output > Exporting internal data structures > Viewer formats

Viewers can take further format specifications by using PetscViewerPushFormat :

PetscViewerPushFormat
   (PETSC_VIEWER_STDOUT_WORLD,
    PETSC_VIEWER_ASCII_INFO_DETAIL);

and afterwards a corresponding PetscViewerPopFormat

38.2.2.3 Commandline option for viewers

crumb trail: > petsc-tools > Program output > Exporting internal data structures > Commandline option for viewers

Petsc objects viewers can be activated by calls such as MatView , but often it is more convenient to do this through commandline options, such as mat_view , vec_view , or ksp_view . By default, these output to stdout in ascii form, but this can be controlled by further option values:

program -mat_view binary:matrix.dat

where binary forces a binary dump ( ascii is the default) and a file name is explicitly given.

If a viewer needs to be triggered at a specific location, calls such as VecViewFromOptions can be used. VecSetOptionsPrefix . These routines all have a similar calling sequence:

#include "petscsys.h"
PetscErrorCode  PetscObjectViewFromOptions(PetscObject obj,PetscObject bobj,const char optionname[])
PetscErrorCode  VecViewFromOptions(Vec A,PetscObject obj,const char name[])

\begin{raggedlist} AOViewFromOptions , DMViewFromOptions , ISViewFromOptions , ISLocalToGlobalMappingViewFromOptions , KSPConvergedReasonViewFromOptions , KSPViewFromOptions , MatPartitioningViewFromOptions , MatCoarsenViewFromOptions , MatViewFromOptions , PetscObjectViewFromOptions , PetscPartitionerViewFromOptions , PetscDrawViewFromOptions , PetscRandomViewFromOptions , PetscDualSpaceViewFromOptions , PetscSFViewFromOptions , PetscFEViewFromOptions , PetscFVViewFromOptions , PetscSectionViewFromOptions , PCViewFromOptions , PetscSpaceViewFromOptions , PFViewFromOptions , PetscLimiterViewFromOptions , PetscLogViewFromOptions , PetscDSViewFromOptions , PetscViewerViewFromOptions , SNESConvergedReasonViewFromOptions , SNESViewFromOptions , TSTrajectoryViewFromOptions , TSViewFromOptions , TaoLineSearchViewFromOptions , TaoViewFromOptions , VecViewFromOptions , VecScatterViewFromOptions , \end{raggedlist}

38.2.2.4 Naming objects

crumb trail: > petsc-tools > Program output > Exporting internal data structures > Naming objects

A helpful facility for viewing is to name an object: that name will then be displayed when the object is viewed.

Vec i_local;
ierr = VecCreate(comm,&i_local); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)i_local,"space local"); CHKERRQ(ierr);

giving:

Vec Object: space local 4 MPI processes
  type: mpi
Process [0]
[ ... et cetera ... ]

38.3 Commandline options

crumb trail: > petsc-tools > Commandline options

PETSc has as large number of commandline options, most of which we will discuss later. For now we only mention -log_summary which will print out profile of the time taken in various routines. For these options to be parsed, it is necessary to pass argc,argv to the PetscInitialize call.

38.3.1 Adding your own options

crumb trail: > petsc-tools > Commandline options > Adding your own options

You can add custom commandline options to your program. Various routines such as PetscOptionsGetInt scan the commandline for options and set parameters accordingly. For instance,

// ksp.c
PetscBool flag;
int domain_size = 100;
ierr = PetscOptionsGetInt
  (NULL,PETSC_NULL,"-n",&domain_size,&flag); CHKERRQ(ierr);
PetscPrintf(comm,"Using domain size %d\n",domain_size);
declares the existence of an option -n to be followed by an integer.

Now executing

mpiexec yourprogram -n 5

will

  1. set the flag to true, and
  2. set the parameter domain_size to the value on the commandline.

Omitting the -n option will leave the default value of domain_size unaltered.

For flags, use PetscOptionsHasName .

Python note

In Python, do not specify the initial hyphen of an option name. Also, the functions such as getInt do not return the boolean flag; if you need to test for the existence of the commandline option, use:

hasn = PETSc.Options().hasName("n")

There is a related mechanism using PetscOptionsBegin  / PetscOptionsEnd :

// optionsbegin.c
ierr = PetscOptionsBegin(comm,NULL,"Parameters",NULL); CHKERRQ(ierr);
ierr = PetscOptionsInt("-i","i value",__FILE__,i_value,&i_value,&i_flag); CHKERRQ(ierr);
ierr = PetscOptionsInt("-j","j value",__FILE__,j_value,&j_value,&j_flag); CHKERRQ(ierr);
ierr = PetscOptionsEnd(); CHKERRQ(ierr);
if (i_flag)
  PetscPrintf(comm,"Option `-i' was used\n");
if (j_flag)
  PetscPrintf(comm,"Option `-j' was used\n");

The selling point for this approach is that running your code with

mpiexec yourprogram -help

will display these options as a block. Together with a ton of other options, unfortunately.

38.3.2 Options prefix

crumb trail: > petsc-tools > Commandline options > Options prefix

In many cases, your code will have only one KSP solver object, so specifying ksp_view or ksp_monitor will display / trace that one. However, you may have multiple solvers, or nested solvers. You may then not want to display all of them.

As an example of the nest solver case, consider the case of a block jacobi preconditioner , where the block is itself solved with an iterative method. You can trace that one with -sub_ksp_monitor .

The sub_ is an option prefix , and you can defined your own with KSPSetOptionsPrefix . (There are similar routines for other PETSc object types.)

Example:

KSPCreate(comm,&time_solver);
KSPCreate(comm,&space_solver);
KSPSetOptionsPrefix(time_solver,"time_");
KSPSetOptionsPrefix(space_solver,"space_");

You can then use options -time_ksp_monitor and such. Note that the prefix does not have a leading dash, but it does have the trailing underscore.

\begin{raggedlist} Similar routines: MatSetOptionsPrefix , PCSetOptionsPrefix , PetscObjectSetOptionsPrefix , PetscViewerSetOptionsPrefix , SNESSetOptionsPrefix , TSSetOptionsPrefix , VecSetOptionsPrefix , and some more obscure ones. \end{raggedlist}

38.3.3 Where to specify options

crumb trail: > petsc-tools > Commandline options > Where to specify options

Commandline options can obviously go on the commandline. However, there are more places where they can be specified.

Options can be specified programmatically with PetscOptionsSetValue :

PetscOptionsSetValue( NULL, // for global options
  "-some_option","value_as_string");

Options can be specified in a file .petscrc in the user's home directory or the current directory.

Finally, an environment variable PETSC_OPTIONS can be set.

The rc file is processed first, then the environment variable, then any commandline arguments. This parsing is done in PetscInitialize , so any values from PetscOptionsSetValue override this.

38.4 Timing and profiling

crumb trail: > petsc-tools > Timing and profiling

PETSc has a number of timing routines that make it unnecessary to use system routines such as getrusage or MPI routines such as MPI_Wtime . The main (wall clock) timer is

Synopsis
Returns the CPU time in seconds used by the process.

#include "petscsys.h"
#include "petsctime.h"
PetscErrorCode PetscGetCPUTime(PetscLogDouble *t)
PetscErrorCode PetscTime(PetscLogDouble *v)
PetscTime . Note the return type of PetscLogDouble which can have a different precision from PetscReal .

The routine PetscGetCPUTime is less useful, since it measures only time spent in computation, and ignores things such as communication.

38.4.1 Logging

crumb trail: > petsc-tools > Timing and profiling > Logging

Petsc does a lot of logging on its own operations. Additionally, you can introduce your own routines into this log.

The simplest way to display statistics is to run with an option log_view . This takes an optional file name argument:

mpiexec -n 10 yourprogram -log_view :statistics.txt

The corresponding routine is PetscLogView .

38.5 Memory management

crumb trail: > petsc-tools > Memory management

Allocate the memory for a given pointer: PetscNew , allocate arbitrary memory with PetscMalloc , allocate a number of objects with

Synopsis
Allocates an array of memory aligned to PETSC_MEMALIGN

C:
#include 
PetscErrorCode PetscMalloc1(size_t m1,type **r1)

Input Parameter:
m1 - number of elements to allocate (may be zero)

Output Parameter:
r1 - memory allocated
PetscMalloc1 (this does not zero the memory allocated, use PetscCalloc1 to obtain memory that has been zeroed); use
Synopsis
Frees memory, not collective

C:
#include 
PetscErrorCode PetscFree(void *memory)

Input Parameter:
memory - memory to free (the pointer is ALWAYS set to NULL upon sucess)
PetscFree to free.

PetscInt *idxs;
PetscMalloc1(10,&idxs);
// better than:
// PetscMalloc(10*sizeof(PetscInt),&idxs);
for (PetscInt i=0; i<10; i++)
  idxs[i] = f(i);
PetscFree(idxs);

Allocated memory is aligned to PETSC_MEMALIGN .

The state of memory allocation can be written to file or standard out with PetscMallocDump . The commandline option -malloc_dump outputs all not-freed memory during PetscFinalize .

38.5.1 GPU allocation

crumb trail: > petsc-tools > Memory management > GPU allocation

The memories of a CPU and GPU are not coherent. This means that routines such as PetscMalloc1 can not immediately be used for GPU allocation. See section  37.4 for details.

Back to Table of Contents