一、问题描述
矩阵LU分解是一种重要的数值线性代数运算方法，广泛应用于科学计算、工程问题和数据分析等领域。LU分解的背景可以追溯到高斯消元法，它是一种常用的线性方程组求解方法。高斯消元法将线性方程组的系数矩阵通过行变换转化为上三角形式，然后通过回代求解得到方程组的解。然而，高斯消元法的计算量较大，尤其对于大规模的线性方程组，计算效率较低。
为了提高计算效率，LU分解应运而生。它将线性方程组的系数矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A = LU。
LU分解的优势在于一次分解后可以重复使用，从而在解决多个线性方程组时提高计算效率。此外，LU分解还可以简化解线性方程组、矩阵求逆、特征值计算和奇异值分解等问题的求解过程。
在实际应用中，LU分解通常与其他数值算法结合使用。例如，在求解线性方程组时，可以通过LU分解将系数矩阵分解为L和U，并利用这个分解形式来快速求解多个右端向量的解。对于矩阵求逆，可以通过LU分解将矩阵分解为L和U，并利用这个分解形式来快速求解多个右端向量的解。对于矩阵求逆，可以通过LU分解将矩阵分解为L和U，并利用这个分解形式来快速计算矩阵的逆。
二、平台说明
1.硬件
使用的服务型号为：HUAWEI TaiShan 200服务器（型号2280）。TaiShan 200服务器是基于华为鲲鹏920处理器的数据中心服务器，是2U 2路机架服务器，支持2路处理器，处理器规格为32核2.6GHz。
以12块硬盘配置为例的外观图如图2-1所示。服务器整体外观图如图2-2所示。
以12盘配置为例的服务器的物理结构如图2-3所示，部件说明如表2-1。
服务器支持Hi1711 BMC插卡，可外出VGA、管理网口、串口、USB Type-C等管理接口，其逻辑结构如图2-4所示。
该服务器支持华为自研的、面向服务器领域的64 bits高性能多核鲲鹏920处理器，内部集成了DDR4、PCIe 4.0、GE等接口，提供完整的SOC功能。单台服务器支持2个处理器，单个处理器最大支持32个内核，能够最大限度地提高多线程应用的并发执行能力。支持多种灵活的硬盘配置方案，提供了弹性的、可扩展的存储容量空间，满足不同存储容量的需求和升级要求。支持灵活插卡，可提供多种以太网卡接口能力。最多可支持8个PCIe 4.0 x8的标准扩展槽位。
擅于面向互联网、分布式存储、云计算、大数据、企业业务等领域，具有高性能计算、大容量存储、低能耗、易管理、易部署等优点。
2.软件
操作系统内核信息：
操作系统版本信息：
查看CPU详细信息：
查看服务器内存：
三、算法流程图
四、代码中的关键并行性
1、MPI并行性：通过MPI_Init()和MPI_Comm_size()函数初始化进程，并使用MPI_Comm_rank()函数获取进程的排名。然后，在主循环中，使用MPI_Bcast()函数将主行元素广播给其他进程，以实现并行计算和通信。
2、OpenMP并行性：通过#pragma omp parallel和#pragma omp for指令实现了OpenMP并行化。在主循环中，使用#pragma omp parallel创建多个线程，然后使用#pragma omp for并行执行对于矩阵进行初等变换的二重循环。这样可以在每个线程上并行处理行操作，加速LU分解的计算过程。
3、NEON指令集并行性：在使用NEON指令集进行计算的部分，通过使用SIMD指令（如vld1q_f32、vmlsq_f32和vst1q_f32），可以在每个迭代步骤中同时处理多个浮点数值，从而实现向量化并行计算。
通过MPI并行性、OpenMP并行性和NEON指令集并行性的结合，该代码能够同时充分利用多个进程、多个线程和向量化处理多个数据，从而加速LU分解的计算过程。
五、测试方案
首先随机生成大小规模为500×500的矩阵，运行速度快数据量小便于测试，观察运行生成的两个下、上三角矩阵文件，发现数据结果符合LU分解结果，为对角线全1的下三角矩阵和对角线不全为1的上三角矩阵，选取首1个和末1个数值进行计算比对发现生成结果与计算结果相同，代码正确性得到验证。
其次随机生成2w*2w、3w*3w、3.3w*3.3w、5w*5w用于分别测试代码的性能，并进行并行性运算优化调增，分别在服务器和PC机进行运行观察。
六、实验结果
1.服务器
采取MPI、OpenMP和NEON指令集的结合，实现对于超大规模矩阵的LU分解，详细服务器端运行代码见附录1-服务器端代码
首先使用了小规模数据（随机生成500×500矩阵）进行代码正确性验证，结果正常输出两个矩阵文件L_0.out和_0.out，打开后显示内容分别为下、上三角矩阵，表明结果正确，代码运行正常。
其次增大至20000乘20000矩阵进行验证，结果正常运行和输出，矩阵变换耗时131.094770 秒；矩阵分离耗时131.862488 秒；结果输出耗时357.272186 秒，并正常输出两个下、上三角矩阵文件。
继续增大矩阵至50000×50000时，程序报错：内存分配失败（为初始化矩阵时自行编写增加的判断输出语句）。原因为无法分配5w×5w大小的连续的float型数据的空间用以初始化矩阵，因此程序终止。
降低矩阵至40000×40000，程序恢复正常运行。（此时由于生成文件耗时过长，预计为矩阵变换所需时间的十倍以上，故将生成文件部分暂时进行注释，仅想先观察程序算法部分是否能够正常运行，因此结果输出和矩阵分离显示时间几乎相同）
补40000×40000规模矩阵LU分解程序运行图，启用结果模块以验证生成结果是否正确，生成程序正常运行。矩阵变换用时2068s，矩阵分离用时2079s，结果输出用时4622s。
2.PC
PC机由于并非ARM指令集，因此本机代码中不包含使用NEON技术进行向量化并行优化计算的部分，而仅用OMP与MPI进行并行加速计算，其余地方均无改动，详细代码见附录2-PC端代码。
首先进行30000×30000规模矩阵的LU分解计算，正常输出结果，矩阵变换和分离用时分别为：77s和97s，生成两个下、上三角矩阵文件L_0.out和U_0.out。
其次增大矩阵规模至33000×33000，同样可以正常输出结果，矩阵变换和分离用时分别为：95s和120s，结果输出时间为2460s，生成两个下、上三角矩阵文件L_0.out和U_0.out。
继续扩大矩阵规模是，程序在运行过程中中断，在变换完成后无法进行LU的正常分解，经查阅异常信息资料和断点位置调试，以及设置测试输出部分矩阵信息代码，发现由于本机内存空间不足，无法继续连续分配空间而导致内存访问出界，无法继续实现更大规模矩阵的LU分解计算。
七、实验总结
尝试解决超大矩阵（规模为五万乘五万）的LU分解问题是一项非常具有挑战性的任务，包括：数学算法、算法代码优化、硬件资源的利用等。在这个过程中锻炼了我分析和解决问题的能力，同时让我回顾和熟悉了OpenMP、MPI和NEON技术以及华为泰山服务器平台的使用方式，我也深刻地领悟到了通过结合OpenMP、MPI和NEON技术，能搞非常明显的提高计算性能。
首先是计算时间的降低，即性能提升：使用并行计算技术可以显著提高LU分解的性能，特别是在处理大规模矩阵时。OpenMP和MPI可以实现多线程和多进程的并行计算，而NEON指令集可以在单个线程中利用SIMD（单指令多数据）并行性。这种组合可以有效地利用计算资源，并加速矩阵的LU分解过程。
其次就是内存空间的有限，要使用任务划分的办法：在使用OpenMP进行并行计算时，可以将矩阵按照适当的块大小进行划分，并将每个线程分配到不同的块上进行计算。这样可以实现数据的并行处理，减少了线程之间的数据依赖，提高了计算效率，降低了单线程内的内存空间资源要求。同时可以考虑进行算法的优化，即本代码中将最开始使用的二维矩阵存储数据改为使用一位数组表示二维矩阵的方法，能够很有限的节省和高效率空间资源，从而部分解决内存空间不足的问题。（仅解决至4w×4w规模的矩阵，目前仍未解决5w×5w矩阵LU分解）
同时要考虑到进程通信的问题：即使用MPI实现多进程并行计算时，需要注意进程之间的通信和数据同步。在LU分解的过程中，主行元素需要广播到所有进程中的所有线程，以保证计算的一致性。本代码中使用了MPI提供的广播操作（MPI_Bcast）等通信机制，用于实现数据的共享和同步。同时考虑到多进程并行时，广播机制可能出现多个线程重复进行的情况，故采用#pragma omp single实现单线程广播至其他线程。
最后由于在服务器进行运行，因此还可以使用NEON优化技术：通过使用NEON指令集，可以在单个线程中对矩阵的行进行并行计算。NEON提供了一组SIMD指令，可以同时对多个数据执行相同的操作，提高了计算效率。在LU分解过程中，可以利用NEON指令对矩阵的行进行并行计算，减少了循环迭代的次数，提高了性能。PC机无法使用NEON指令集（其为ARN指令集下的向量并行优化计算技术）同时服务器硬件资源优于PC机，于是实现了服务器中实现40000×40000规模矩阵的LU分解，而PC机仅最高实现33000×33000矩阵LU分解。
未采取并行优化技术时，程序容易内存爆炸，如下图7-1；且用时过长，如下图7-2（未采用NEON技术优化时计算2w×2w矩阵规模）；优化后如图7-3。
综上所述，结合OpenMP、MPI和NEON技术可以有效地解决超大矩阵的LU分解问题，并提高计算性能。但是实现并行算法时需要注意考虑到考虑数据划分、通信、同步、采用线程数等问题，同时也需要根据具体硬件和任务的特点进行优化，以获得最佳的性能。
当矩阵大小较小且采用的分块过少时，如500×500大小，仅需分块数为一块时，运行线程数不能大于1，否则多余线程将没有需要处理的分块矩阵对象，而导致内存访问越界（EXIT CODE 139）。如下图7-4所示。
八、附录
1.服务端代码
1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <time.h>
4. #include <mpi.h>
5. #include <omp.h>
6. #include <arm_neon.h>
7.
8. #define BLOCK_SIZE 1000
9.
10. // 方阵
11. float* A;
12. int N = 20000;
13. float* L, * U;
14.
15. // 初始化矩阵
16. void initialize_matrix(int myid, int n_threads) {
17. int local_N = N / n_threads;
18. A = (float*)malloc(local_N * N * sizeof(float));
19. if (A == NULL) {
20. perror("内存分配失败");
21. exit(-1);
22. }
23.
24. srand(time(NULL) + myid);
25. #pragma omp parallel for
26. for (int i = 0; i < local_N * N; i++) {
27. A[i] = (float)rand() / (float)(RAND_MAX / 100);
28. }
29. }
30.
31. // 释放矩阵内存
32. void free_matrix(float* matrix) {
33. if (matrix != NULL) {
34. free(matrix);
35. }
36. }
37.
38. // 从计算好的A中分离、整理，得到最终L、U结果
39. void getLU(int myid, int n_threads) {
40. int local_N = N / n_threads;
41. // 生成和初始化 L、U 矩阵
42. L = (float*)calloc(local_N * N, sizeof(float));
43. U = (float*)calloc(local_N * N, sizeof(float));
44.
45. // 从计算好的 A 中，分离得到 L、U
46. for (int i = 0; i < local_N; i++) {
47. for (int j = 0; j < N; j++) {
48. if (i + myid * local_N < j) {
49. U[i * N + j] = A[i * N + j];
50. L[i * N + j] = 0;
51. }
52. else if (i + myid * local_N == j) {
53. U[i * N + j] = A[i * N + j];
54. L[i * N + j] = 1;
55. }
56. else {
57. U[i * N + j] = 0;
58. L[i * N + j] = A[i * N + j];
59. }
60. }
61. }
62. }
63.
64. int main(int argc, char* argv[]) {
65. int myid = 0;
66. int n_threads = 10;
67. double te, ts;
68.
69. MPI_Init(&argc, &argv);
70. MPI_Comm_rank(MPI_COMM_WORLD, &myid);
71. MPI_Comm_size(MPI_COMM_WORLD, &n_threads);
72.
73. // 初始化矩阵A
74. if (argc >= 2) {
75. // 通过运行参数，指定输入文件
76. //init(argv[1]);
77. exit(0);
78. }
79. else {
80. // 或者：随机数初始化
81. initialize_matrix(myid, n_threads);
82. }
83.
84. // 计时开始
85. ts = MPI_Wtime();
86.
87. // 分块处理
88. for (int i = 0; i < N; i += BLOCK_SIZE) {
89. int row_end = i + BLOCK_SIZE;
90. if (row_end > N) {
91. row_end = N;
92. }
93.
94. // 逐行作为主行元素，进行初等行变换
95. #pragma omp parallel shared(A)
96. {
97. for (int j = i; j < row_end; j++) {
98. // 同步当前的主行元素到所有线程
99. #pragma omp barrier
100. #pragma omp single
101. {
102. MPI_Bcast(&A[j * N], N, MPI_FLOAT, j % n_threads, MPI_COMM_WORLD);
103. }
104.
105. // 各线程将自己所负责的行完成行初等变换
106. #pragma omp for
107. for (int k = j + 1; k < N; k++) {
108. A[j * N + k] = A[j * N + k] / A[j * N + j]; // 初等行变换的第一步，第一个元素作除法（也相当于计算出了 L[j][j]）
109. }
110.
111. #pragma omp for
112. for (int k = j + 1; k < N; k++) {
113. float32x4_t ajk_vec = vdupq_n_f32(A[j * N + k]);
114. for (int l = j + 1; l < row_end; l += 4) {
115. float32x4_t alj_vec = vld1q_f32(&A[l * N + j]);
116. float32x4_t alk_vec = vld1q_f32(&A[l * N + k]);
117. alk_vec = vmlsq_f32(alk_vec, alj_vec, ajk_vec);
118. vst1q_f32(&A[l * N + k], alk_vec);// 初等行变化第二步，该行其他元素进行初等行变化
119. }
120. }
121. }
122. }
123. }
124. // 计时结束
125. te = MPI_Wtime();
126. if (myid == 0) {
127. printf("矩阵变换完成，耗时：%f 秒。\n", te - ts);
128. }
129.
130. // 从计算好的 A 中，分离得到 L、U
131. getLU(myid, n_threads);
132.
133. // 计时结束
134. te = MPI_Wtime();
135. if (myid == 0) {
136. printf("矩阵分离完成，耗时：%f 秒。\n", te - ts);
137. }
138.
139. //输出结果
140. if (argc >= 3) {
141. // 输出到自定文件
142. // output(argv[2]);
143. }
144. else {
145. // 输出到默认文件
146. // 输出L、U矩阵到文件中
147. char filenameL[20], filenameU[20];
148. sprintf(filenameL, "L_%d.out", myid);
149. sprintf(filenameU, "U_%d.out", myid);
150. FILE* fpL = fopen(filenameL, "w");
151. FILE* fpU = fopen(filenameU, "w");
152. for (int i = 0; i < N / n_threads * N; i++) {
153. fprintf(fpL, "%f ", L[i]);
154. fprintf(fpU, "%f ", U[i]);
155. if (i % N == 0 && i > 0) {
156. fprintf(fpL, "\n");
157. fprintf(fpU, "\n");
158. }
159. }
160. fclose(fpL);
161. fclose(fpU);
162. }
163.
164. // 计时结束
165. te = MPI_Wtime();
166. if (myid == 0) {
167. printf("结果输出完成，耗时：%f 秒。\n", te - ts);
168. }
169.
170. // 释放内存
171. free_matrix(A);
172. free_matrix(L);
173. free_matrix(U);
174.
175. MPI_Finalize();
176.
177. return 0;
178. }
2.PC端代码
1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <time.h>
4. #include <mpi.h>
5. #include <omp.h>
6.
7. #define BLOCK_SIZE 100
8.
9. // 方阵
10. float* A;
11. int N = 34000;
12. float* L, * U;
13.
14. // 初始化矩阵
15. void initialize_matrix(int myid, int n_threads) {
16. int local_N = N / n_threads;
17. A = (float*)malloc(local_N * N * sizeof(float));
18. if (A == NULL) {
19. perror("内存分配失败");
20. exit(-1);
21. }
22.
23. srand(time(NULL) + myid);
24. #pragma omp parallel for
25. for (int i = 0; i < local_N * N; i++) {
26. A[i] = (float)rand() / (float)(RAND_MAX / 100);
27. }
28. }
29.
30. // 释放矩阵内存
31. void free_matrix(float* matrix) {
32. if (matrix != NULL) {
33. free(matrix);
34. }
35. }
36.
37. // 从计算好的A中分离、整理，得到最终L、U结果
38. void getLU(int myid, int n_threads) {
39. int local_N = N / n_threads;
40. //printf("local_N: %d, N: %d\n", local_N, N);
41.
42. // 生成和初始化 L、U 矩阵
43. L = (float*)calloc(local_N * N, sizeof(float));
44. U = (float*)calloc(local_N * N, sizeof(float));
45.
46. // 从计算好的 A 中，分离得到 L、U
47. for (int i = 0; i < local_N; i++) {
48. for (int j = 0; j < N; j++) {
49. //printf("i: %d, j: %d, A: %p, U: %p\n", i, j, A, U);
50.
51. if (i + myid * local_N < j) {
52. U[i * N + j] = A[i * N + j];
53. L[i * N + j] = 0;
54. }
55. else if (i + myid * local_N == j) {
56. U[i * N + j] = A[i * N + j];
57. L[i * N + j] = 1;
58. }
59. else {
60. U[i * N + j] = 0;
61. L[i * N + j] = A[i * N + j];
62. }
63. }
64. }
65. }
66.
67. int main(int argc, char* argv[]) {
68. int myid = 0;
69. int n_threads = 10;
70. double te, ts;
71.
72. MPI_Init(&argc, &argv);
73. MPI_Comm_rank(MPI_COMM_WORLD, &myid);
74. MPI_Comm_size(MPI_COMM_WORLD, &n_threads);
75.
76. // 初始化矩阵A
77. if (argc >= 2) {
78. // 通过运行参数，指定输入文件
79. //init(argv[1]);
80. exit(0);
81. }
82. else {
83. // 或者：随机数初始化
84. initialize_matrix(myid, n_threads);
85. }
86.
87. // 计时开始
88. ts = MPI_Wtime();
89.
90. // 分块处理
91. for (int i = 0; i < N; i += BLOCK_SIZE) {
92. int row_end = i + BLOCK_SIZE;
93. if (row_end > N) {
94. row_end = N;
95. }
96.
97. // 逐行作为主行元素，进行初等行变换
98. #pragma omp parallel shared(A)
99. {
100. for (int j = i; j < row_end; j++) {
101. // 同步当前的主行元素到所有线程
102. #pragma omp barrier
103. #pragma omp single
104. {
105. MPI_Bcast(&A[j * N], N, MPI_FLOAT, j % n_threads, MPI_COMM_WORLD);
106. }
107.
108. // 各线程将自己所负责的行完成行初等变换
109. #pragma omp for
110. for (int k = j + 1; k < N; k++) {
111. A[j * N + k] = A[j * N + k] / A[j * N + j]; // 初等行变换的第一步，第一个元素作除法（也相当于计算出了 L[j][j]）
112. }
113.
114. #pragma omp for collapse(2)
115. for (int k = j + 1; k < N; k++) {
116. for (int l = j + 1; l < row_end; l++) {
117. A[l * N + k] = A[l * N + k] - A[l * N + j] * A[j * N + k]; // 初等行变化第二步，该行其他元素进行初等行变化
118. }
119. }
120. }
121. }
122. }
123. // 计时结束
124. te = MPI_Wtime();
125. if (myid == 0) {
126. printf("矩阵变换完成，耗时：%f 秒。\n", te - ts);
127. }
128.
129. // 从计算好的 A 中，分离得到 L、U
130. getLU(myid, n_threads);
131.
132. // 计时结束
133. te = MPI_Wtime();
134. if (myid == 0) {
135. printf("矩阵分离完成，耗时：%f 秒。\n", te - ts);
136. }
137.
138. //输出结果
139. if (argc >= 3) {
140. // 输出到自定文件
141. // output(argv[2]);
142. }
143. else {
144. // 输出到默认文件
145. // 输出L、U矩阵到文件中
146. char filenameL[20], filenameU[20];
147. sprintf(filenameL, "L_%d.out", myid);
148. sprintf(filenameU, "U_%d.out", myid);
149. FILE* fpL = fopen(filenameL, "w");
150. FILE* fpU = fopen(filenameU, "w");
151. for (int i = 0; i < N / n_threads * N; i++) {
152. fprintf(fpL, "%f ", L[i]);
153. fprintf(fpU, "%f ", U[i]);
154. if (i % N == 0 && i > 0) {
155. fprintf(fpL, "\n");
156. fprintf(fpU, "\n");
157. }
158. }
159. fclose(fpL);
160. fclose(fpU);
161. }
162.
163. // 计时结束
164. te = MPI_Wtime();
165. if (myid == 0) {
166. printf("结果输出完成，耗时：%f 秒。\n", te - ts);
167. }
168.
169. // 释放内存
170. free_matrix(A);
171. free_matrix(L);
172. free_matrix(U);
173.
174. MPI_Finalize();
175.
176. return 0;
177. }
836
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
