PETSc的使用

运行 PETSc 程序

mpirun -n 8 ./petsc_program_name petsc_options

PETSc 程序支持以下选项：

log_view : summarize the program’s performance
fp_trap : stop on floating-point exceptions; for example divide by zero
malloc_dump : enable memory tracing; dump list of unfreed memory at conclusion of the run, see Detecting Memory Allocation Problems,
malloc_debug : enable memory debugging (by default this is activated for the debugging version of PETSc), see Detecting Memory Allocation Problems,
start_in_debugger [noxterm,gdb,lldb] [-display name] : start all processes in debugger. See Debugging, for more information on debugging PETSc programs.
on_error_attach_debugger [noxterm,gdb,lldb] [-display name] : start debugger only on encountering an error
info : print a great deal of information about what the program is doing as it runs

编写 PETSc 程序

1
2
3


ierr = PetscInitialize(int *argc,char ***argv,char *file,char *help);
if (ierr)
    return ierr;

初始化PETSc和MPI。默认会将PETSC_COMM_WORLD设为MPI_COMM_WORLD。所有的PETSc操作都会返回一个PetscErrorCode，该值为0代表检测没有错误，非0代表错误发生。

Vectors and Parallel Data

VecCreateSeq(PETSC_COMM_SELF,PetscInt m,Vec *x);

创建序列化的向量

VecCreateMPI(MPI_Comm comm,PetscInt m,PetscInt M,Vec *x);

在通信域comm中创建分布式的向量,m代表了向量存在local的大小，M代表总的向量大小。可以将m``M其中一个设置为PETSC_DECIDE或PETSC_DETERMINE。

VecSet(Vec x,PetscScalar value);

将向量中的所有值设为value

1
2
3
4


VecSetValues(Vec x,PetscInt n,PetscInt *indices,PetscScalar *values,INSERT_VALUES);
VecAssemblyBegin(Vec x);
...
VecAssemblyEnd(Vec x);

n代表插入的数据的数目，indices代表插入向量的全局位置，values代表插入的值

VecView(Vec x,PetscViewer v);

为了打印到屏幕上，可以使用viewerPETSC_VIEWER_STDOUT_WORLD

VecDuplicate(Vec old,Vec *new);

复制一个相同的向量

VecDuplicateVecs(Vec old,PetscInt n,Vec **new);

复制多个相同的向量

VecDestroy(Vec *x);

销毁向量

VecDestroyVecs(PetscInt n,Vec **vecs);

销毁向量数组

VecCreateSeqWithArray(PETSC_COMM_SELF,PetscInt bs,PetscInt n,PetscScalar *array,Vec *V);

VecCreateMPIWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,PetscInt N,PetscScalar *array,Vec *vv);

基于用户提供的数组创建向量

VecGetOwnershipRange(Vec vec,PetscInt *low,PetscInt *high);

得到process所拥有的向量部分的范围,low代表起始位置， high代表结束位置

VecGetArray(Vec v,PetscScalar **array);

获取向量中的元素

VecRestoreArray(Vec v, PetscScalar **array);

当获取的array不再需要

VecGetArrayRead(Vec v, const PetscScalar **array);

VecRestoreArrayRead(Vec v, const PetscScalar **array);

只读模式获取元素

VecGetLocalSize(Vec v,PetscInt *size);

获取本地元素个数

VecGetSize(Vec v,PetscInt *size);

获取全局元素个数

1
2
3
4


VecScatterCreate(Vec x,IS ix,Vec y,IS iy,VecScatter *ctx);
VecScatterBegin(VecScatter ctx,Vec x,Vec y,INSERT_VALUES,SCATTER_FORWARD);
VecScatterEnd(VecScatter ctx,Vec x,Vec y,INSERT_VALUES,SCATTER_FORWARD);
VecScatterDestroy(VecScatter *ctx);

Scatter 和 gathers

映射问题

本地标识和全局标识的映射

1
2
3
4


ISLocalToGlobalMappingCreate(MPI_Comm comm,PetscInt bs,PetscInt N,PetscInt* globalnum,PetscCopyMode mode,ISLocalToGlobalMapping* ctx);
ISLocalToGlobalMappingApply(ISLocalToGlobalMapping ctx,PetscInt n,PetscInt *in,PetscInt *out);
ISLocalToGlobalMappingApplyIS(ISLocalToGlobalMapping ctx,IS isin,IS* isout);
ISLocalToGlobalMappingDestroy(ISLocalToGlobalMapping *ctx);

基本的向量操作函数

Matrices

创建和收集

创建矩阵

1
2


MatCreate(MPI_Comm comm,Mat *A)
MatSetSizes(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N)

插值

MatSetValues(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[], const PetscScalar values[],INSERT_VALUES);

MatSetValues(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar values[],ADD_VALUES);

矩阵处理需要调用

MatAssemblyBegin(Mat A,MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(Mat A,MAT_FINAL_ASSEMBLY);

参考

PETSc manual