随便写写
看到隔壁 HPC 组又开始办比赛了。比赛时大部分时间都在外面玩，作为懒狗的我本来没准备参加。但是后来看到奖金给的还挺多的，就随便打了打玩玩。最后在不超过 15h 工作量的情况下拿到了二等名次（差一点就一等了，不过要卷的话还是要花不少时间了），还不错。
不过我对 HPC 也不是很懂，只有很少的了解。在此处要特别感谢上个学期上《并行与分布式导论》课程的 lgj 老师，以及 ChatGPT 老师（同时感谢一下某位 Host 了一个 Bot 的群友，有一部分题目相关性比较弱的问题是薅的他的 GPT4）。
做题主要靠直接搜索，问 GPT，以及看上学期的课件。WriteUp 就随便写写了，因为我也不懂。
比赛的评测环境不是很稳。不过这种比赛确实运维太麻烦了，所以还是要感谢主办方同学。
直接提交就好了！不过这题写了的测试平台地址及账号密码我居然一直没看到。做了好几道题后感觉配不动环境了才仔细找了找发现。
照着 slurm 的文档，以及题目给的命令随便搞一下就好了。就是这题必须要在题目给的集群上做。
直接把题目描述复制给 GPT 老师。
CC = g++
MPICC = mpic++
NVCC = nvcc
CFLAGS = -fopenmp
MPIFLAGS =
CUDAFLAGS =
all: hello_omp hello_mpi hello_cuda
hello_omp: hello_omp.cpp
$(CC) $(CFLAGS) -o hello_omp hello_omp.cpp
hello_mpi: hello_mpi.cpp
$(MPICC) $(MPIFLAGS) -o hello_mpi hello_mpi.cpp
hello_cuda: hello_cuda.cu
$(NVCC) $(CUDAFLAGS) -o hello_cuda hello_cuda.cu
clean:
rm -f hello_omp hello_mpi hello_cuda
cmake_minimum_required(VERSION 3.10)
project(HelloWorld)
set(CMAKE_CXX_STANDARD 14)
find_package(CUDA REQUIRED)
find_package(MPI REQUIRED)
find_package(OpenMP REQUIRED)
if(OpenMP_CXX_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
endif()
cuda_add_executable(hello_cuda hello_cuda.cu)
add_executable(hello_mpi hello_mpi.cpp)
target_link_libraries(hello_mpi MPI::MPI_CXX)
add_executable(hello_omp hello_omp.cpp)
继续问 GPT 老师，以及上网查一些东西。基本都不难，虽然 GPT 老师会错几道。但反正每次提交会返回具体哪道题错了，多试几遍就好了。
t1 65.396
t2 105.26 107.14
t3 iv
t4 ii i i
t5 8 i ii
t6 iv
t7 iii,iv
t8 iv
t9 ii
t10 i
第五题这个 Intel 和 AMD 的区别还是挺经典的。经常能在某处辩 AVX 指令集作用的时候能看到这俩说法。
第七题我完全没搞懂题目选项到底指的是什么，反正第一个肯定错，第三个肯定对，剩下的只能根据推测多试几遍了。
最后做的题。最开始本来都没准备做，因为看大家环境经常炸。后来发现不做这题都拿不了二等了，就在最后一天做了一下。测试的时候开本地就能测正确性了。用示例脚本开的话也是有大于 10% 的成功率的，多试试就好。
让 GPT 老师写了份示例代码，然后改了改。这块我是手工去做 pipeline 了，虽然这个做的也不是很均匀。不过因为 Ray 自带了一定的 pipeline 特性（看那个 .get 的设计不难猜出），所以效率肯定是够了。
因为我比较懒，所以调度都是在主节点做的。实际上让各节点做自己的调度（以及数据的存取）应该会更快一点。
这题要喷一下出题人。似乎数据给的是 float 和 double 混着来的。使用 float 就能拿满分，然而什么都不管就会被直接提升到 double，似乎怎么搞都拿不了满分。显然我也没满。
import os
import ray
import numpy as np
import time
ray.init(address=f"{os.environ['RAY_CLUSTER_ADDR']}")
# ray.init(num_cpus=16)
print("Init Ready")
pg = ray.util.placement_group([{"CPU": 4},{"CPU": 4},{"CPU": 4},{"CPU": 4}], strategy="PACK")
@ray.remote(num_cpus=4)
class Layer:
def __init__(self):
pass
def loadfile(self,fn):
self.weight = np.load(fn)
def calc(self, inx):
return np.maximum(0,np.dot(inx,self.weight))
ray.get(pg.ready())
layers = [Layer.options(placement_group=pg).remote() for i in range(4)]
print("Ray Ready")
sttb=time.time()
for i in range(4):
layers[i].loadfile.remote(f"weights/weight_{i}.npy")
sttw=time.time()
print("Weight ready "+str(sttw-sttb))
os.makedirs("outputs",exist_ok=True)
pipelined=[None]*4
MAXN=100
for inputs in range(MAXN+5):
if pipelined[0] is not None:
idx,obj = pipelined[0]
result=ray.get(obj)
#save result
np.save(f"outputs/output_{idx}.npy",result)
print(f"saved {idx} {time.time()-sttb}")
for i in range(3):
pipelined[i]=pipelined[i+1]
if pipelined[i] is not None:
pipelined[i]=(pipelined[i][0],layers[3-i].calc.remote(pipelined[i][1]))
#load next input
if inputs<MAXN:
pipelined[3]=(inputs,layers[0].calc.remote(np.load(f"inputs/input_{inputs}.npy")))
print(f"load {inputs} {time.time()-sttb}")
ray.shutdown()
比较好玩的题。看到这题我就是奔着满分去了（最后是几乎满分了）。
显然，这题卡时间需要从两个地方卡：IO 时间，计算时间。
题目都说了这是可并行的校验。对 SHA256 的 Merkle–Damgård 结构熟悉的同学都知道，对于 SHA256，我们完全可以流式地输入，并提前做一些计算（好像这题不熟悉的也能看出来）。因此我们不难对计算部分进行并行：不同块的主要计算部分（即 1M 大小的文件内容部分）由不同的进程完成。当所有的主要部分都完成后再来从第一个块到最后一个块做链式的计算，即把上一个块的哈希值放到末尾，再计算当前块的哈希值。
为了提高并行效率，减少进程间通讯，我们让一段连续的块由一个进程完成。让每个进程都分到尽量等长的连续段。
接下来发现，时间中有很大一部分都被硬盘 IO 占据了。题目测评环境中着重强调了测试环境是没有 cache 的，因此我们要考虑对这块做优化。
我们的环境是多机多进程。一个显然的想法是让每一个进程做独立 IO，指读取自己需要的部分。这么操作后发现时间还是很慢，因为从硬盘复制一大坨的耗时和延迟都太大了。既然有这些问题，传统的做法是使用异步 IO，使得在计算的同时后台进行硬盘的读取。这个应该确实是没问题的，只可惜我不会，试了下发现好像有错，于是就摸了。
于是再尝试传统 IO 优化手段：mmap。测试发现在单进程环境下确实快了一点，但是多进程环境下慢成狗，不如直接大段读到内存。猜测是因为每次获取内容的时候才去陷入 kernel 然后调用硬盘 IO 操作，然后又有多个进程发送多个请求，导致 IO 操作碎片化太严重了。
最后，我猜想，既然题目环境说了硬盘没 cache，那我们能不能强行让他 cache 到内存。感觉是可以直接对操作系统发号施令的。于是问了问 GPT 老师有没有能搞 file IO cache hint 的，发现果然有个 fadvise 函数可以干这个。用了一下后发现飞快！具体流程也很简单：在 main 函数开始的时候直接用这个函数对内核发号施令，之后在每次要计算一个块的时候，再使用 pread 直接序列读取文件。
这个方法拿到了 99 分。其实不难发现这个操作还是可以优化的。比如说多个块 pread 合并一下，减少调用次数。以及 fadvise 不要一次全发完，可以逐渐 pipeline 地发。应该是能更好的。
#include <algorithm>
#include <chrono>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <iostream>
#include <mpi.h>
#include <openssl/evp.h>
#include <openssl/sha.h>
#include <sys/mman.h>
#include <sys/types.h>
#include <fcntl.h>
#include <unistd.h>
#include <sys/stat.h>
#include <assert.h>
#include <fcntl.h>
#include <stdio.h>
#include <unistd.h>
namespace fs = std::filesystem;
constexpr size_t BLOCK_SIZE = 1024 * 1024;
void checksum(uint8_t *data, size_t len, uint8_t *obuf);
void checksum_p(uint8_t *data, size_t len, uint8_t *obuf);
void print_checksum(std::ostream &os, uint8_t *md, size_t len);
int rank, nprocs;
int fd;
int main(int argc, char *argv[]) {
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
if (rank == 0) {
if (argc < 3) {
std::cout << "Usage: " << argv[0] << " <input_file> <output_file>"
<< std::endl;
MPI_Abort(MPI_COMM_WORLD, 1);
}
}
std::cout<<"rank "<<rank<<std::endl;
fs::path input_path = argv[1];
fs::path output_path = argv[2];
auto total_begin_time = std::chrono::high_resolution_clock::now();
auto file_size = fs::file_size(input_path);
std::cout << input_path << " size: " << file_size << std::endl;
uint8_t *buffer = nullptr;
if (file_size != 0) {
// buffer = new uint8_t[file_size];
// read the file content in binary format
// std::ifstream istrm(input_path, std::ios::binary);
// istrm.read(reinterpret_cast<char *>(buffer), file_size);
//Open file
fd = open(argv[1], O_RDONLY, 0);
assert(fd != -1);
//Execute mmap
// void* mmappedData = mmap(NULL, file_size, PROT_READ, MAP_PRIVATE | MAP_POPULATE, fd, 0);
// assert(mmappedData != MAP_FAILED);
// buffer=(uint8_t *)mmappedData;
//Write the mmapped data to stdout (= FD #1)
//Cleanup
}
size_t num_block = (file_size + BLOCK_SIZE - 1) / BLOCK_SIZE;
size_t mysb=num_block*rank/nprocs;
size_t myeb=num_block*(rank+1)/nprocs;
size_t rdoffset=mysb*BLOCK_SIZE;
size_t rdlen=(myeb-mysb)*BLOCK_SIZE;
if (posix_fadvise(fd, rdoffset, rdlen, POSIX_FADV_WILLNEED) == -1) {
perror("posix_fadvise");
close(fd);
return 1;
}
// buffer = new uint8_t[(myeb-mysb)*BLOCK_SIZE];
// pread(fd,buffer,rdlen,rdoffset);
// record begin time
auto begin_time = std::chrono::high_resolution_clock::now();
// calculate the checksum
uint8_t obuf[SHA512_DIGEST_LENGTH];
checksum_p(buffer, file_size, obuf);
// record end time
auto end_time = std::chrono::high_resolution_clock::now();
if (rank == nprocs-1) {
// print debug information
std::cout << "checksum: ";
print_checksum(std::cout, obuf, SHA512_DIGEST_LENGTH);
std::cout << std::endl;
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(
end_time - begin_time);
std::cout << "checksum time cost: " << std::dec << duration.count() << "ms"
<< std::endl;
// write checksum to output file
std::ofstream output_file(output_path);
print_checksum(output_file, obuf, SHA512_DIGEST_LENGTH);
auto total_end_time = std::chrono::high_resolution_clock::now();
auto total_duration = std::chrono::duration_cast<std::chrono::milliseconds>(
total_end_time - total_begin_time);
std::cout << "total time cost: " << total_duration.count() << "ms"
<< std::endl;
}
// delete[] buffer;
MPI_Finalize();
return 0;
}
typedef uint8_t shablock[SHA512_DIGEST_LENGTH];
void checksum_p(uint8_t *data, size_t len, uint8_t *obuf) {
size_t num_block = (len + BLOCK_SIZE - 1) / BLOCK_SIZE;
size_t mysb=num_block*rank/nprocs;
size_t myeb=num_block*(rank+1)/nprocs;
size_t mycnt=myeb-mysb;
shablock* prev_md= new shablock[num_block+1];
EVP_MD_CTX **ctxs=new EVP_MD_CTX *[num_block+1];
MPI_Status status;
// EVP_MD_CTX *ctx = EVP_MD_CTX_new();
EVP_MD *sha512 = EVP_MD_fetch(nullptr, "SHA512", nullptr);
for (int i = mysb; i < myeb; i++) {
ctxs[i]=EVP_MD_CTX_new();
uint8_t buffer[BLOCK_SIZE]{};
EVP_DigestInit_ex(ctxs[i], sha512, nullptr);
pread(fd,buffer,std::min(BLOCK_SIZE, len - i * BLOCK_SIZE),i* BLOCK_SIZE);
// std::memcpy(buffer, data + (i-mysb) * BLOCK_SIZE,
// std::min(BLOCK_SIZE, len - i * BLOCK_SIZE));
EVP_DigestUpdate(ctxs[i], buffer, BLOCK_SIZE);
}
if(mysb==0)
{
SHA512(nullptr, 0, prev_md[mysb]);
}
else
{
MPI_Recv((void*)prev_md[mysb],SHA512_DIGEST_LENGTH,MPI_BYTE,rank-1,0,MPI_COMM_WORLD,&status);
}
for (int i = mysb; i < myeb; i++) {
EVP_DigestUpdate(ctxs[i], prev_md[i], SHA512_DIGEST_LENGTH);
unsigned int len = 0;
EVP_DigestFinal_ex(ctxs[i], prev_md[i+1], &len);
}
if(myeb==num_block)
{
std::memcpy(obuf, prev_md[myeb], SHA512_DIGEST_LENGTH);
}
else
{
MPI_Send((void*)prev_md[myeb],SHA512_DIGEST_LENGTH,MPI_BYTE,rank+1,0,MPI_COMM_WORLD);
}
// EVP_MD_CTX_free(ctx);
// EVP_MD_free(sha512);
}
void print_checksum(std::ostream &os, uint8_t *md, size_t len) {
for (int i = 0; i < len; i++) {
os << std::setw(2) << std::setfill('0') << std::hex
<< static_cast<int>(md[i]);
}
}
没做。看着就不是能短时间内拿到分的题。
一看，这不是经典并行计算题目吗。直接上网搜 openmp avx512 matrix mul，抄到个 GitHub 代码。
这个代码直接拉下来发现有好几个问题：行列存储方式不一致；只针对了行列相等的情况。还有一些小问题。
对于行列存储方式不一致，一个清真的方法是直接对矩阵做个转置再传入，答案再来转置一下。不过我比较懒，就用了不清真的做法：直接将原矩阵当作转置后的结果，让他计算 ，这样传入是正的，算出来的也是正的。
对于只针对行列相等的情况，当然是一个个分析每个变量的含义（痛苦面具），把 n 替换成 n1 之类的（感觉不如重写一遍）。
交上去后发现只有一点分。本地测了一下，发现在某个规模后（印象里是 2048）时间增速超过了复杂度增速。上过 ICS 的同学不难猜出，这是因为 cache 失效了。而应对 cache 失效的方法也很简单，就是让每列不要按 2^k 对齐就好了。因此我在每个矩阵开的数组中，加了一点虚拟的行和列，使得实际同列的位置不对齐。注意 AVX512 最好让内存地址以 64 对齐。具体分块大小以及添加的行列大小就多试试，随便调一下参。
最后拿了 70 多分。vtune 看了一下，瓶颈主要在内存，感觉要提高的话得去精细地优化内存布局了。不会。
#include <iostream>
#include <chrono>
#include <x86intrin.h>
using namespace std;
/*
Optimized version of DGEMM using C intrinsics to generate the AVX subword-parallel instructions
for the x86, loop unrolling and blocking to create more opportunities for instruction-level parallelism.
*/
static constexpr uint32_t BLOCKSIZE = 64;
static inline void do_block(const uint32_t n1, const uint32_t n2, const uint32_t n3, const uint32_t si, const uint32_t sj, const uint32_t sk, const double *A, const double*B, double *C)
{
constexpr uint32_t UNROLL = 8;
//set unroll to 8
/*
b(A) n2*n3
a(B) n1*n2
c(C) n1*n3
*/
for( uint32_t i = si; i < si + BLOCKSIZE; i += UNROLL * 8)
{
for( uint32_t j = sj; j < sj + BLOCKSIZE; ++j)
{
__m512d c[UNROLL];
#pragma GCC unroll 8
for( uint32_t r = 0; r < UNROLL; r++)
{
c[r] = _mm512_loadu_pd(C + i + r * 8 + j * n3); //[ UNROLL];
}
for( uint32_t k = sk; k < sk + BLOCKSIZE; k++ )
{
__m512d bb = _mm512_broadcastsd_pd(_mm_load_sd(B + j * n2 + k));
#pragma GCC unroll 8
for( uint32_t r = 0; r < UNROLL; r++)
{
c[r] = _mm512_fmadd_pd(_mm512_loadu_pd(A + n3 * k + r * 8 + i), bb, c[r]);
}
}
#pragma GCC unroll 8
for( uint32_t r = 0; r < UNROLL; r++)
{
_mm512_storeu_pd(C + i + r * 8 + j * n3, c[r]);
}
}
}
}
void dgemm_openmp(const double* A, const double* B, double* C,const uint32_t n1, const uint32_t n2, const uint32_t n3)
{
#pragma omp parallel for
for( uint32_t sj = 0; sj < n1; sj += BLOCKSIZE )
{
for( uint32_t si = 0; si < n3; si += BLOCKSIZE )
{
for( uint32_t sk = 0; sk < n2; sk += BLOCKSIZE )
{
// c[i][j] += A[i][k] * B[k][j]
do_block(n1,n2,n3, si, sj, sk, A, B, C);
}
}
}
}
void mul(double* a, double* b, double* c, uint64_t n1, uint64_t n2, uint64_t n3) {
for (int i = 0; i < n1; i++) {
for (int j = 0; j < n2; j++) {
for (int k = 0; k < n3; k++) {
c[i * n3 + k] += a[i * n2 + j] * b[j * n3 + k];
}
}
}
}
int main() {
uint64_t n1, n2, n3;
uint64_t nn1, nn2, nn3;
FILE* fi;
fi = fopen("conf.data", "rb");
fread(&n1, 1, 8, fi);
fread(&n2, 1, 8, fi);
fread(&n3, 1, 8, fi);
nn1=n1+BLOCKSIZE;
nn2=n2+BLOCKSIZE;
nn3=n3+BLOCKSIZE;
// double* a = (double*)malloc(n1 * n2 * 8+4*64);
// double* b = (double*)malloc(n2 * n3 * 8+4*64);
// double* c = (double*)malloc(n1 * n3 * 8+4*64);
double* a = (double*)malloc(nn1 * nn2 * 8+4*64);
double* b = (double*)malloc(nn2 * nn3 * 8+4*64);
double* c = (double*)malloc(nn1 * nn3 * 8+4*64);
a = (double*)(((uint64_t)a + 63) & ~63);
b = (double*)(((uint64_t)b + 63) & ~63);
c = (double*)(((uint64_t)c + 63) & ~63);
// fread(a, 1, n1 * n2 * 8, fi);
// fread(b, 1, n2 * n3 * 8, fi);
for(int i=0;i<n1;i++)
fread(a+i*nn2, 1, n2 * 8, fi);
for(int i=0;i<n2;i++)
fread(b+i*nn3, 1, n3 * 8, fi);
fclose(fi);
for (uint64_t i = 0; i < nn1; i++) {
for (uint64_t k = 0; k < nn3; k++) {
c[i * nn3 + k] = 0;
}
}
auto t1 = std::chrono::steady_clock::now();
// mul(a, b, c, n1, n2, n3);
dgemm_openmp(b,a, c, nn1, nn2, nn3);
// dgemm_avx512(b,a, c, n1, n2, n3);
//B.T*A.T=C.T lol
auto t2 = std::chrono::steady_clock::now();
int d1 = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
printf("%d %d %d\n", nn1,nn2,nn3);
printf("%d\n", d1);
fi = fopen("out.data", "wb");
// fwrite(c, 1, n1 * n3 * 8, fi);
for(int i=0;i<n1;i++)
fwrite(c+i*nn3, 1, n3 * 8, fi);
fclose(fi);
return 0;
}
入门 AVX512 优化题。
首先我们不难随手改一个 OpenMP 的多数据并行版本。也不难随手改一个 AVX512 版本，每步对多个数据同时做迭代。但是最后发现只能拿到七十多分。更恐怖的是，程序的时间是满分时间好几倍。事实上，我们去看这个版本的代码，不难发现，每一次循环的 AVX 指令都是依赖上一次循环的计算结果的。而 AVX 指令的吞吐量很大，延迟也很大。这造成了严重的数据依赖，下一次循环要使用 AVX 计算时，上一次循环的结果还没出来。发现这点后，我们不难再做优化：每一次循环我们同时执行八条 AVX 指令（因为评测用的 CPU 的 AVX 指令延迟为 4 cycle，吞吐为 0.5 cycle）。这样 CPU 连续发射的八条 AVX 指令之间就是没有数据依赖的了。之后当首先发射的指令计算结果出来后，我们正好开始发射下一次的八条指令。这样我们就把 CPU 的吞吐打满了。
#include <immintrin.h>
#include <iostream>
#include <chrono>
#include <omp.h>
#include <string.h>
// double it(double r, double x, int64_t itn) {
// for (int64_t i = 0; i < itn; i++) {
// x = r * x * (1.0 - x);
// }
// return x;
// }
// void itv(double r, double* x, int64_t n, int64_t itn) {
// #pragma omp parallel for
// for (int64_t i = 0; i < n; i++) {
// x[i] = it(r, x[i], itn);
// }
// }
void itv(double r, double* x, int64_t n, int itn) {
__m512d vr = _mm512_set1_pd(r);
__m512d vone = _mm512_set1_pd(1.0);
#pragma omp parallel for
for (int64_t i = 0; i < n; i += 8*8) {
__m512d vx1 = _mm512_loadu_pd(x + i);
__m512d vx2 = _mm512_loadu_pd(x + i + 8);
__m512d vx3 = _mm512_loadu_pd(x + i + 16);
__m512d vx4 = _mm512_loadu_pd(x + i + 24);
__m512d vx5 = _mm512_loadu_pd(x + i + 32);
__m512d vx6 = _mm512_loadu_pd(x + i + 40);
__m512d vx7 = _mm512_loadu_pd(x + i + 48);
__m512d vx8 = _mm512_loadu_pd(x + i + 56);
for (int j = 0; j < itn; j+=64) {
#pragma GCC unroll 64
for (int k = 0; k < 64; k++) {
// vx = _mm512_mul_pd(_mm512_mul_pd(vr, vx), _mm512_sub_pd(vone, vx));
vx1 = _mm512_mul_pd(_mm512_mul_pd(vr, vx1), _mm512_sub_pd(vone, vx1));
vx2 = _mm512_mul_pd(_mm512_mul_pd(vr, vx2), _mm512_sub_pd(vone, vx2));
vx3 = _mm512_mul_pd(_mm512_mul_pd(vr, vx3), _mm512_sub_pd(vone, vx3));
vx4 = _mm512_mul_pd(_mm512_mul_pd(vr, vx4), _mm512_sub_pd(vone, vx4));
vx5 = _mm512_mul_pd(_mm512_mul_pd(vr, vx5), _mm512_sub_pd(vone, vx5));
vx6 = _mm512_mul_pd(_mm512_mul_pd(vr, vx6), _mm512_sub_pd(vone, vx6));
vx7 = _mm512_mul_pd(_mm512_mul_pd(vr, vx7), _mm512_sub_pd(vone, vx7));
vx8 = _mm512_mul_pd(_mm512_mul_pd(vr, vx8), _mm512_sub_pd(vone, vx8));
}
}
_mm512_storeu_pd(x + i, vx1);
_mm512_storeu_pd(x + i + 8, vx2);
_mm512_storeu_pd(x + i + 16, vx3);
_mm512_storeu_pd(x + i + 24, vx4);
_mm512_storeu_pd(x + i + 32, vx5);
_mm512_storeu_pd(x + i + 40, vx6);
_mm512_storeu_pd(x + i + 48, vx7);
_mm512_storeu_pd(x + i + 56, vx8);
}
}
int main(){
FILE* fi;
fi = fopen("conf.data", "rb");
int64_t itn;
double r;
int64_t n;
double* x;
fread(&itn, 1, 8, fi);
fread(&r, 1, 8, fi);
fread(&n, 1, 8, fi);
x = (double*)malloc((n+256) * 8);
// align to 64 bytes
x = (double*)(((uint64_t)x + 63) & ~63);
fread(x, 1, n * 8, fi);
fclose(fi);
//set 0 to the rest
memset(x + n, 0, 64 * 8);
auto t1 = std::chrono::steady_clock::now();
itv(r, x, n, itn);
auto t2 = std::chrono::steady_clock::now();
int d1 = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
printf("%d\n", d1);
fi = fopen("out.data", "wb");
fwrite(x, 1, n * 8, fi);
fclose(fi);
return 0;
}
比较好玩的题。我是在题目给的示例代码上做的优化。能优化的点还挺多的。
不难发现代码中最耗时的部分在第三步。我们先做一些显然的优化。比如说多加一些 openmp parallel for。注意我们需要先改写原来的 vector 迭代器枚举的方式。另外有些地方需要 atomic。经过进一步的分析也可以发现，最耗时的是第三步的 mmv 函数（也是需要 atomic 的函数）。这一步做完之后就能拿到接近 20 分了。
之后我们再做一些可能没那么显著的优化，这要对代码进行大量的重构，但这也使得后续的高级优化更方便了。事实上测试发现确实不显著，做完这些后分数就高了不到 5 分。
一个是不要写 cpp 了。先把绝大多数的 vector 都改成 malloc 出的数组，以便于后续实现内存上的优化（比如对齐什么的）。还可以小上一个 AVX512（不过这块我没仔细搞，我怀疑没有把吞吐打满，不过看 vtune 的话瓶颈应该在内存上，优化这个不知道有没有用），这个让 GPT 老师去做还挺快的。还做了一个我以为重量级级的优化：把 mmv 的 atomic 去掉。具体思路很简单，就是把那个 vector 内的内容按 row 排序，之后每个 row 只由一个线程负责就行。由于每个 row 的任务量不同，最开始我是让每个线程平均分配到相同的 row 数量，但这可能导致实际任务量差距比较大。可能需要调一下调度器才能使得效率最优，我最后是强行手分了任务。关于排序，我也是使用的基数排序，这样效率会比较快（多线程的效果也会好一点）。后来发现，这个优化只在线程数比较少的比较明显，线程多的情况下不是很明显，不清楚为什么，所以就只快了一点。
做完这些优化后，只快了一点。但是这时候我们发现，尽管耗时最大的过程在第三步，但是完全串行第二步的 act 也占了接近一半的时间。我们看看这个能不能优化掉。注意到对于原来的代码，我们完全没法让 for 循环并行。在 if 内修改状态的代码加入 critical 之类的只会更慢（使用第三方库 atomic 的 vector 可能会好一些）。试图手动维护一个全局的 atomic cnt 也会导致每次写入事件都需要锁一下变量，很慢。而我们后面又要做基数排序，因此我们不妨考虑能不能不提前写入数组，直接在计算的过程中做基数排序。基数排序分为两个步骤，第一步获取每个关键字的元素数量，并在最后求前缀和获取每个关键字对应的下标范围；第二步将元素依次填入到对应的下标范围中。我们可以做两遍参考代码的大循环，分别完成第一步和第二步。只要这两个循环能并行，额外多算一遍的开销也是可以承受的。我们从上面的经验就可以发现，对于冲突很少的 atomic 操作，效率是很高的。而在基数排序的第一步和第二步中，我们每次只会访问和当前元素关键字相关的位置，这意味着不同线程的操作是很难发生冲突的（因为关键字很多）。当然我们还有个求前缀和的操作，不过因为关键字的个数只有不到元素个数的 1/10，因此串行求就好了（其实也可以大力并行的，比如分两遍什么的，不过这只能优化几十毫秒最多几百毫秒的时间吧就没做）。
做完这个优化之后，就有接近 60 分了。之后又做了一些小优化，拿了六十多分。比如对于 getsp 函数中调用的某些可以合并的操作，我们就合并一下，减少内存访问次数。
#include <vector>
#include <unordered_map>
#include <stdint.h>
#include <math.h>
#include <omp.h>
#include <chrono>
#include <immintrin.h>
#ifdef _WIN32
#include <windows.h>
#define popcnt __popcnt64
#else
#define popcnt __builtin_popcountll
#endif
#include <iostream>
#include<cstring>
using namespace std;
typedef std::unordered_map<uint64_t, size_t> map_t;
struct restrict_t {
int offset, range, minocc, maxocc;
int occ;
uint64_t substate;
};
template <typename VT>
struct term_t {
VT value;
uint64_t an, cr, signmask, sign;
};
struct sparse_t {
std::vector<size_t> row;
std::vector<size_t> col;
std::vector<double> data;
};
struct newsparse
{
//sorts
std::vector<size_t> rowix;
std::vector<size_t> rowcnt;
std::vector<size_t> rowchecks;
// std::vector<size_t> nextp;
// std::vector<size_t> colv;
// std::vector<double> datav;
};
int itrest(restrict_t& rest) {
if (!rest.substate) {
goto next;
}
{
uint64_t x = rest.substate & (-(int64_t)rest.substate);
uint64_t y = x + rest.substate;
rest.substate = y + (y ^ rest.substate) / x / 4;
}
if (rest.substate >> rest.range) {
next:
if (rest.occ == rest.maxocc) {
rest.occ = rest.minocc;
rest.substate = (uint64_t(1) << rest.occ) - 1;
return 1;
}
rest.occ++;
rest.substate = (uint64_t(1) << rest.occ) - 1;
return 0;
}
return 0;
}
int itrest(std::vector<restrict_t>& rest) {
for (restrict_t& re : rest) {
if (!itrest(re)) {
return 0;
}
}
return 1;
}
uint64_t getstate(const std::vector<restrict_t>& rest) {
uint64_t state = 0;
for (const restrict_t& re : rest) {
state |= re.substate << re.offset;
}
return state;
}
int generatetable(std::vector<uint64_t>& table, map_t& map, std::vector<restrict_t>& rest) {
for (restrict_t& re : rest) {
re.occ = re.minocc;
re.substate = (uint64_t(1) << re.occ) - 1;
}
size_t index = 0;
do {
uint64_t state = getstate(rest);
table.push_back(state);
map.insert(std::make_pair(state, index));
index++;
} while (!itrest(rest));
return 0;
}
template <typename VT>
term_t<VT> getterm(VT value, const std::vector<int>& cr, const std::vector<int>& an) {
term_t<VT> term;
term.value = value;
term.an = 0;
term.cr = 0;
term.signmask = 0;
uint64_t signinit = 0;
for (int x : an) {
uint64_t mark = uint64_t(1) << x;
term.signmask ^= (mark - 1) & (~term.an);
term.an |= mark;
}
for (int x : cr) {
uint64_t mark = uint64_t(1) << x;
signinit ^= (mark - 1) & term.cr;
term.signmask ^= (mark - 1) & (~term.an) & (~term.cr);
term.cr |= mark;
}
term.sign = popcnt(signinit ^ (term.signmask & term.an));
term.signmask = term.signmask & (~term.an) & (~term.cr);
return term;
}
template <typename VT>
int act(std::vector<size_t>& row, std::vector<size_t>& col, std::vector<VT>& data, const std::vector<term_t<VT>>& op, const std::vector<uint64_t>& table, const map_t& map,newsparse& nopm) {
int64_t n = table.size();
nopm.rowix.resize(n);
nopm.rowcnt.resize(n);
int mxt=omp_get_max_threads();
nopm.rowchecks.resize(mxt+1);
#pragma omp parallel for schedule(guided)
// #pragma omp parallel for
for (int64_t i = 0; i < n; i++) {
uint64_t srcstate = table[i];
for (const term_t<VT>& term : op) {
if ((srcstate & term.an) == term.an) {
uint64_t dststate = srcstate ^ term.an;
if ((dststate & term.cr) == 0) {
dststate ^= term.cr;
auto it = map.find(dststate);
if (it != map.end()) {
uint64_t sign = term.sign + popcnt(srcstate & term.signmask);
VT v = term.value;
if (sign & 1) {
v = -v;
}
// #pragma omp critical
{
#pragma omp atomic
nopm.rowcnt[it->second]++;
}
}
}
}
}
}
nopm.rowix[0]=nopm.rowcnt[0];
for (int64_t i = 1; i < n; i++)
{
nopm.rowix[i]=nopm.rowix[i-1]+nopm.rowcnt[i];
}
const int multiper=5;
size_t avgc=(nopm.rowix[n-1]+n*multiper)/mxt;
#pragma omp parallel for
for (int64_t i = 1; i < n-1; i++)
{
if((nopm.rowix[i]+i*multiper)/avgc!=(nopm.rowix[i+1]+(i+1)*multiper)/avgc)
{
size_t lxp=(nopm.rowix[i+1]+(i+1)*multiper)/avgc;
nopm.rowchecks[lxp]=i+1;
}
}
nopm.rowchecks[mxt]=n;
//sanity check
for(int i=1;i<=mxt;i++)
{
if(nopm.rowchecks[i]==0 || nopm.rowchecks[i-1]>=nopm.rowchecks[i])
{
printf("%d %d %d %d\n",i,mxt,nopm.rowchecks[i-1],nopm.rowchecks[i]);
return -1;
}
}
size_t nn=nopm.rowix[n-1];
row.resize(nn);
col.resize(nn);
data.resize(nn);
// #pragma omp parallel for schedule(dynamic)
#pragma omp parallel for schedule(guided)
// #pragma omp parallel for
for (int64_t i = 0; i < n; i++) {
uint64_t srcstate = table[i];
for (const term_t<VT>& term : op) {
if ((srcstate & term.an) == term.an) {
uint64_t dststate = srcstate ^ term.an;
if ((dststate & term.cr) == 0) {
dststate ^= term.cr;
auto it = map.find(dststate);
if (it != map.end()) {
uint64_t sign = term.sign + popcnt(srcstate & term.signmask);
VT v = term.value;
if (sign & 1) {
v = -v;
}
// #pragma omp critical
{
int64_t idx;
#pragma omp atomic capture
idx=--nopm.rowix[it->second];
// row[idx]=it->second;
col[idx]=i;
data[idx]=v;
}
}
}
}
}
}
printf("OK\n");
return 0;
}
int readss(FILE* fi, std::vector<uint64_t>& table, map_t& map) {
int n;
fread(&n, 1, 4, fi);
std::vector<restrict_t> restv(n);
for (auto& rest : restv) {
fread(&rest, 1, 16, fi);
}
generatetable(table, map, restv);
return 0;
}
int readop(FILE* fi, std::vector<term_t<double>>& op) {
int n, order;
fread(&n, 1, 4, fi);
fread(&order, 1, 4, fi);
std::vector<double> v(n);
fread(v.data(), 1, 8 * n, fi);
std::vector<int> rawterm(order);
std::vector<int> cr, an;
for (int i = 0; i < n; i++) {
fread(rawterm.data(), 1, 4 * order, fi);
int tn = rawterm[0];
for (int j = 0; j < tn; j++) {
int type = rawterm[tn * 2 - 1 - j * 2];
if (type) {
cr.push_back(rawterm[tn * 2 - j * 2]);
}
else {
an.push_back(rawterm[tn * 2 - j * 2]);
}
}
op.push_back(getterm(v[i], cr, an));
cr.clear();
an.clear();
}
return 0;
}
double* genvec(size_t n)
{
double* v=(double*)calloc(n*8+64*(2+8)+8,1);
v = (double*)(((uint64_t)v + 63+8) & ~63);
*((size_t*)((uint64_t)v-8))=n;
return v;
}
size_t veclen(double* v)
{
return *((size_t*)((uint64_t)v-8));
}
#define BLOCKSIZE 64
const int UNROLL=8;
//out=m*v;
double mmv(double* out, const sparse_t& m, double* v,const newsparse& nopm) {
size_t n=veclen(out);
size_t nn=m.data.size();
//faster a bit
double s=0;
#pragma omp parallel reduction(+: s)
{
int idx=omp_get_thread_num();
size_t tod=nopm.rowchecks[idx+1];
for (size_t r = nopm.rowchecks[idx]; r < tod; r++) {
double st=0;
size_t bs=nopm.rowix[r];
size_t rc=nopm.rowcnt[r];
#pragma GCC ivdep
for(size_t p=0;p<rc;p++)
{
st+=m.data[bs+p]*v[m.col[bs+p]];
// out[m.row[bs+p]]+=m.data[bs+p]*v[m.col[bs+p]];
}
out[r]=st;
s+=st*v[r];
}
}
return s;
}
//v*=s;
void msv(const double s, double* v) {
size_t n=veclen(v);
__m512d scalar = _mm512_set1_pd(s);
#pragma omp parallel for
for (size_t i = 0; i < n; i += BLOCKSIZE) {
#pragma GCC unroll UNROLL
for(int j=0;j<UNROLL;j++)
{
__m512d vec = _mm512_load_pd(&v[i+8*j]);
__m512d result = _mm512_mul_pd(vec, scalar);
_mm512_store_pd(&v[i+8*j], result);
}
}
}
//v'*v;
double norm2(double* v) {
double s = 0;
size_t n=veclen(v);
#pragma omp parallel for reduction(+: s)
for (size_t i = 0; i < n; i += BLOCKSIZE) {
#pragma GCC unroll UNROLL
for(int j=0;j<UNROLL;j++)
{
__m512d vec = _mm512_load_pd(&v[i+8*j]);
__m512d result = _mm512_mul_pd(vec, vec);
s += _mm512_reduce_add_pd(result);
}
}
return s;
}
double avv1(double* v1, const double s, double* v2) {
size_t n=veclen(v1);
double rt = 0;
__m512d scalar = _mm512_set1_pd(s);
#pragma omp parallel for reduction(+: rt)
for (size_t i = 0; i < n; i += BLOCKSIZE) {
#pragma GCC unroll UNROLL
for(int j=0;j<UNROLL;j++)
{
__m512d vec1 = _mm512_load_pd(&v1[i+8*j]);
__m512d vec2 = _mm512_load_pd(&v2[i+8*j]);
__m512d result = _mm512_fmadd_pd(scalar, vec2, vec1);
_mm512_store_pd(&v1[i+8*j], result);
rt+=_mm512_reduce_add_pd(_mm512_mul_pd(result,result));
}
}
return rt;
}
double avv2(double* v1, const double sa, double* v2, const double sb, double* v3) {
size_t n=veclen(v1);
double rt = 0;
__m512d scalara = _mm512_set1_pd(sa);
__m512d scalarb = _mm512_set1_pd(sb);
#pragma omp parallel for reduction(+: rt)
for (size_t i = 0; i < n; i += BLOCKSIZE) {
#pragma GCC unroll UNROLL
for(int j=0;j<UNROLL;j++)
{
__m512d vec1 = _mm512_load_pd(&v1[i+8*j]);
__m512d vec2 = _mm512_load_pd(&v2[i+8*j]);
__m512d vec3 = _mm512_load_pd(&v3[i+8*j]);
__m512d result = _mm512_fmadd_pd(scalara, vec2, vec1);
result = _mm512_fmadd_pd(scalarb, vec3, result);
_mm512_store_pd(&v1[i+8*j], result);
rt+=_mm512_reduce_add_pd(_mm512_mul_pd(result,result));
}
}
return rt;
}
void getsp(double* out, int itn, const sparse_t m, double* v,const newsparse& nopm) {
out[0] = sqrt(norm2(v));
msv(1.0 / out[0], v);
double *a=genvec(itn),*b=genvec(itn-1);
double *v_=genvec(veclen(v)),*v__=genvec(veclen(v));
for (int i = 0; i < itn; i++) {
swap(v__,v_);
swap(v_,v);
double tt=mmv(v, m, v_,nopm);
a[i]=tt;
// a[i] = dot(v, v_);
double rt=0;
if (i < itn - 1) {
if (i == 0) {
// avv(v, -a[i], v_);
rt=avv1(v,-a[i],v_);
}
else {
// avv(v, -a[i], v_);
// avv(v, -b[i - 1], v__);
rt=avv2(v,-a[i],v_,-b[i-1],v__);
}
// b[i] = sqrt(norm2(v));
b[i] = sqrt(rt);
msv(1.0 / b[i], v);
}
}
#pragma omp parallel for
for (int i = 0; i < itn; i++) {
out[1 + i] = a[i];
}
#pragma omp parallel for
for (int i = 0; i < itn - 1; i++) {
out[1 + itn + i] = b[i];
}
}
int main()
{
FILE* fi;
std::vector<uint64_t> table;
map_t map;
std::vector<term_t<double>> op;
fi = fopen("conf.data", "rb");
auto t1 = std::chrono::steady_clock::now();
readss(fi, table, map);
auto t2 = std::chrono::steady_clock::now();
readop(fi, op);
sparse_t opm;
newsparse nopm;
act(opm.row, opm.col, opm.data, op, table, map,nopm);
auto t3 = std::chrono::steady_clock::now();
int d1 = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();
int d2 = std::chrono::duration_cast<std::chrono::milliseconds>(t3 - t2).count();
printf("Tt: %d,%d\n", d1, d2);
int itn;
fread(&itn, 1, 4, fi);
double* vv=genvec(table.size());
fread(vv, 1, table.size() * 8, fi);
fclose(fi);
printf("%d %d %lld %d\n", table.size(), itn,opm.data.size(),op.size());
double* result=genvec(itn*2);
getsp(result, itn, opm, vv,nopm);
auto t4 = std::chrono::steady_clock::now();
fi = fopen("out.data", "wb");
fwrite(result, 1, 16 * itn, fi);
fclose(fi);
int d3 = std::chrono::duration_cast<std::chrono::milliseconds>(t4 - t3).count();
printf("%d,%d,%d\n", d1, d2, d3);
std::cout << "Hello World!\n";
}
不会。一看就是要写清真光栅化。图形学半路退课的人表示完全不懂（虽然好像也不讲怎么手写光栅化）。不知道是哪个学了图形学的出题人出的。
打开一看，题面怎么这么复杂。看了一下代码，怎么这么简单。但是又思考了一下，这个代码加并行不难，就是 atomic 冲突问题可能很难解决。不知道怎么这么多人满分的。
后来随手算了一下加速比：怎么 64 核要求这么点加速比。于是就随便加了个 parallel for 和 atomic，就过了。
#include <array>
#include <fstream>
#include <iostream>
#include <omp.h>
#include <vector>
#include <cmath>
#include <tuple>
using std::vector, std::array, std::tuple, std::string;
void particle2grid(int resolution, int numparticle,
const vector<double> &particle_position,
const vector<double> &particle_velocity,
vector<double> &velocityu, vector<double> &velocityv,
vector<double> &weightu, vector<double> &weightv) {
double grid_spacing = 1.0 / resolution;
double inv_grid_spacing = 1.0 / grid_spacing;
auto get_frac = [&inv_grid_spacing](double x, double y) {
int xidx = floor(x * inv_grid_spacing);
int yidx = floor(y * inv_grid_spacing);
double fracx = x * inv_grid_spacing - xidx;
double fracy = y * inv_grid_spacing - yidx;
return tuple(array<int, 2>{xidx, yidx},
array<double, 4>{fracx * fracy, (1 - fracx) * fracy,
fracx * (1 - fracy),
(1 - fracx) * (1 - fracy)});
};
#pragma omp parallel for
for (int i = 0; i < numparticle; i++) {
array<int, 4> offsetx = {0, 1, 0, 1};
array<int, 4> offsety = {0, 0, 1, 1};
auto [idxu, fracu] =
get_frac(particle_position[i * 2 + 0],
particle_position[i * 2 + 1] - 0.5 * grid_spacing);
auto [idxv, fracv] =
get_frac(particle_position[i * 2 + 0] - 0.5 * grid_spacing,
particle_position[i * 2 + 1]);
#pragma GCC unroll 4
for (int j = 0; j < 4; j++) {
int tmpidx = 0;
tmpidx =
(idxu[0] + offsetx[j]) * resolution + (idxu[1] + offsety[j]);
#pragma omp atomic
velocityu[tmpidx] += particle_velocity[i * 2 + 0] * fracu[j];
#pragma omp atomic
weightu[tmpidx] += fracu[j];
tmpidx = (idxv[0] + offsetx[j]) * (resolution + 1) +
(idxv[1] + offsety[j]);
#pragma omp atomic
velocityv[tmpidx] += particle_velocity[i * 2 + 1] * fracv[j];
#pragma omp atomic
weightv[tmpidx] += fracv[j];
}
}
}
int main(int argc, char *argv[]) {
if (argc < 2) {
printf("Usage: %s inputfile\n", argv[0]);
return -1;
}
string inputfile(argv[1]);
std::ifstream fin(inputfile, std::ios::binary);
if (!fin) {
printf("Error opening file");
return -1;
}
int resolution;
int numparticle;
vector<double> particle_position;
vector<double> particle_velocity;
fin.read((char *)(&resolution), sizeof(int));
fin.read((char *)(&numparticle), sizeof(int));
particle_position.resize(numparticle * 2);
particle_velocity.resize(numparticle * 2);
printf("resolution: %d\n", resolution);
printf("numparticle: %d\n", numparticle);
fin.read((char *)(particle_position.data()),
sizeof(double) * particle_position.size());
fin.read((char *)(particle_velocity.data()),
sizeof(double) * particle_velocity.size());
vector<double> velocityu((resolution + 1) * resolution, 0.0);
vector<double> velocityv((resolution + 1) * resolution, 0.0);
vector<double> weightu((resolution + 1) * resolution, 0.0);
vector<double> weightv((resolution + 1) * resolution, 0.0);
string outputfile;
particle2grid(resolution, numparticle, particle_position,
particle_velocity, velocityu, velocityv, weightu,
weightv);
outputfile = "output.dat";
std::ofstream fout(outputfile, std::ios::binary);
if (!fout) {
printf("Error output file");
return -1;
}
fout.write((char *)(&resolution), sizeof(int));
fout.write(reinterpret_cast<char *>(velocityu.data()),
sizeof(double) * velocityu.size());
fout.write(reinterpret_cast<char *>(velocityv.data()),
sizeof(double) * velocityv.size());
fout.write(reinterpret_cast<char *>(weightu.data()),
sizeof(double) * weightu.size());
fout.write(reinterpret_cast<char *>(weightv.data()),
sizeof(double) * weightv.size());
return 0;
}
看了一眼。怎么有比赛出题让选手修库 bug 的。这题还没发实时评测，直接就鸽了。
看都没看，直接就鸽了。
日期： 2024-02-03