Example program structure

The following is an example of the global structure of a program using ParPre preconditioners in an iterative method. Many details have been omitted here.

First of all, the main program declares a coefficient matrix and a preconditioner (see section 2.2):

#include "petsc.h"
#include "parpre.h"
int main(int Argc,char **Argv)
{
  Mat       A;
  PC        B;

Since the ParPre library is based on Petsc, we have to initialise Petsc, and create the data structures for the matrix and the preconditioner:

  PetscInitialize(&Argc,&Args,PETSC_NULL,PETSC_NULL);
  ierr = MatCreateMPIAIJ
         (MPI_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,
          xdim*ydim,xdim*ydim,0,0,0,0,&A); 
  ierr = PCCreate(MPI_COMM_WORLD,&the_pc);

The three main user routines are now to construct the matrix and preconditioner, and to use them in an iterative method:

  ierr = make_mat(A, [other arguments] );
  ierr = prec_setup(A,&B);
  ierr = cg(A,B, [other arguments] );

At the end of the program, the data structures have to be released again:

  MatDestroy(A);
  PCDestroy(B);
  PetscFinalize();
}

The matrix construction routine calculates or retrieves the elements of the coefficient matrix, and stores them in a Petsc version of the matrix. This doubling of the matrix storage is unfortunately unavoidable, but see section 3.3 for details.

int make_mat(Mat A, [other arguments] )
{
  for ( i=0; i<m; i++ ) { 
    for ( j=jlo; j<jhi; j++ ) {
      /* compute or retrieve element value */
      ierr = MatSetValues(A,1,&I,1,&I,&val,INSERT_VALUES);
    }
  }

In the conjugate gradient method, vectors are allocated normally as double*; the Petsc vectors subsequently created from them are basically wrappers for these pointers; see section 3.2.

int cg(Mat A, PC B, [other arguments] )
{
  double *z,*r; /* allocatable vectors */
  Vec r_vec,z_vec; /* Petsc equivalents */

  ierr = VecCreateMPIWithArray(MPI_COMM_WORLD,lsize,r,&r_vec);
  ierr = VecCreateMPIWithArray(MPI_COMM_WORLD,lsize,z,&z_vec);

Application of the preconditioner is a simple call in the iterative loop (see section 2.3):

  for (it=1; rho/rho_0>tolerance; it++) {
    /* solve Mz = r */
    ierr = PCApply(B,r_vec,z_vec);
    /* other calls in the iterative method */
  }
}

Creating the preconditioner involves several standard calls, and possibly a number of method-specific auxiliaries; see section 2.2. The following is an example setup for the domain decomposition method:

int prec_setup(
  ierr = PCSetType(the_pc,PCDomainDecomp);
  {
    PC intface_pc;
    ierr = PCDomainDecompGetInterfacePC(the_pc,&intface_pc);
    ierr = PCSetType(intface_pc,PCJACOBI);
  }
  /* Declare local/global aspects */
  {
    PC local_pc;
    ierr = PCParallelGetLocalPC(the_pc,&local_pc);
    ierr = PCSetType(local_pc,PCLU);
  }
  ierr = ParPreSetup(MPI_COMM_WORLD,A,the_pc);
  return 0;
}

After this example, we will give detailed information on how to declare and construct preconditioners in the rest of this report.