Xloudmax

Note

0%

petsc_solver优化

发表于 分类于

petsc_solver优化
编译petsc_solver

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
#!/bin/bash

# 启用错误捕捉，遇到错误立即退出
set -e

# 设置编译器为 Intel 的 mpiicc 和 mpiicpc
export CC=mpiicx
export CXX=mpiicpx
export PETSC_DIR=/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/petsc-3.19.3
export PETSC_ARCH=petsc_install
# 设置 LAPACK 的头文件和库路径
export CPATH=/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/lapack-3.11/CBLAS/include:/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/lapack-3.11/LAPACKE/include:${CPATH}
export LD_LIBRARY_PATH=/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/lapack-3.11:${LD_LIBRARY_PATH}

# 清理并创建新的构建目录
rm -rf build 
mkdir build
cd build

# 配置项目并生成构建文件
cmake ..

# 使用并行编译，利用所有可用核心
make -j$(nproc)

# 将生成的库移动到指定目录
mv libpetsc_solver.a ../lib/

cmakelists.txt

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
cmake_minimum_required(VERSION 3.1)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

# 设置为 Release 模式，并启用优化选项
set(CMAKE_CONFIGURATION_TYPES "Release;Debug" CACHE STRING "Configs" FORCE)
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 -march=native")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -g")

if(NOT CMAKE_BUILD_TYPE)
    set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Build type" FORCE)
endif()
message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")

# 如果使用 Intel 编译器，请添加额外的优化标志
if (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM" OR CMAKE_CXX_COMPILER_ID STREQUAL "Intel")
    set(CMAKE_C_FLAGS_RELEASE "${CMAKE_C_FLAGS_RELEASE} -O3 -xHost -qopt-report=5")
    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 -xHost -qopt-report=5")
endif()

# 设置 PETSc 库路径
set(PETSC_DIR "/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/petsc-3.19.3/")
set(PETSC_ARCH "petsc_install") 
set(PETSC_LIBRARIES "/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/petsc-3.19.3/petsc_install/lib/libpetsc.so")
set(PETSC_INCLUDE_DIRS "/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/petsc-3.19.3/include" "/home/shenyufei/OpenCAEPoro/OpenCAEPoro_ASC2024/petsc-3.19.3/petsc_install/include")

project(petsc_solver)

include_directories(${PROJECT_SOURCE_DIR}/src/)
include_directories(${PROJECT_SOURCE_DIR}/include/)
include_directories(${PETSC_INCLUDE_DIRS})

file(GLOB_RECURSE SRC_FILES ${PROJECT_SOURCE_DIR}/src/*.cpp)

add_library(petsc_solver STATIC ${SRC_FILES})
target_link_libraries(petsc_solver ${PETSC_LIBRARIES})

# 安装目标
set(CMAKE_INSTALL_LIBDIR "${PROJECT_SOURCE_DIR}/lib/")
install(TARGETS petsc_solver
        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})

优化一：

export CC=mpiicx
export CXX=mpiicpx
显示指定了编译器为mpiicx和mpiicpx

优化二：

1
2
3
4
# 设置为 Release 模式，并启用优化选项
set(CMAKE_CONFIGURATION_TYPES "Release;Debug" CACHE STRING "Configs" FORCE)
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3 -march=native")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -Wall -g")

优化三：

使用PetscMalloc1和PetscFree等petsc自带的方法管理内存，确保与 PETSc 的内存管理机制兼容。

源代码

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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
start = MPI_Wtime(); // 开始记录解耦部分的时间
// 调用 decoup 函数进行解耦操作，处理稀疏矩阵的值，依据不同的解耦类型对矩阵进行处理。
// 其中 val 表示矩阵的值，rhs 表示右手边向量，rpt 和 cpt 分别表示行和列的指针，
// nb 表示块大小，nBlockRows 是当前处理的块行数，DecoupType 表示解耦的类型，is_thermal 表示是否是热问题。
decoup(val, rhs, rpt, cpt, nb, nBlockRows, DecoupType, is_thermal);

//----------------------------------------------------------------------
// Initializing MAT and VEC
// 初始化矩阵 A 和向量 x, b

Vec x, b; /* x 表示近似解，b 表示右手边向量 */
Mat A;    /* A 表示线性系统矩阵 */
KSP ksp;  /* KSP 表示线性求解器上下文 */
PC pc;    /* PC 表示预处理器上下文 */
PetscInt Ii, iters; // PETSc 的整型变量，用于矩阵行数和迭代次数
PetscErrorCode ierr; // PETSc 的错误代码

// 计算对角块和非对角块的元素数
for (int i = 0; i < nBlockRows; i++) {
    nDCount[i] = 0; // 初始化对角块元素数为 0
    // 遍历每一行，统计每一行中属于对角块的非零元素个数
    for (int j = rpt[i]; j < rpt[i + 1]; j++) {
        if (cpt[j] >= Istart && cpt[j] <= Iend) {
            nDCount[i]++;
        }
    }
}

// 计算每一行的非对角块的元素个数
for (int i = 0; i < nBlockRows; i++) {
    nNDCount[i] = rpt[i + 1] - rpt[i] - nDCount[i]; // 非对角块元素数 = 总元素数 - 对角块元素数
}

// 创建块压缩稀疏矩阵（Block AIJ 格式，使用并行通信 PETSC_COMM_WORLD）
ierr = MatCreateBAIJ(PETSC_COMM_WORLD, blockSize, rowWidth, rowWidth, matrixDim, matrixDim, 0, nDCount, 0, nNDCount, &A);
CHKERRQ(ierr); // 检查是否有错误

// 从命令行选项设置矩阵，并对矩阵进行初始化
ierr = MatSetFromOptions(A); // 从选项中读取设置
CHKERRQ(ierr);
ierr = MatSetUp(A); // 初始化矩阵，分配必要的内存
CHKERRQ(ierr);

// 将值填入矩阵 A
double *valpt = val; // 指向矩阵的值
int *globalx = (int *)malloc(blockSize * sizeof(int)); // 用于存储全局 x 方向的索引
int *globaly = (int *)malloc(blockSize * sizeof(int)); // 用于存储全局 y 方向的索引

// 遍历每一行的块并设置矩阵的值
for (Ii = 0; Ii < nBlockRows; Ii++) {
    for (i = 0; i < blockSize; i++) {
        globalx[i] = (Ii + Istart) * blockSize + i; // 计算每一行在全局矩阵中的索引
    }

    // 遍历每一列，插入矩阵值
    for (int i = rpt[Ii]; i < rpt[Ii + 1]; i++) {
        for (int j = 0; j < blockSize; j++) {
            globaly[j] = cpt[i] * blockSize + j; // 计算全局列索引
        }
        ierr = MatSetValues(A, blockSize, globalx, blockSize, globaly, valpt, INSERT_VALUES); // 设置矩阵值
        CHKERRQ(ierr); // 检查错误
        valpt += blockSize * blockSize; // 移动指针以处理下一块
    }
}

// 完成矩阵的装配，标志装配的开始和结束
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);

// 创建右手边向量 b，并设置其大小
ierr = VecCreate(PETSC_COMM_WORLD, &b);
CHKERRQ(ierr);
ierr = VecSetSizes(b, rowWidth, matrixDim); // 设置 b 的大小
CHKERRQ(ierr);
ierr = VecSetFromOptions(b); // 从命令行选项设置向量
CHKERRQ(ierr);
ierr = VecDuplicate(b, &x); // 创建与 b 相同大小的解向量 x
CHKERRQ(ierr);

// 设置 b 向量的值
int *bIndex = (int *)malloc(rowWidth * sizeof(int)); // 动态分配数组存储索引
for (i = 0; i < rowWidth; i++) {
    bIndex[i] = Istart * blockSize + i; // 计算索引
}
VecSetValues(b, rowWidth, bIndex, rhs, INSERT_VALUES); // 设置 b 向量的值
VecAssemblyBegin(b); // 开始装配 b 向量
VecAssemblyEnd(b);   // 结束装配 b 向量

// 创建用于进一步计算的 dBSR 矩阵对象
dBSRmat_ BMat; // 定义一个 dBSRmat 结构，用于存储矩阵信息
BMat.ROW = nBlockRows; // 设置矩阵行数
BMat.COL = nrow_global; // 设置矩阵列数
BMat.NNZ = nnz; // 设置非零元素个数
BMat.nb = nb; // 设置块大小
BMat.IA = rpt; // 设置行指针
BMat.JA = cpt; // 设置列指针
BMat.val = val; // 设置矩阵值

Mat App; // App 用于存储子矩阵
// Mat Ass; // Ass 是另一个子矩阵，但当前没有使用

finish = MPI_Wtime(); // 结束时间
Matrix_Conversion_Time += finish - start; //! 记录矩阵转换的时间（全局变量）

优化后的代码如下

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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
////////////////////////////////////////////////////////////////////
// Modifiable Area  |
// Modifiable Area  v
/////////////////////////////////////////////////////////////////////
start = MPI_Wtime();
decoup(val, rhs, rpt, cpt, nb, nBlockRows, DecoupType, is_thermal);

//----------------------------------------------------------------------
// Initializing MAT and VEC

Vec x, b; /* approx solution, RHS, exact solution */
Mat A;    /* linear system matrix */
KSP ksp;  /* linear solver context */
PC pc;
// PetscReal      norm;     /* norm of solution error */
PetscInt Ii, iters;
PetscErrorCode ierr;

// 计算 nDCount 和 nNDCount
for (int i = 0; i < nBlockRows; i++)
{
    nDCount[i] = 0;
    for (int j = rpt[i]; j < rpt[i + 1]; j++)
    {
        if (cpt[j] >= Istart && cpt[j] <= Iend)
        {
            nDCount[i]++;
        }
    }
}

for (int i = 0; i < nBlockRows; i++)
{
    nNDCount[i] = rpt[i + 1] - rpt[i] - nDCount[i];
}

// 创建矩阵 A
ierr = MatCreateBAIJ(PETSC_COMM_WORLD, blockSize, rowWidth, rowWidth, matrixDim, matrixDim, 0, nDCount, 0, nNDCount, &A);
CHKERRQ(ierr);
ierr = MatSetFromOptions(A);
CHKERRQ(ierr);
ierr = MatSetUp(A);
CHKERRQ(ierr);

// 使用 PETSc 进行内存分配
double *valpt = val;
int *globalx;
ierr = PetscMalloc1(blockSize, &globalx);
CHKERRQ(ierr);

int *globaly;
ierr = PetscMalloc1(blockSize, &globaly);
CHKERRQ(ierr);

for (Ii = 0; Ii < nBlockRows; Ii++)
{
    // 预计算 rowBase 以减少重复计算
    int rowBase = (Ii + Istart) * blockSize;

    // 填充 globalx
    for (i = 0; i < blockSize; i++)
    {
        globalx[i] = rowBase + i;
    }

    // 填充 globaly 和设置矩阵值
    for (int i = rpt[Ii]; i < rpt[Ii + 1]; i++)
    {
        for (int j = 0; j < blockSize; j++)
        {
            globaly[j] = cpt[i] * blockSize + j;
        }
        ierr = MatSetValues(A, blockSize, globalx, blockSize, globaly, valpt, INSERT_VALUES);
        CHKERRQ(ierr);
        valpt += blockSize * blockSize; // 更新 valpt
    }
}

// 矩阵组装
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);

// 创建和设置向量 b
ierr = VecCreate(PETSC_COMM_WORLD, &b);
CHKERRQ(ierr);
ierr = VecSetSizes(b, rowWidth, matrixDim);
CHKERRQ(ierr);
ierr = VecSetFromOptions(b);
CHKERRQ(ierr);
ierr = VecDuplicate(b, &x);
CHKERRQ(ierr);

// 使用 PETSc 进行内存分配
int *bIndex;
ierr = PetscMalloc1(rowWidth, &bIndex);
CHKERRQ(ierr);

// 填充 bIndex
for (i = 0; i < rowWidth; i++)
{
    bIndex[i] = Istart * blockSize + i;
}

// 设置向量 b 的值
VecSetValues(b, rowWidth, bIndex, rhs, INSERT_VALUES);
VecAssemblyBegin(b);
VecAssemblyEnd(b);

// 定义和初始化 dBSRmat_ 结构体
dBSRmat_ BMat;
BMat.ROW = nBlockRows;
BMat.COL = nrow_global;
BMat.NNZ = nnz;
BMat.nb = nb;
BMat.IA = rpt;
BMat.JA = cpt;
BMat.val = val;
Mat App;
// 释放内存
ierr = PetscFree(globalx);
CHKERRQ(ierr);
ierr = PetscFree(globaly);
CHKERRQ(ierr);
ierr = PetscFree(bIndex);
CHKERRQ(ierr);

finish = MPI_Wtime();
Matrix_Conversion_Time += finish - start; //! matrix conversion time (Global variables)
    
////////////////////////////////////////////////////////////////////
// Modifiable Area  ^
// Modifiable Area  |
/////////////////////////////////////////////////////////////////////

  1. 移除了重复内存分配的部分：

    • 在原代码中，对于 globalx 和 globaly 这两个指针的内存分配是通过 malloc 实现的，而在优化后代码中，这些内存是使用 PetscMalloc1 分配的，这样可以更好地整合 PETSc 自己的内存管理，提高内存分配和释放的性能。

    • 原代码中使用了 int *globalx = (int *)malloc(blockSize * sizeof(int)) 和 int *globaly = (int *)malloc(blockSize * sizeof(int))，而在第二段代码中替换为 ierr = PetscMalloc1(blockSize, &globalx) 和 ierr = PetscMalloc1(blockSize, &globaly)，并且这些内存通过 PetscFree 进行释放，确保与 PETSc 的内存管理机制兼容。

  2. 重复计算的减少：

    • 在优化后的代码中，循环中预计算了 rowBase 变量，这减少了在每次循环中重复进行的乘法运算：

      1
      int rowBase = (Ii + Istart) * blockSize;
      • 这一优化避免了在循环内多次计算 rowBase，提高了性能，特别是在循环次数较多时。

  3. 代码的组织更加结构化：

    • 优化后的代码将 globalx 和 globaly 的赋值与矩阵 A 的值的填充更加结构化，使代码逻辑清晰，更易于维护和理解。

  4. 内存释放管理的改进：

    • 优化后的代码使用 PetscFree 函数显式释放内存，这样可以与 PETSc 的内存管理体系保持一致，减少内存泄漏的风险。

    • 例如：

      1
      2
      3
      4
      5
      6
      ierr = PetscFree(globalx);
      CHKERRQ(ierr);
      ierr = PetscFree(globaly);
      CHKERRQ(ierr);
      ierr = PetscFree(bIndex);
      CHKERRQ(ierr);

      而在原代码中则没有对 globalx 和 globaly 进行显式的释放。

  5. 代码的错误检查：

    • 优化后的代码在执行内存分配、矩阵和向量创建等操作后，通过 CHKERRQ(ierr) 检查错误，提高了代码的鲁棒性，使得在出现错误时更容易定位问题。

可以继续尝试的可能

矩阵和向量的重用：对于多次重复使用的矩阵 A 和向量 b，可以考虑将它们保留并复用，而不是每次都重新创建，这样可以减少初始化的开销。

并行化内存初始化：在对 nDCount 和 nNDCount 进行初始化时，可以考虑并行化这部分代码，利用 OpenMP 或 MPI 实现并行加速，从而加快初始化速度。

进一步的内存对齐和缓存优化：考虑到 globalx 和 globaly 数组可以频繁访问，可以使用内存对齐来优化它们在 CPU 缓存中的使用，利用指令集提供的对齐内存分配，例如使用 _mm_malloc()。

批量设置矩阵值：使用 PETSc 的批量操作 MatSetValuesBlocked 来设置矩阵值，可以减少函数调用次数，提高整体性能。

同时PETSc中集成了与GPU加速线性求解器的接口

在编译PETSc时，使用 ./configure --with-cuda=1 --with-cudac=/usr/local/cuda/bin/nvcc，这样可以让PETSc支持CUDA，并且与GPU进行集成。

代码修改：

  • 在PETSc代码中，将向量类型设置为CUDA：

    1
    2
    cpp复制代码VecCreate(PETSC_COMM_WORLD, &x);
    VecSetType(x, VECCUDA);

  • 设置矩阵类型为CUDA：

    1
    2
    cpp复制代码MatCreate(PETSC_COMM_WORLD, &A);
    MatSetType(A, MATAIJCUDA);

  • 这些修改将确保PETSc将矩阵和向量的操作放在GPU上执行。

预处理器和求解器设置：

  • 设置求解器使用GPU加速：

    1
    2
    3
    4
    5
    bash


    复制代码
    -pc_type hypre -pc_hypre_type boomeramg -pc_hypre_boomeramg_relax_type_all symmetric-SOR/Jacobi -pc_hypre_boomeramg_interp_type ext+i -ksp_type gmres

  • 这些选项会使用Hypre的BoomerAMG预处理器，并在GPU上进行加速。

开启cuda加速

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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
////////////////////////////////////////////////////////////////////
// Modifiable Area  |
// Modifiable Area  v
/////////////////////////////////////////////////////////////////////
start = MPI_Wtime();
decoup(val, rhs, rpt, cpt, nb, nBlockRows, DecoupType, is_thermal);

//----------------------------------------------------------------------
// Initializing MAT and VEC

Vec x, b; /* approx solution, RHS, exact solution */
Mat A;    /* linear system matrix */
KSP ksp;  /* linear solver context */
PC pc;
// PetscReal      norm;     /* norm of solution error */
PetscInt Ii, iters;
PetscErrorCode ierr;

// 计算 nDCount 和 nNDCount
for (int i = 0; i < nBlockRows; i++)
{
    nDCount[i] = 0;
    for (int j = rpt[i]; j < rpt[i + 1]; j++)
    {
        if (cpt[j] >= Istart && cpt[j] <= Iend)
        {
            nDCount[i]++;
        }
    }
}

for (int i = 0; i < nBlockRows; i++)
{
    nNDCount[i] = rpt[i + 1] - rpt[i] - nDCount[i];
}

// 创建矩阵 A，并启用 GPU 支持
ierr = MatCreateBAIJ(PETSC_COMM_WORLD, blockSize, rowWidth, rowWidth, matrixDim, matrixDim, 0, nDCount, 0, nNDCount, &A);
CHKERRQ(ierr);

// 使用 GPU 进行计算，设置矩阵类型为 cuSPARSE
ierr = MatSetType(A, MATAIJCUSPARSE); CHKERRQ(ierr);
ierr = MatSetOption(A, MAT_USE_GPU, PETSC_TRUE); CHKERRQ(ierr);

ierr = MatSetFromOptions(A); CHKERRQ(ierr);
ierr = MatSetUp(A); CHKERRQ(ierr);

// 使用 PETSc 进行内存分配
double *valpt = val;
int *globalx;
ierr = PetscMalloc1(blockSize, &globalx);
CHKERRQ(ierr);

int *globaly;
ierr = PetscMalloc1(blockSize, &globaly);
CHKERRQ(ierr);

for (Ii = 0; Ii < nBlockRows; Ii++)
{
    // 预计算 rowBase 以减少重复计算
    int rowBase = (Ii + Istart) * blockSize;

    // 填充 globalx
    for (int i = 0; i < blockSize; i++)
    {
        globalx[i] = rowBase + i;
    }

    // 填充 globaly 和设置矩阵值
    for (int i = rpt[Ii]; i < rpt[Ii + 1]; i++)
    {
        for (int j = 0; j < blockSize; j++)
        {
            globaly[j] = cpt[i] * blockSize + j;
        }
        ierr = MatSetValues(A, blockSize, globalx, blockSize, globaly, valpt, INSERT_VALUES);
        CHKERRQ(ierr);
        valpt += blockSize * blockSize; // 更新 valpt
    }
}

// 矩阵装配
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

// 创建和设置向量 b，启用 GPU 支持
ierr = VecCreate(PETSC_COMM_WORLD, &b); CHKERRQ(ierr);
ierr = VecSetSizes(b, rowWidth, matrixDim); CHKERRQ(ierr);
ierr = VecSetType(b, VECCUDA); CHKERRQ(ierr); // 设置向量类型为 CUDA
ierr = VecSetFromOptions(b); CHKERRQ(ierr);
ierr = VecDuplicate(b, &x); CHKERRQ(ierr);
ierr = VecSetType(x, VECCUDA); CHKERRQ(ierr); // 设置 x 的类型为 CUDA

// 使用 PETSc 进行内存分配
int *bIndex;
ierr = PetscMalloc1(rowWidth, &bIndex);
CHKERRQ(ierr);

// 填充 bIndex
for (int i = 0; i < rowWidth; i++)
{
    bIndex[i] = Istart * blockSize + i;
}

// 设置向量 b 的值
VecSetValues(b, rowWidth, bIndex, rhs, INSERT_VALUES);
VecAssemblyBegin(b);
VecAssemblyEnd(b);

// 定义和初始化 dBSRmat_ 结构体
dBSRmat_ BMat;
BMat.ROW = nBlockRows;
BMat.COL = nrow_global;
BMat.NNZ = nnz;
BMat.nb = nb;
BMat.IA = rpt;
BMat.JA = cpt;
BMat.val = val;
Mat App;

// 释放内存
ierr = PetscFree(globalx);
CHKERRQ(ierr);
ierr = PetscFree(globaly);
CHKERRQ(ierr);
ierr = PetscFree(bIndex);
CHKERRQ(ierr);

finish = MPI_Wtime();
Matrix_Conversion_Time += finish - start; //! matrix conversion time (Global variables)

////////////////////////////////////////////////////////////////////
// Modifiable Area  ^
// Modifiable Area  |
/////////////////////////////////////////////////////////////////////

在查阅PETSc文档后，确定了Mat, Vector, Ksp等PETSc变量可以被重复利用的前提下，可以将其改变为静态类型变量，得到进一步的性能提升，但是可惜的是并没有观察到有效的变化，估计是编译器做了一些神奇的优化掩盖掉了。（o3魅力时刻了说是）

-------------本文结束感谢您的阅读-------------