PETSc源码分析：Algebraic Objects

本文结合PETSc源代码，旨在分析PETSc中的基础数据对象。

注1：限于研究水平，分析难免不当，欢迎批评指正。

注2：文章内容会不定期更新。

一、Design Rationales in PETSc

1.1 OOD: Object-Oriented Design

PETSc使用C语言实现了封装、继承、多态等部分面向对象编程特性。

• Encapsulation

在PETSc中，将PetscObject定义为描述类型信息的元对象数据指针，可以看到PetscObject主要包含两部分内容，一是类型相关成员变量；二是类型相关的成员函数，藉此完成了面向对象中‘封装’特性。

// include/petscsystypes.h
typedef struct _p_PetscObject *PetscObject;

// include/petsc/private/petscimpl.h
/*
   All PETSc objects begin with the fields defined in PETSCHEADER.
   The PetscObject is a way of examining these fields regardless of
   the specific object. In C++ this could be a base abstract class
   from which all objects are derived.
*/
#define PETSC_MAX_OPTIONS_HANDLER 5
typedef struct _p_PetscObject {
  PetscOps      bops[1];
  PetscClassId  classid;
  MPI_Comm      comm;
  PetscObjectId id; /* this is used to compare object for identity that may no longer exist since memory addresses get recycled for new objects */
  PetscInt      refct;
  PetscErrorCode (*non_cyclic_references)(PetscObject, PetscInt *);
  PetscInt64        cidx;
  PetscMPIInt       tag;
  PetscFunctionList qlist;
  PetscObjectList   olist;
  char             *class_name; /*  for example, "Vec" */
  char             *description;
  char             *mansec;
  char             *type_name; /*  this is the subclass, for example VECSEQ which equals "seq" */
  char             *name;
  char             *prefix;
  PetscInt          tablevel;
  void             *cpp;
  PetscObjectState  state;
  PetscInt          int_idmax, intstar_idmax;
  PetscObjectState *intcomposedstate, *intstarcomposedstate;
  PetscInt         *intcomposeddata, **intstarcomposeddata;
  PetscInt          real_idmax, realstar_idmax;
  PetscObjectState *realcomposedstate, *realstarcomposedstate;
  PetscReal        *realcomposeddata, **realstarcomposeddata;
#if PetscDefined(USE_COMPLEX)
  PetscInt          scalar_idmax, scalarstar_idmax;
  PetscObjectState *scalarcomposedstate, *scalarstarcomposedstate;
  PetscScalar      *scalarcomposeddata, **scalarstarcomposeddata;
#endif
  void (**fortran_func_pointers)(void);             /* used by Fortran interface functions to stash user provided Fortran functions */
  PetscFortranCallbackId num_fortran_func_pointers; /* number of Fortran function pointers allocated */
  PetscFortranCallback  *fortrancallback[PETSC_FORTRAN_CALLBACK_MAXTYPE];
  PetscFortranCallbackId num_fortrancallback[PETSC_FORTRAN_CALLBACK_MAXTYPE];
  void                  *python_context;
  PetscErrorCode (*python_destroy)(void *);

  PetscInt noptionhandler;
  PetscErrorCode (*optionhandler[PETSC_MAX_OPTIONS_HANDLER])(PetscObject, PetscOptionItems *, void *);
  PetscErrorCode (*optiondestroy[PETSC_MAX_OPTIONS_HANDLER])(PetscObject, void *);
  void *optionctx[PETSC_MAX_OPTIONS_HANDLER];
#if defined(PETSC_HAVE_SAWS)
  PetscBool amsmem;          /* if PETSC_TRUE then this object is registered with SAWs and visible to clients */
  PetscBool amspublishblock; /* if PETSC_TRUE and publishing objects then will block at PetscObjectSAWsBlock() */
#endif
  PetscOptions options; /* options database used, NULL means default */
  PetscBool    optionsprinted;
  PetscBool    donotPetscObjectPrintClassNamePrefixType;
} _p_PetscObject;

#define PETSCHEADER(ObjectOps) \
  _p_PetscObject hdr; \
  ObjectOps      ops[1]

• Inheritance

PETSc提供了PETSCHEADER宏用于在PetscObject基础扩展类型定义，这样，在逻辑上，PETSc中大多数类型均派生于PetscObject。

对于Vec、Mat等类型描述，参见下文。

• Composite

PetscObject使用_p_PetscObject::olist来记录子对象，进而实现对象复合。同时提供了PetscObjectCompose接口构建对象间的复合关系。

// src/sys/objects/inherit.c
PetscErrorCode PetscObjectCompose(PetscObject obj, const char name[], PetscObject ptr)
{
  PetscFunctionBegin;
  PetscValidHeader(obj, 1);
  PetscAssertPointer(name, 2);
  if (ptr) PetscValidHeader(ptr, 3);
  PetscCheck(obj != ptr, PetscObjectComm((PetscObject)obj), PETSC_ERR_SUP, "Cannot compose object with itself");
  if (ptr) {
    char     *tname;
    PetscBool skipreference;

    PetscCall(PetscObjectListReverseFind(ptr->olist, obj, &tname, &skipreference));
    if (tname) PetscCheck(skipreference, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "An object cannot be composed with an object that was composed with it");
  }
  PetscCall(PetscObjectListAdd(&obj->olist, name, ptr));
  PetscFunctionReturn(PETSC_SUCCESS);
}

• Polymorphism

基于上述封装、继承机制，通过修改成员变量函数指针，便可在形式上实现"多态"机制。

1.2 MPI

二、Vectors

2.1 Class Outline

Vec在逻辑上继承自PetscObject，增加了向量相关的内存管理、缓存、向量操作等内容。

// include/petscvec.h
typedef struct _p_Vec *Vec;

// include/petsc/private/vecimpl.h
struct _p_Vec {
  PETSCHEADER(struct _VecOps);
  PetscLayout map;
  void       *data; /* implementation-specific data */
  PetscBool   array_gotten;
  VecStash    stash, bstash; /* used for storing off-proc values during assembly */
  PetscBool   petscnative;   /* means the ->data starts with VECHEADER and can use VecGetArrayFast()*/
  PetscInt    lock;          /* lock state. vector can be free (=0), locked for read (>0) or locked for write(<0) */
#if PetscDefined(USE_DEBUG)
  PetscStack lockstack; /* the file,func,line of where locks are added */
#endif
  PetscOffloadMask offloadmask; /* a mask which indicates where the valid vector data is (GPU, CPU or both) */
#if defined(PETSC_HAVE_DEVICE)
  void     *spptr; /* this is the special pointer to the array on the GPU */
  PetscBool boundtocpu;
  PetscBool bindingpropagates;
  size_t    minimum_bytes_pinned_memory; /* minimum data size in bytes for which pinned memory will be allocated */
  PetscBool pinned_memory;               /* PETSC_TRUE if the current host allocation has been made from pinned memory. */
#endif
  char *defaultrandtype;
};

// include/petsc/private/vecimpl.h
typedef struct _VecOps *VecOps;
struct _VecOps {
  PetscErrorCode (*duplicate)(Vec, Vec *);                                          /* get single vector */
  PetscErrorCode (*duplicatevecs)(Vec, PetscInt, Vec **);                           /* get array of vectors */
  PetscErrorCode (*destroyvecs)(PetscInt, Vec[]);                                   /* free array of vectors */
  PetscErrorCode (*dot)(Vec, Vec, PetscScalar *);                                   /* z = x^H * y */
  PetscErrorCode (*mdot)(Vec, PetscInt, const Vec[], PetscScalar *);                /* z[j] = x dot y[j] */
  PetscErrorCode (*norm)(Vec, NormType, PetscReal *);                               /* z = sqrt(x^H * x) */
  PetscErrorCode (*tdot)(Vec, Vec, PetscScalar *);                                  /* x'*y */
  PetscErrorCode (*mtdot)(Vec, PetscInt, const Vec[], PetscScalar *);               /* z[j] = x dot y[j] */
  PetscErrorCode (*scale)(Vec, PetscScalar);                                        /* x = alpha * x   */
  PetscErrorCode (*copy)(Vec, Vec);                                                 /* y = x */
  PetscErrorCode (*set)(Vec, PetscScalar);