HPC-Game

HPC Game 0th. 作为HPC的小练习来完成，因为0th是一年前了，所以记录自己的解题应该没有问题。因为提供的超算平台现在不可以用了，相关文档现在也无法访问了，所以是直接做的题目。注册还可以看到题目。

主页：https://hpcgame.pku.edu.cn/

2024.1.16更新：

才只做了一点，现在主页进去已经是HPC Game 1th了，开始时间是下周，这次可以跟着参加尝试一下了，本篇也不用更新了。

一.高性能计算简介

1.实验室的新机器

本题的要求是编写一个脚本在集群上运行程序。

Slurm是一个开源的集群资源管理和作业调度系统，可以在Linux集群上运行各种并行程序，如MPI等。一些相关的信息可查看上交超算平台手册-Slurm作业调度系统或快速入门。

根据要求编写一个脚本启动程序，因为题目中说明不会使用sbatch，所以使用srun运行。

1
2
3

#!/bin/bash
srun -p compute -N1 -n1 -c1 ./program $1 > output.dat
seff $(cat job_id.dat) > seff.dat

如果用sbatch启动：

#!/bin/bash
#SBATCH -J program           
#SBATCH -p compute
#SBATCH -N 1                 
#SBATCH --cpus-per-task=1    
#SBATCH -t x:00       #任务最大运行时间
#SBATCH -o output.dat #将输出结果保存

./program $1
seff $(cat job_id.dat) > seff.dat

2.基础知识问答

TOP500中能效最高的超算：能效最高的应该是Henri，性能最高的是Frontier。
戈登贝尔奖是并行计算领域的最高奖项，设立于1987年。中国人第一次获得戈登贝尔奖是2016年，第一作者是杨超，现工作于计算所。
天河二号采用的是Xeon E5和Xeon Phi处理器，太湖之光采用申威众核处理器。TOP500曾排名最高的是太湖之光。
一般来说，服务器使用的 CPU 相比于同代同等级的游戏 CPU，核心更多，多核性能更强。单核频率更低，单核性能更弱。

5-7题：

假设单位内存是1，内存写为直写，全相联，块大小为1，替换策略为LRU，对于如下的矩阵乘法：

for m in [1, N]:
    for n in [1,N]:
        for k in [1, N]:
            C[m,n] = C[m,n] + A[m, k] * B[k, n]

5.对ABC均进行了N^3次读，对C进行了N^3次写。

6.如果有一个一级cache处理器，cache大小为3*N^2，运算开始时ABC在内存，结束后强制将cache写回内存(直写)。则会对内存进行3N^2次读运算，N^3+3N^2次写运算。

7.如果6中的cache大小为2N+1，则会对内存进行3N^3次读，N^3+2N+1次写。

3.简单题

本题的要求是使用P个核心，把N个 32 位有符号整数数组中所有元素 + 1，并且求原来（未加一）的数字之和(MOD 100001651)并输出。

这里复习一下mmap和文件操作，fstream的文件如果要mmap需要int fd = fs.rdbuf()->fd()转为整型，转完就可以把文件流关掉了，然后就可以mmap了。

#include <iostream>
#include <sys/mman.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <unistd.h>
#include <omp.h>
using namespace std;

int main() {
    int fd1 = open("input1.bin", O_RDONLY);
    int fd2 = open("output.bin", O_RDWR | O_CREAT, 0666);
    if (fd1 == -1 || fd2 == -1) {
        cout << "Open file error." << endl;
        return 0;
    }
    struct stat st;
    int r = fstat(fd1, &st);
    if (r == -1) {
        cout << "Get file size error." << endl;
        close(fd1);
        close(fd2);
        return 0;
    }
    int size = st.st_size;
    ftruncate(fd2, size-4);
    void *p1 = mmap(NULL, size, PROT_READ, MAP_PRIVATE, fd1, 0);
    if (p1 == MAP_FAILED ) {
        cout << "Map file error." << endl;
        close(fd1);
        close(fd2);
        return 0;
    }
    void *p2 = mmap(NULL, size-4, PROT_READ | PROT_WRITE, MAP_SHARED, fd2, 0);
    if (p2 == MAP_FAILED ) {
        cout << "Map file error." << endl;
        close(fd1);
        close(fd2);
        return 0;
    }

    int *input = (int *)p1;
    int *output = (int*)p2;
    int p = *input;
    int n = *(input + 1);
    cout << "p = " << p << ", n = " << n << endl;
    int sum = 0;
    #pragma omp prallel for reduction(+: sum) num_threads(p)
    for(int i=0;i<n;i++){
        int num = *(input+2+i);
        sum = (sum+num)%100001651;
        *(output+1+i) = num;
    }
    *output = sum%100001651;
    cout<<"sum = "<<sum<<endl;

    msync(p2, size-4, MS_SYNC);
    munmap(p1, size);
    munmap(p2, size-4);
    close(fd1);
    close(fd2);
    return 0;
}

二.并行与大规模

1.求积分

#include "rhs.h"
#include <omp.h>
#include <fstream>
#include <stdlib.h>
#include <iostream>
#include <iomanip>
using namespace std;

const int g_ncore = omp_get_num_procs(); 
int dtn(int n, int min_iteration){
    int max_threads = n/min_iteration;
    return max_threads>g_ncore?g_ncore:max(max_threads,1);
}

int main(int argc, char *argv[]){
    if(argc<2){
        cout << "usage: main <N>" <<endl;
        return 1;
    }
    int N = atoi(argv[1]);
    ofstream fs("output.dat");
    double sum = 0;
    #pragma omp parallel for reduction(+:sum) num_threads(dtn(N, 1e3))
    for(int i=1;i<=N;i++){
        double x = (i-0.5)/N;
        double res;
        rhs(x,res);
        sum += res;
    }
    fs << fixed << setprecision(12) << sum << "\n";
    fs.close();
    return 0;
}

2.乘一乘

写的时候遇到C++小问题，rhs提供的函数参数是unsigned int的引用，如果传int进去，会被转换为unsigned，变成临时值，会变成右值。

typedef unsigned int uint;
const int MAXN = 1e5+5;
const int g_ncore = omp_get_num_procs();
int dtn(int n, int min_iteration){
    int max_threads = n/min_iteration;
    return max_threads>g_ncore?g_ncore:max(max_threads,1);
}

void matrix_mul_parallel(double **result, int N1, int N2, int N3){
    clock_t start, end;
    start = clock();
    #pragma omp parallel for num_threads(dtn(N1,10))
    for(uint i=0;i<(uint)N1;i++){
        for(uint k=0;k<(uint)N2;k++){
            for(uint j=0;j<(uint)N3;j++){
                double va, vb;
                matA(i,k,va);
                matB(k,j,vb);
                result[i][j] += va*vb;  
            }
        }
    }
    end = clock();
    cout << "matrix_mul_parallel time = " << double(end - start) / CLOCKS_PER_SEC << "s" << endl;
}

这次还尝试了一下collapse和只对第一层循环并行的比较，时间结果如下：

./omp_cmul 100 1000 1000
matrix_mul time = 7.333s
matrix_mul_parallel time = 0.91s
matrix_mul_collapse time = 7.479s
PS C:\Users\Desktop\tmp\hpcgame\q2b> ./omp_cmul 1000 1000 1000
matrix_mul time = 73.78s
matrix_mul_parallel time = 8.02s
matrix_mul_collapse time = 9.282s

矩阵乘是有访问的空间局部性的，并且每一层循环内的负载是均衡的，因此将循环第一层并行化更快，collapse的效果和矩阵的形状有很大关系，但尝试的一些形状都比将第一层循环并行化更慢。

3.解方程

其实f()的值在每轮迭代都是一样的，可以先处理出来，不用每次算，但是反正也没办法交到平台测性能，这里就不优化了。遇到的一个小tip是reduction操作会算上原来的值，因此要注意进入并行循环之前的值。

void possion(double mat[MAXN][MAXN], int N){
    for(int i=0;i<=N;i++){
        for(int j=0;j<=N;j++) mat[i][j] = 0.0;
    }
    double eps;
    do{
        eps = 1e-16;
        #pragma omp parallel for reduction(max:eps) num_threads(dtn(N*N, 20))
        for(int i=1;i<N;i++){
            for(int j=1;j<N;j++){
                double xi, yi, va, tmp;
                xi = i/(double)N;
                yi = j/(double)N;
                rhs(xi, yi, va); 
                tmp = mat[i][j];
                mat[i][j] = (mat[i-1][j]+mat[i+1][j]+mat[i][j-1]+mat[i][j+1]+va/(double)N/(double)N)/4.0;
                eps = (eps<1e-15)?abs(mat[i][j]-tmp):max(abs(mat[i][j]-tmp),eps);
            }
        }
    }while(eps>1e-15);
}