Skip to content
This repository has been archived by the owner on Nov 27, 2024. It is now read-only.

[spec] Solver interface implementation

Jump to bottom
kynan edited this page Feb 23, 2013 · 1 revision

Note: this blueprint is out of date, the Solver interface has since been implemented using petsc4py

As discussed in Issue #68 we should provide an interface to solvers that meets the following constraints:

  1. The interface should allow the use of different solver backends
  2. There must be a way to pass options to the solver.
  3. The interface should be exposed to Fluidity users using a similar interface to that provided by Dolfin

Planned user interface

The user instantiates a PyOP2 solver object as:

solver = Solver()

The Solver object instantiates an instance of the backend library solver in its constructor. For now we only have PETSc, but we should ensure that it is easy to add other linear algebra backends. It is expected that the solver library may either be chosen in the OP2 Config, although the choice of OP2 backend may restrict the choice of linear algebra backends.

Solver options can be set by updating the parameters, for example:

solver.parameters.update({'type' : 'cg', 'preconditioner' : ..., ...})

When it comes to solving a system, the user uses:

# A is an op2.Mat object, x and b are op2.Dat objects
solver.solve(A, x, b)

Subsequently, the user might enquire about the solve:

solver.get_iteration_number() # Return the number of iterations executed
solver.get_converged_reason() # Return the reason the iteration stopped

These are modelled on PETSc's convergence reasons, which is probably adequate for now, but could be revisited in the future.

Planned implementation

For each solver backend, a Python class will be created. The correct class for the chosen backend will be instantiated by the op2.Solver class when it is constructed. The interface for the backend objects involves defining the following methods:

  • __init__(self): Performs any backend-specific initialisation
  • __dealloc__(self): Performs backend-specific cleanup
  • set_parameters(self, parameters): Takes a dict of parameters and sets these parameters in the linear algebra backend accordingly. Parameters include 'linear_solver' (the solver algorithm, e.g. 'cg' or 'gmres'), 'preconditioner' (the preconditioner, e.g. 'jacobi' or 'ilu'), tolerances including 'relative_tolerance', 'absolute_tolerance', and 'divergence_tolerance', and 'maximum_iterations', the number of iterations to stop at if convergence is not achieved.
  • solve(self, A, x, b): Solves the system Ax = b using the backend.
  • get_converged_reason(self): Returns the reason the iteration terminated in the previous solve. It is an error to call this method if there was no previous solve.
  • get_iteration_number(self): Returns the number of iterations executed in the previous solve. It is an error to call this method if there was no previous solve.

PETSc backend example

The class ksp_solver is an example implementation of a backend that uses the PETSc KSP libraries. It is implemented as a Cython interface over a C glue layer that abstracts away a lot of the data types used in PETSc. The C layer shields the Cython implementation from having to be aware of the structure of PETSc data structures, which would otherwise need duplicating into a Cython header file.

A bare minimum C interface would contain the following functions:

void* create_ksp()
{
  KSP ksp;
  KSPCreate(PETSC_COMM_WORLD,&ksp);
  KSPSetFromOptions(ksp);
  return (void*)ksp;
}

void destroy_ksp(void *ksp)
{
#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2)
  KSPDestroy((KSP*)(&ksp));
#else
  KSPDestroy((KSP)ksp);
#endif
}

int ksp_get_converged_reason(void *ksp)
{
  KSPConvergedReason reason;
  KSPGetConvergedReason((KSP)ksp, &reason);
  return (int)reason;
}

int ksp_get_iteration_number(void *ksp)
{
  int its;
  KSPGetIterationNumber((KSP)ksp, &its);
  return its;
}

void ksp_set_type(void *ksp, const char *type)
{
  KSPSetType((KSP)ksp, type);
}

void* ksp_get_pc(void *ksp)
{
  PC pc;
  KSPGetPC((KSP)ksp, &pc);
  return (void*)pc;
}

void pc_set_type(void *pc, const char *type)
{
  PCSetType((PC)pc, type);
}

void ksp_set_tolerances(void *ksp, double rtol, double atol, double dtol, int maxits)
{
  KSPSetTolerances((KSP)ksp, rtol, atol, dtol, maxits);
}

void ksp_solve(void *ksp, op_mat mat, op_dat x, op_dat b)
{
  assert( mat && b && x );

  Vec p_b = op_create_vec(b);
  Vec p_x = op_create_vec(x);
  Mat A = (Mat) mat->mat;

  KSPSetOperators((KSP)ksp,A,A,DIFFERENT_NONZERO_PATTERN);
  KSPSolve((KSP)ksp,p_b,p_x);

  op_destroy_vec(p_b, b);
  op_destroy_vec(p_x, x);
}

Assuming we then create a Cython header file called _op_lib_petsc.pxd for all these functions, we can then implement the ksp_solver object as:

cimport _op_lib_petsc as petsc

class ksp_solver:
    cdef void * _ksp
    cdef void * _pc

    def __cinit__(self):
        self._ksp = petsc.create_ksp()
        self._pc = petsc.ksp_get_pc(self._ksp)

    def __dealloc__(self):
        petsc.destroy_ksp(self._ksp)

    def set_parameters(self, parameters):
        ksptype = parameters['linear_solver']
        pctype  = parameters['preconditioner']
        rtol    = parameters['relative_tolerance']
        atol    = parameters['absolute_tolerance']
        dtol    = parameters['divergence_tolerance']
        maxits  = parameters['maximum_iterations']

        petsc.ksp_set_type(self._ksp, ksptype)
        petsc.pc_set_type(self._pc, pctype)
        petsc.ksp_set_tolerances(self._ksp, rtol, atol, dtol, maxits)

    def solve(self, A, x, b):
        cdef op_mat cA
        cdef op_dat cx, cb
        cA = A._c_handle
        cx = x._c_handle
        cb = b._c_handle
        petsc.ksp_solve(self._ksp, cA._handle, cx._handle, cb._handle)

    def get_converged_reason(self):
        return petsc.ksp_get_converged_reason(self._ksp)

    def get_iteration_number(self):
        return petsc.ksp_get_iteration_number(self._ksp)

An op2.Solver object that makes use of the ksp_solver class would look like the following (noting that this doesn't include a method for configuring different backends, but illustrates how the Solver object would interface with the backend solver object):

_DEFAULT_SOLVER_PARAMETERS = {'linear_solver':      'cg',
                              'preconditioner':     'jacobi',
                              'relative_tolerance': 1.0e-5,
                              'absolute_tolerance': 1.0e-50,
                              'divergence_tolerance': 1.0e+4,
                              'maximum_iterations': 10000 }

class Solver(object):
    def __init__(self, parameters=None):
        self.parameters = parameters or _DEFAULT_SOLVER_PARAMETERS.copy()
        self._ksp_solver = core.ksp_solver()

    def solve(self, A, x, b):
        self._ksp_solver.set_parameters(self.parameters)
        self._ksp_solver.solve(A, x, b)
        reason = self._ksp_solver.get_converged_reason()
        its = self._ksp_solver.get_iteration_number()
        print "Converged reason", reason
        print "Iterations", its

Demo implementation

An implementation of the PETSc backend in this form for the sequential OP2 backend is in the following two branches:

Note that this implementation is presently tangled up in the cpp-matrices implementation, but there is no definite need for all the code in cpp-matrices - this was the easiest place to build this experimental implementation quickly.

Why not PETSc4py?

  • For our use case, we only need to access a small subset of PETSc's functionality.
  • PETSc4py covers a much larger set of PETSc, but it is much more complicated than the simple PETSc implementation proposed here.
  • Additionally, PETSc4py doesn't provide a way to directly access the matrix data, which is a requirement for PyOP2.
  • Thus in order to make use of PETSc4py, it would have to be forked. This is slightly scary, since there is little documentation for PETSc4py, and it mirrors all of PETSc's data structures, so we'd be restricted to a specific PETSc version in our fork, and having to keep maintaining our fork for newer versions of PETSc.
  • Summary: using PETSc4py looks like much harder work for our use case.

Implementation using master branch of OP2-Common

Since we have decided not to integrate the matrix support into the OP2-Common master branch, we need to move all the matrix functionality into PyOP2. The above PETSc implementation theoretically has dependencies on code that is only in the cpp-matrices branch, but this could all be moved into PyOP2 so that we could make this work with the master branch.

The main change that we would need to be aware of is the removal of the op_sparsity and op_mat data types. The change in PyOP2 would probably involve the linear algebra backend interface being used to construct sparsity and matrix objects that it supports instead.

Present interface

At the moment, the core interface implements the core.op_sparsity with the following members and methods:

  • cdef core.op_sparsity _handle. This is a handle to a C-layer op_sparsity object, which has the following structure:
typedef struct {
  op_map *rowmaps;              /* row maps */
  op_map *colmaps;              /* column maps */
  int nmaps;                    /* number of maps */
  int dim[2];
  size_t nrows;       /* number of rows */
  size_t ncols;       /* number of columns */
  int    *d_nnz;      /* vector of number of nonzeros per row in
                       * diagonal subblock */
  int    *o_nnz;      /* vector of number of nonzeros per row in
                       * off-diagonal subblock */
  int    total_nz;    /* total nonzeros in the pattern sum(d_nnz + o_nnz) */
  int    *rowptr;     /* csr row pointer (accumulation of d_nnz and o_nnz) */
  int    *colidx;     /* csr column indices for each row */
  char const *name;   /* name of dataset */
} op_sparsity_core;

typedef op_sparsity_core * op_sparsity;
  • __cinit__(self,sparsity): Takes an op2.Sparsity object and uses it to instantiate the core op_sparsity to be stored in self._handle. This involves picking out the pointers to the C-Layer maps referenced in the sparsity, and then calls the C-layer op_decl_sparsity_core. The main function of op_decl_sparsity_core is to perform some argument checking, and then construct the sparsity pattern by calling op_build_sparsity_pattern - this constructs the rowptr and colidx arrays.
  • total_nz(self): returns the total number of nonzero elements in the sparsity pattern.
  • rowptr(self): returns a numpy array wrapping the rowptr member of the _handle
  • colidx(self): returns a numpy array wrapping the colidx member of the _handle

The core.op_mat is implemented with the following methods and members:

  • cdef core.op_mat _handle. This is a handle to a C-layer op_mat object, which has the following structure:
typedef struct
{
  int         index;  /* index */
  int         dim[2], /* dimension of data */
              size;   /* size of each element in dataset */
  void*       mat;    /* handle to matrix object from 3rd party library */
  void*       mat_array;        /* handle to matrix value from 3rd party library */
  char const *type,   /* datatype */
             *name;   /* name of matrix */
  op_sparsity sparsity;         /* sparsity structure */
  char       *data;   /* matrix entries (device) */
  char       *lma_data;         /* matrix entries in LMA form (device) */
} op_mat_core;

typedef op_mat_core * op_mat;
  • _nnzeros: This stores the total number of nonzeroes in the matrix - this is a multiple of the sparsity total_nz property, since the matrix entries may have dimension greater than 1
  • __cinit__(self, mat): Takes an op2.Mat object and initialises a C-layer op_mat object from it by calling op_decl_mat. This creates either an MPIAIJ or SeqAIJ PETSc matrix object with the given sparsity, depending on whether there is more than one MPI process.
  • zero(self): Zeroes all entries in the matrix
  • zero_rows(self, rows, v): Zeroes all the rows in the rows list and sets the diagonal entry to v.
  • assemble(self): Calls the PETSc matrix assembly begin and end routines to assemble the matrix throigh core.op_mat_assemble. This is called by the OpenCL backend to put the matrix in the assembled statement, even though it does not call MatSetValues by any path.
  • array(self): property that returns the value array of the matrix wrapped in a numpy array. As far as I can determine, the semantics of this mean that PETSc no longer considers itself the owner of the array: http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetArray.html#MatGetArray
  • restore_array(self): Relinquishes ownership of the array data after array() has been called and there is no longer any need to access the array data.
  • cptr(self): Returns a C pointer to _handle (an op_mat) so that the sequential backend can pass it to op_mat_addto
  • values(self): Returns a dense numpy array containing the matrix values.

In addition to the objects presented to the Python interface, there are also addto functions that are exposed such that the sequential backend can call them in the generated code. These are:

  • op_mat_addto_scalar( op_mat mat, const void* value, int row, int col ): Add the value into the entry at (row,col) in the matrix mat. Calls the PETSc MatSetValues interface.
  • op_mat_addto( op_mat mat, const void* values, int nrows, const int *irows, int ncols, const int *icols ): Adds a vector of values to the entries pointed to by irows and icols into mat. Calls the PETSc MatSetValues interface.
  • op_mat_assemble( op_mat mat ): Calls the Petsc MatAssemblyBegin and MatAssemblyEnd interfaces to assemble the matrix after the addtos have been performed.

These functions are exposed to the generated code through mat_utils.cxx, which Instant compiles and links with the generated code:

#include "op_lib_mat.h"
#include "mat_utils.h"

void addto_scalar(op_mat mat, const void *value, int row, int col)
{
  op_mat_addto_scalar(mat, value, row, col);
}

void addto_vector(op_mat mat, const void *values,
                  int nrows, const int *irows,
                  int ncols, const int *icols)
{
    op_mat_addto(mat, values, nrows, irows, ncols, icols);
}

void assemble_mat(op_mat mat)
{
    op_mat_assemble(mat);
}

The PyOP2 Mat object constructs an underlying core.op_mat object when its _c_handle is accessed:

    @property
    def _c_handle(self):
        if self._lib_handle is None:
            self._lib_handle = core.op_mat(self)
        return self._lib_handle

The PyOP2 Sparsity object behaves similarly:

    @property
    def _c_handle(self):
        if self._lib_handle is None:
            key = (self._rmaps, self._cmaps, self._dims)
            self._lib_handle = _sparsity_cache.get(key) or core.op_sparsity(self)
            _sparsity_cache[key] = self._lib_handle
        return self._lib_handle

Proposed changes

The op_sparsity and op_mat declarations in op_lib_core.pyx need to be replaced with linear algebra backend-specific interfaces that contain code necessary to interact with the library without depending on any code specific to the cpp-matrices branch. As with the solver interface, the correct interface object will need to be instantiated depending on the OP2 configuration.

For example: we may begin by replacing op_mat and op_sparsity with a generic sparsity Python class, and a PETSc-specific petsc_mat class. The Python sparsity class could be generic since the current op_sparsity used in the cpp-matrices branch doesn't have any PETSc-specific functionality. The petsc_mat class will present the same interface as the current core.op_mat, but will call methods to create the matrix and manipulate it through a small C stub, in a similar fashion to the stub for accessing PETSc's KSP* methods in the example above.

The mat_utils.cxx stub will need changing so that instead of operating on op_mat objects, it accepts and works with PETSc Mat objects instead. Similarly, the cptr method of petsc_mat should return a pointer to the PETSc Mat rather than an op_mat. The op_mat_addto* functions will also take a PETSc Mat instead of an op_mat. These should be relatively simple changes, since the current implementations that take op_mat objects have to unbox the Mat objects from the op_mats.

Additionally, the PyOP2 Mat object will need to construct a petsc_mat instead of a core.op_mat:

    @property
    def _c_handle(self):
        if self._lib_handle is None:
            self._lib_handle = petsc.petsc_mat(self)
        return self._lib_handle

However, the PyOP2 Sparsity object may not need another object at the linear algebra backend layer, since its members may be generic across backends. In this case, the sparsity-building functions that are presently written in C in the cpp-matrices branch may be ported into the Python Sparsity class implementation.

Clone this wiki locally