本文写于ASC24结束后一天，可能存在诸多错误，能力有限。同时这是我第一次接触MPI程序和HPC，还请谅解，之前学习的领域有重合但不多，非常欢迎大家指出错误斧正。
OpenCAEPoro的背景我想不需要赘述。开门见山的说，大家都发现了dgemm函数是一个热点，如果优化了这个函数，加速效果会非常明显，这也是我们所优化的地方。此外，CAE对左右邻居查找的算法时间复杂度特别高，基本上都是O(N)或O(N^2)的时间复杂度，这也是我们所优化的重点之一。最后PETSc有非常多可以优化调节的点，通过优化内存访问，decoup函数优化，我们可以达到更快的速度。同时合理利用异步计算，可以使程序同时榨干CPU运算性能和远端访问性能。
接下来，我们给出我们在决赛中的成绩作为参考，后续我会将初赛提交的文件上传到Github：
以下是初赛的加速比：
dgemm算法
首先是dgemm算法，ASC非官方交流群有大佬给出了应该是更优解，是Intel黑科技，我们不敢说我的代码更快，直达连接：https://www.intel.com/content/www/us/en/developer/articles/technical/onemkl-improved-small-matrix-performance-using-just-in-time-jit-code.html。相比之下我们的利用block缓存和cpp temple编译期间计算的方法就显得有些无力。这里可能会有人好奇，为什么不使用MKL的dgemm算法呢？经过我们的测试，发现MKL在对小矩阵中性能不佳，CAE中testcase1所计算的矩阵是相当小的，20x20以内（比赛中的case不得而知）。同时，我们还给出我们参考的代码，是来自于英伟达工程师的小教程，简直就是HPC入门最好的指南：https://github.com/yzhaiustc/Optimizing-DGEMM-on-Intel-CPUs-with-AVX512F
在我们的这份代码，小矩阵乘法相比作者的原版，能提供约30%的性能提升（因为部分缩放操作不会被编译）。
#define A(i, j) A[(i) * LDA + (j)]
#define B(i, j) B[(i) * LDB + (j)]
#define C(i, j) C[(i) * LDC + (j)]
#include "immintrin.h"
#include <stdint.h>
#include <cstring>
void scale_c_k4(double *C, int M, int N, int LDC, double scalar)
{
int i, j;
for (i = 0; i < M; i++)
{
for (j = 0; j < N; j++)
{
C(i, j) *= scalar;
}
}
}
void mydgemm_cpu_opt_k4(int M, int N, int K, double alpha, double *A, int LDA, double *B, int LDB, double beta, double *C, int LDC)
{
int i, j, k;
if (beta != 1.0)
scale_c_k4(C, M, N, LDC, beta);
for (i = 0; i < M; i++)
{
for (j = 0; j < N; j++)
{
double tmp = C(i, j);
for (k = 0; k < K; k++)
{
tmp += alpha * A(i, k) * B(k, j);
}
C(i, j) = tmp;
}
}
}
template <typename BETA, typename ALPHA>
void mydgemm_cpu_v4(int M, int N, int K, ALPHA alpha, double *A, int LDA, double *B, int LDB, BETA beta, double *C, int LDC)
{
int i, j, k;
if constexpr (beta == 0) // optimize for beta == 0
memset(C, 0, sizeof(double) * M * N);
else
{
if constexpr (beta != 1.0)
scale_c_k4(C, M, N, LDC, beta);
}
int M4 = M & -4, N4 = N & -4;
for (i = 0; i < M4; i += 4)
{
for (j = 0; j < N4; j += 4)
{
double c00 = C(i, j);
double c01 = C(i, j + 1);
double c02 = C(i, j + 2);
double c03 = C(i, j + 3);
double c10 = C(i + 1, j);
double c11 = C(i + 1, j + 1);
double c12 = C(i + 1, j + 2);
double c13 = C(i + 1, j + 3);
double c20 = C(i + 2, j);
double c21 = C(i + 2, j + 1);
double c22 = C(i + 2, j + 2);
double c23 = C(i + 2, j + 3);
double c30 = C(i + 3, j);
double c31 = C(i + 3, j + 1);
double c32 = C(i + 3, j + 2);
double c33 = C(i + 3, j + 3);
for (k = 0; k < K; k++)
{
double a0, a1, a2, a3;
if constexpr (alpha != 1.0)
{
a0 = alpha * A(i, k);
a1 = alpha * A(i + 1, k);
a2 = alpha * A(i + 2, k);
a3 = alpha * A(i + 3, k);
}
else
{
a0 = A(i, k);
a1 = A(i + 1, k);
a2 = A(i + 2, k);
a3 = A(i + 3, k);
}
double b0 = B(k, j);
double b1 = B(k, j + 1);
double b2 = B(k, j + 2);
double b3 = B(k, j + 3);
c00 += a0 * b0;
c01 += a0 * b1;
c02 += a0 * b2;
c03 += a0 * b3;
c10 += a1 * b0;
c11 += a1 * b1;
c12 += a1 * b2;
c13 += a1 * b3;
c20 += a2 * b0;
c21 += a2 * b1;
c22 += a2 * b2;
c23 += a2 * b3;
c30 += a3 * b0;
c31 += a3 * b1;
c32 += a3 * b2;
c33 += a3 * b3;
}
C(i, j) = c00;
C(i, j + 1) = c01;
C(i, j + 2) = c02;
C(i, j + 3) = c03;
C(i + 1, j) = c10;
C(i + 1, j + 1) = c11;
C(i + 1, j + 2) = c12;
C(i + 1, j + 3) = c13;
C(i + 2, j) = c20;
C(i + 2, j + 1) = c21;
C(i + 2, j + 2) = c22;
C(i + 2, j + 3) = c23;
C(i + 3, j) = c30;
C(i + 3, j + 1) = c31;
C(i + 3, j + 2) = c32;
C(i + 3, j + 3) = c33;
}
}
if (M4 == M && N4 == N)
return;
// boundary conditions
if (M4 != M)
mydgemm_cpu_opt_k4(M - M4, N, K, alpha, A + M4, LDA, B, LDB, 1.0, &C(M4, 0), LDC);
if (N4 != N)
mydgemm_cpu_opt_k4(M4, N - N4, K, alpha, A, LDA, &B(0, N4), LDB, 1.0, &C(0, N4), LDC);
}
查找左右邻接组分的优化
哎呀，接下来其实Vtune Profile就很难给出什么有用的信息了，好在程序代码不多，我们可以从最花费时间的地方优化起。Partion部分，我们人肉顶针出了OpenCAEPoro奇怪的地方，不同组分直接的连接是通过***遍历数组***实现的。这就导致查找左右孩子（不同组分的连接），在2x2x2的组分中，就有可能有24个链接，获得连接信息就要跑24次。有没有可能用哈希表啊？来自O(1)时间算法的自信。
然后，我们两位负责这个赛题的人，对着源代码肉眼检查了两周，终于发现可以简易优化的地方，位于Partition.cpp文件中的SetDistribution函数，以下仅为代码摘要，能这样修改的一共有两处，在这个部分，性能提高了7倍。
/// Get elements and their neighbors belonging to current process,
/// Also, process of them are needed.
/// Setup elements those will be sent to corresponding process
// target process, num of elements, num of edges, global index of elements, xadj, adjncy, adjproc
vector<vector<idx_t>> send_buffer;
vector<vector<idx_t>>& recv_buffer = elementCSR;
recv_buffer.push_back(vector<idx_t> {myrank, 0, 0});
// euphoria
unordered_map<idx_t, idx_t> ump_send_buffer;
// Calculate necessary memory first
for (idx_t i = 0; i < numElementLocal; i++) {
if (part[i] == myrank) {
recv_buffer[0][1] ++;
recv_buffer[0][2] += (xadj[i + 1] - xadj[i]);
}
else {
// euphoria
if(ump_send_buffer[part[i]]){
idx_t k = ump_send_buffer[part[i]] - 1;
send_buffer[k][1] ++;
send_buffer[k][2] += (xadj[i + 1] - xadj[i]);
}
else{
send_buffer.push_back(vector<idx_t>{part[i], 1, xadj[i + 1] - xadj[i]});
ump_send_buffer[part[i]] = send_buffer.size();
}
// idx_t k;
// for (k = 0; k < send_buffer.size(); k++) {
// if (part[i] == send_buffer[k][0]) {
// // existing process
// send_buffer[k][1] ++;
// send_buffer[k][2] += (xadj[i + 1] - xadj[i]);
// break;
// }
// }
// if (k == send_buffer.size()) {
// // new process
// send_buffer.push_back(vector<idx_t>{part[i], 1, xadj[i + 1] - xadj[i]});
// }
}
}
PETSc LinerSovler的优化
重写decoup函数是我们首先应用的优化，这是其中的一个子函数，先把重复释放的内存操作给修改了，改为静态内存，避免反复申请释放对内存池的压力，同时申请对齐的内存，提高访问效率。
void inverse1(double *A, int N)
{
static int *IPIV = static_cast<int *>(aligned_alloc(64, N * sizeof(int)));
int LWORK = N * N;
static int last_N = N;
static double *WORK = static_cast<double *>(aligned_alloc(64, LWORK * sizeof(double)));
int INFO;
if (last_N < N)
{
last_N = N;
std::free(IPIV);
std::free(WORK);
IPIV = static_cast<int *>(aligned_alloc(64, N * sizeof(int)));
WORK = static_cast<double *>(aligned_alloc(64, LWORK * sizeof(double)));
}
dgetrf_(&N, &N, A, &N, IPIV, &INFO);
if (INFO > 0)
{
printf("WARNING: The factorization has been completed, but the factor U is exactly singular");
}
else if (INFO < 0)
{
printf("WARNING: There is an illegal value");
}
dgetri_(&N, A, &N, IPIV, WORK, &LWORK, &INFO);
}
然后，我们想到同样的操作也也可以用在整个LS求解函数中，从微小的args，到相对较大的nNDCount等变量开始改起，这些操作约可以带来了10%的Object Time提升。在查阅文档的时候，还意外发现PETSc自己实现了一个内存池的操作，这个又可以带来10%的性能提升。
在查阅PETSc文档后，确定了Mat, Vector, Ksp等PETSc变量可以被重复利用的前提下，我们将其改变为静态类型变量，企图得到进一步的性能提升，但是可惜的是并没有观察到有效的变化，估计是编译器做了一些神奇的优化掩盖掉了。
在求解函数中，还可以发现存在MatAssemblyBegin和MatAssemblyBegin的函数，这些字眼有效暗示了我们这个是异步操作函数，查阅PETSc文档后得到了肯定。我们在其中塞入了一些计算过程，但是Object Time也没有肉眼可以观察到的变化，或许case太小了没法体现出来。同样，我们复用了这一思想，在结果拷贝的过程中使用C++ 20 feature特性，异步拷贝了数据（但这不在Object time内，只是为了比赛求解更快）。
接下来我要介绍一些我们失败的优化，尝试过更改矩阵搜引，矩阵类型，修改预求解器的方法，加速LS求解过程，亦或者使用Cuda加速线性求解过程（即使它并不在Object Time其中，但是在比赛时可以有效减少题目运行时间），都因为没有明显的性能提升或者求解错误而作罢。在这里面，我们花费了超过一半以上的备赛时间，却很遗憾没有收获。
修改汇总时的文件读写IO
重点：我这个实现的代码有Bug，思路没问题，但是最后主线程输出改炸了，会导致SUMMAY.out文件输出到虚空里面。这个修改是我自己完成的，在比赛时才发现，我本人难逃其咎，已经向队友祈求原谅了，很遗憾因为这个错误，导致几乎所有的case在决赛时都跑了两次，在此处公开道歉了QAQ。
在OpenCAEPoro中，汇总结果阶段会向磁盘中写入SUMMAY_proc_xxx和Fast_proc_xxx文件，再由主线程读取，整理出最终的数据，这样会导致可拓展性不高，对于上百核的并行，可能会造成性能损失，在研读代码后，我决定修改这一部分，不写入磁盘而通过MPI通讯汇总数据，最终得到10%的性能提升（线程越多性能提升越高，上千线程估计能跑到几百倍？）。
原来的过程像是这样：
我改成每个线程都发送给主线程（先发送长度，在发送内容），初次编写MPI不太会，应该用Gather方法会比较合适。
因为改炸了，实在不想说什么，居然能犯下这种错误，就说到这里吧。
其他七七八八的小优化
接下来的优化就开始玄学了，来自OI选手的快乐。
- 输入输出加速
取消掉与C语言库输入输出的同步，可以使得输入输出流达到Printf的水平（当然还有更快的）：
std::ios::sync_with_stdio(false);
std::cin.tie(0);
std::cout.tie(0);
给vector提前分配空间
代码中，不喜欢给vector模板reserve提前分配空间，但是有部分是可以被计算出来的，可以提前reserve分配好合适的空间。但是算不好会有反效果，默认下vector自动扩增的大小是当前大小的两倍。更好看的加速比
经过我们的多次测试，核心数越少，加速比越高（通讯的开销大了），这个代码在单机48核心的加速比约为1.68。手动展开
我们把自己修改的部分，手动展开了一部分循环。从汇编上对比手动展开的和没有手动展开的，手动展开的代码会更长，但是性能会有1-2%的提升（可以认为是误差了），求大佬解释，我们自己也没有一个很清晰的定论。
自我总结
如果没有改错的输入输出的话，我对我自己的修改是满意的。其他再厉害的方法，已经超出了我的能力范围，是还有待学习的领域。同时，这也是我们团队五个人齐心协力的结果，自问已经做到了最好，但是我们一致相信肯定还有更好的方法存在，也还希望大佬点评斧正。
或许可以优化的点？
- PetSc矩阵类型
- PetSc矩阵搜引
- 为什么OpenCAEPoro中要组装两次矩阵？（不确定是否正确的点）
Hi，此外我有幸作为打杂&后勤队员，能参加这个比赛倍感荣幸，我还在寻找RA，PhD或者其他工作机会，希望聚焦于控制，航空航天，机器人学，如果您愿意向我伸出橄榄枝，我感激不尽。