Petsc入门基础篇

欧拉常数的计算

本人之前讲述过在虚拟机环境中安装Petsc，下面都默认在已经有虚拟机的情况运行。
$e={\sum }_{n=0}\frac{1}{n!}$
下面我们先通过一个简单的代码来介绍一些基础知识：

#include <petsc.h>

int main(int argc,char **argv)
{
    PetscErrorCode ierr;
    PetscMPIInt rank;
    PetscMPIInt i;
    PetscReal localval,globalsum;

    PetscInitialize(&argc,&argv,NULL,"compute e in parallel in petsc.\n");
    ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);

    localval = 1.0;
    for(i = 2; i < rank + 1; i++)
	localval = localval/i;

    ierr = MPI_Allreduce(&localval,&globalsum,1,MPIU_REAL,MPIU_SUM,
	PETSC_COMM_WORLD);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"e is %17.15f\n",globalsum);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_SELF,"rank %d did %d flops \n",
	rank,(rank > 0)? rank - 1:0);CHKERRQ(ierr);
    return PetscFinalize();
}

上面有一些换行涉及空格，本人复制过程容易混乱，读者自己注意
上来我们第一个看见的就是头文件petsc.h，这个就和C语言中的头文件stdio.h一样。
在main()函数里面有两个参数，第一个是整型变量，第二个是字符串型数组。所有的代码的参数都是传入PetscInitialize()里面，其中PetscInitialize()一般在主程序开头，用于启动Petsc环境，PetscFinalize()在主程序结尾，用于终止Petsc环境。误差捕捉宏CHKERRQ(ierr)。
在初始化环境之前，我们定义了一些变量，这些我们可以直接按照C语言风格来理解，PetscInt=int,PetscReal=double。下面把makefile文件内容也放上去。
上来我们第一个看见的就是头文件petsc.h，这个就和C语言中的头文件stdio.h一样。
然后定义一个main函数，在main函数里面有两个参数，第一个是整型int变量，第二个是字符串型数组char变量，这个框架和c语言中编写hello world一模一样。
然后开始申明变量，变量的类型分别是PetscErrorCode, PetscMPIInt, PetscReal，主要介绍后面两个变量类型，其中PetscMPIInt相当于c语言中的整型int，PetscReal相当于c语言中的双精度类型double。说明rank,i都是整数，localval,globalsum都是双精度浮点数。这里rank表示运行这个代码调用的进程数目，这里的rank表示的其实也是参与计算欧拉常数的整数数目。这里的i表示的是索引，localval在后面会发现这个变量其实记录的是1/n!，globalsum记录的是每一个进程存储的localval的求和。PetscErrorCode是误差捕捉宏。
申明完需要用到的变量以后，现在开始启动pestc环境，这里用的是PetscInitialize()，代码里面的是PetscInitialize()的标准用法，不管接下来的petsc代码如何编写，PetscInitialize()始终可以保持不变放在主程序的开头，启动petsc环境以后就可以运行代码，运行完以后需要调用PetscFinalize()终止petsc环境，PetscFinalize()只需要放在代码末尾就可以终止，这里需要注意的是PetscInitialize()和PetscFinalize()一定是成对出现的。
PETSC_COMM_WORLD表示MPI进程的通信器，是进程的集合，ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);这一句是在告诉程序调用多少个进程来执行这个代码，也就是告诉程序这里一共用了多少个求和项。代码的主要部分是是那个for循环，这里特别说明，每一个进程都会执行一个for循环：
localval = 1.0;
for(i = 2; i < rank + 1; i++)
localval = localval/i;
这里是在计算每个进程存储的localval数值，我们可以把一个进程看成一个独立的处理器，这个处理器可以独立做简单的四则运算，并且这个进程还有一定的存储能力，可以存储少量数据，这里详细介绍一下这个for循环：
首先需要说明，rank从0开始，如果我们传入的n=5，那么调用的rank只有0，1，2，3，4.
当rank=0时，无法进入循环，localval=1。
当rank=1时，仍让无法进入循环，localval=1。
当rank=2时，此时进入循环发现i只会取2，localval=1/2。
当rank=3时，此时进入循环发现i可以取2,3，当i=2时，根据localval = localval/i;得到此时localval=1/2，然后i=3，新的localval=1/2/3=1/6。
当rank=4时，此时进入循环发现i可以取2,3,4，根据rank=3的经验，我们会发现localval从1更新到了1/2(i=2)，然后更新到了1/6(i=3)，最后更新到了1/24(i=4)。
现在每一个进程的执行都结束了，我们发现0号进程(rank=0)存储的localval=1，1号进程(rank=1)存储的localval=1，2号进程(rank=2)存储的localval=1/2，3号进程(rank=3)存储的localval=1/6，4号进程(rank=4)存储的localval=1/24。现在每个进程上存储的localval都有了，那么问题就变成了怎么把不同进程上存储的信息做一个求和。
这个求和用的就是MPI里面的一个规约函数Allreduce，具体用法就是下面这行代码。
ierr = MPI_Allreduce(&localval,&globalsum,1,MPIU_REAL,MPIU_SUM,
PETSC_COMM_WORLD);CHKERRQ(ierr);
在这行代码里面，调用的时MPI_Allreduce函数，这个函数接受的参数分别是localval,globalsum的地址，MPIU_REAL表示参与规约运算的都是实数类型变量，MPIU_SUM表示这里的规约操作用的时求和。事实上MPI的规约操作还包括求平均值，求最大值最小值，这里不做具体介绍。这个MPI_Allreduce接受的最后一个参数PETSC_COMM_WORLD表示这里的规约操作是针对所有的进程。
后面两行代码都是负责打印输出，没有什么难点。

include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
CFLAGS += -pedantic -std=c99

e: e.o
	-${CLINKER} -o e e.o  ${PETSC_LIB}
	${RM} e.o

# testing
rune_1:
	-@../testit.sh e "" 1 1

test_e: rune_1

test: test_e

.PHONY: distclean rune_1 test_e test

distclean:
	@rm -f *~ e *tmp

事实上，后面所有的petsc代码都是这个结构。
代码的编译和执行
输入:
make e
./e或者mpiexec -n 5 ./e
rank表示MPI的进程号。第二个mpiexec -n 5 ./e指的是调用5个MPI进程(process)来估算欧拉常数e。
其中flop表示float point operations，浮点数运算。
PETSC_COMM_WORLD表示MPI进程的通信器，是一个进程的集合。每一个MPI进程都只算了一个$\frac{1}{n!}$，其中通过MPI_Comm()函数来返回n的取值。
如果我们使用N个进程，那么代码就调用MPI_Allreduce()函数来计算第N个部分和，然后把结果返回给每个进程。
在打印结果的时候，我们有两个打印过程，第一个打印只会在0号进程工作，所以只有一行输出。第二个打印在每个进程都工作。但是第二个打印PESTC_COMM_SELF的打印是随机的，换句话说，一般不会先打印1号进程，2号进程，顺序是随机的。如果想要按顺序打印，那么需要使用PetscSynchronizedPrintf()。但是这个需要额外在进程之间增加通信成本。
下面这个是mpiexec -n 7 ./e的结果