一.安装/配置PETSc

PETSc需要MPI和BLAS库。还有gcc等基础包。mpich可以直接apt-get安装：

1
sudo apt-get install mpich

如果不确定有没有安装BLAS/mpich，可以apt list —installed | grep blas找一下。

有了这两个之后下载petsc，解压后执行./configure，接下来按照给出的提示完成安装就可以了。如果没有mpi或者blas，可以在后面加--download-mpich --download-fblaslapack。

PETSC有两个可设置的环境变量，分别是：

1
2
PETSC_DIR=/path/to/petsc
PETSC_ARCH=/

PETSC这个环境变量的作用是指定特定的构建和配置，一般来说，在一般机器上会构建多个PETSC版本。默认情况下(按照上面的直接安装)是针对调试使用的，其路径为arch-linux-c-debug。还可以构建一个针对性能优化的版本(例如arch-linux-opt)。

以上已经编译了一个调试使用的版本，接下来构建一个arch-linux-opt配置的版本，和上面是类似的，运行以下命令，然后根据给出的提示完成安装配置：

1
./configure PETSC_ARCH=arch-linux-opt --with-debugging=no COPTFLAGS="-O3 -march=armv8-a" CXXOPTFLAGS="-O3 -march=armv8-a" FOPTFLAGS="-O3 --march=armv8-a"

二.编译/运行PETSc程序

上述工作后，PETSc就编译好了。编写PETSC程序时引入头文件，就可以使用相应的API。编译PETSC程序前，都要保证环境变量中有PETSC_DIR和PETSC_ARCH，可以在编译前进行export，也可以直接写入~/.bashrc。然后使用以下命令编译程序：

1
mpicxx -o main main.cpp -I${PETSC_DIR}/include -I${PETSC_DIR}/${PETSC_ARCH}/include -L${PETSC_DIR}/${PETSC_ARCH}/lib -lpetsc

运行的时候可能会发现报error while loading shared libraries: libpetsc.so.3.20: cannot open shared object file: No such file or directory，PETSC动态库没有链接上，遇到这种情况可以将其添加到动态库路径：

1
2
3
#找不到.so的位置可以find一下：
find / -name libpetsc.so.3.20	
export LD_LIBRARY_PATH=/root/petsc-3.20.5/arch-linux-opt/lib:LD_LIBRARY_PATH

然后就可以运行程序了：

1
mpirun -n 1 ./main

三.PETSc的API

PETSc的API可以在petsc.org找到，文档编写的非常详细。下面以用PETSc解决SpMV问题为例，为了方便没有使用PetscCall()处理可能的错误，而且是编写SpMV测试时使用的，因此值是从定义好的矩阵格式获取再SetValue插入的，PETSc还支持SetValues直接插入数组来填充矩阵值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
void petsc_evaluation(SparseMatrixCOO *ordered_coo, char *output_path){    
    
    Mat A;
    PetscMPIInt rank;
    double *vecx = (double *)malloc(sizeof(double)*ordered_coo->ncolumns);
    struct timespec start, end;
    long long unsigned int elapsed;

    PetscFunctionBeginUser;
    PetscInitialize(NULL, NULL, NULL, NULL);
    MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
    MatCreate(PETSC_COMM_WORLD, &A);
    MatSetType(A, MATMPIAIJ);
    MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, ordered_coo->nrows, ordered_coo->ncolumns);
    MatSetFromOptions(A);
    MatSetUp(A);
    // 填充矩阵的值
    for(int i=0;i<ordered_coo->nnz;i++){
        MatSetValue(A, (PetscInt)ordered_coo->row_index[i], (PetscInt)ordered_coo->col_index[i], (PetscScalar)ordered_coo->values[i], INSERT_VALUES);
    }
    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

    //创建向量x
    Vec x;
    VecCreate(PETSC_COMM_WORLD, &x);
    VecSetSizes(x, PETSC_DECIDE, ordered_coo->ncolumns);
    VecSetFromOptions(x);
    generate_x(vecx, ordered_coo->ncolumns);
    for(int i=0;i<ordered_coo->ncolumns;i++) VecSetValue(x, (PetscInt)i, (PetscScalar)vecx[i], INSERT_VALUES);
    VecAssemblyBegin(x);
    VecAssemblyEnd(x);

    Vec y;
    VecCreate(PETSC_COMM_WORLD, &y);
    VecSetSizes(y, PETSC_DECIDE, ordered_coo->ncolumns);
    VecSetFromOptions(y);
    VecAssemblyBegin(y);
    VecAssemblyEnd(y);

    int iterations;
    if(ordered_coo->nnz>150000) iterations = 5000;
    else if(ordered_coo->nnz>300000) iterations = 1000;
    else iterations = 50000;
    #ifdef DEBUG
        iterations = 1;
    #endif

    //SpMV kernel
    clock_gettime(CLOCK_MONOTONIC_RAW, &start);
    for(int i=0; i<iterations; i++) {
        MatMult(A, x, y);
        PetscBarrier(NULL);
    }
    clock_gettime(CLOCK_MONOTONIC_RAW, &end);
    elapsed = (end.tv_sec - start.tv_sec) * 1000000 + (end.tv_nsec - start.tv_nsec) / 1000;
    //gflops
    clock_gettime(CLOCK_MONOTONIC_RAW, &start);
    MatMult(A, x, y);
    PetscBarrier(NULL);
    clock_gettime(CLOCK_MONOTONIC_RAW, &end);
    double flops_count = ordered_coo->nnz*2;
    double sec = (end.tv_sec - start.tv_sec) + (end.tv_nsec - start.tv_nsec) / 1e9;
    double GFLOPS =  flops_count / (sec*1e9);

    if(rank==0){
        FILE *fp;
        fp = fopen(output_path, "a");
        if (fp == NULL) {
            printf("无法打开文件 %s\n", output_path);
        }
        fprintf(fp, "%s %llu %lf\n", "petsc", elapsed, GFLOPS);
        fclose(fp);
    }
    // 打印结果
    //MatView(A, PETSC_VIEWER_STDOUT_WORLD);
    //VecView(x, PETSC_VIEWER_STDOUT_WORLD);
    //VecView(y, PETSC_VIEWER_STDOUT_WORLD);

    free(vecx);
    VecDestroy(&x);
    VecDestroy(&y);
    MatDestroy(&A);
    PetscFinalize();
}