forrestguang的博客

之前关于flash attention的实现参考lash attention的数学变换：给定三个矩阵Q,K,V，形状都是[N,d]，计算S=QK.T，然后针对dim=1做softmax，然后和V继续做矩阵乘法得到形状为[N,d]的输出矩阵O，即O=softmax(QK.T,dim=1)V。下面本人的实现以下图为基础，从下图也可以看出attention算法的关键在于对Q,K,V的分块，因此我们需要首先考虑清楚编程过程中一个线程块到底处理哪些数据。

回顾pytorch里面搭建的全连接神经网络，我们以最简单的全连接神经网络举例子，这种神经网络涉及向量的传输，矩阵乘法计算，以及激活函数定义，其中的矩阵计算以及激活函数就可以理解为神经网络里面的算子。之前在上一篇文章添加链接描述我们介绍了CUDA版本的矩阵乘法，这里我们不妨以OpenMP为工具搭建一个CPU版本的矩阵乘法算子。

关于CMakeLists.txt的相关介绍，这里不赘诉，本人的出发点是借助于CMakeLists.txt掌握基本的C++构建项目流程，下面是本人根据网络资料以及个人实践掌握的资料。

我们以矩阵乘法C=AB，其中A形状为[M,K],B形状为[K,N]，C形状为[M,N]举例子，下面的分布式算法我们默认以MPI来切分数据，其中每个进程之前的数据都是私有，进程之间的数据交互使用通信来完成。A=(a0,0a0,1a0,2⋯a0,K−1a1,0a1,1a1,2⋯a1,K−1⋮⋮⋱⋮aM−2,0⋯aM−2,K−3aM−2,K−2aM−2,K−1aM−1,0aM−1,1⋯aM−1,K−2aM−1,K−1)\begin{equation} A=\left(\begin{array}{ccccc} a

关于寒武纪编程可以参考本人之前的文章。

在上一篇文章中我们介绍了寒武纪显卡实现基本的softmax代码，这里我们借助于寒武纪的流水并行来编写进一步的策略。

寒武纪实现softmax包括下面5个步骤，我们也采取5个kernel来实现softmax：unionMaxKernel(float* middle, float* source1, int num)，这个kernel使用的任务类型是union1，其中middle的长度为taskDim，达到的目的是middle对应元素存储的是不同taskId处理的那部分数据的局部max。

Queue是用来管理并行操作的一种方式。在同一个Queue中所下发的操作顺序是依次按照下发顺序执行的，在不同Queue中的执行操作则是乱序的。Queue所支持的操作不仅包括Kernel（在设备上的执行程序），还包括Notifier、内存拷贝、设备与主机的数据传输等一系列操作，所有需要按序执行的操作都可以下发至同一个Queue，Queue在当前操作满足执行条件后会依次按照进入顺序执行操作；而不同Queue的操作则不保证顺序，这些操作的顺序是无法预测的。

前面博客的实现逻辑是：以blockIdx.x作为矩阵Q的第i行输入，blockIdx.y × blockDim.y + threadIdx.y作为矩阵V的第j行输入，中间遍历矩阵K,V时使用for循环串行遍历。在for循环里面相当于每次只能获得QK.T的一个元素，我们的优化策略是以threadIdx.y来加速QK.T的矩阵乘法。这种做法的缺点是：如果N很大，d很小，那么会造成大量的线程浪费闲置，尤其是当d=1的时候，此时算法接近串行。

CαopAopBβCCαopAopB))βC，这里的运算符op决定是否转置。这里首先说明，cublasSgemm专门处理float类型的二维矩阵乘法，针对其他类型的矩阵，cublas有类似的函数接口cublasDgemm：处理double类型矩阵cublasCgemm：处理cuComplex类型矩阵cublasZgemm：处理cuDoubleComplex类型矩阵cublasHgemm：处理__half类型的矩阵。

本文参考了，其中oneflow实现softmax的CUDA编程源代码参考链接关于softmax的解读以及CUDA代码实现可以参考本人之前编写的几篇文章。

代码22行这里的i=blockIdx.x表示的输出矩阵的第i行，也是输入矩阵Q的第i行，代码23行的phd表示输出矩阵O的第phd列，也是输入矩阵V的第phd列。这两个索引在全过程不参与运算，只参与标记元素的位置。现在还剩下threadIdx.x没有派上用场，我们将利用threadIdx.x尝试遍历矩阵K的每一行，但是BLOCK_DIM_x是固定的，因此我们需要借助一个额外的循环来遍历K，因此phNumN就是总共要遍历的次数，这里一定要向上取整。

本文参考了链接添加链接描述flash attention介绍flash attention的介绍可以参考论文：FlashAttention: Fast and Memory-Efficient Exact Attentionwith IO-Awareness，具体的数学公式参考下面这个图片：其中注意关于矩阵S有两个维度，softmax的操作维度是dim=1，用pytorch表示就是torch.softmax(S, dim=1)对于flash attention来说，里面有两次矩阵乘法，对于这样的二

在上一篇博客我们介绍了softmax的内容以及相关的编程实现，总结一下softmax的特点如下：（没有特殊说明的情况下，考虑的仍然只是1D向量x1：获得向量的全局最大值M2：针对向量xexpxi​−Mi​求和得到S3：输出向量yxS从常理来看，我们至少需要发动3次kernel才能获得数值结果，但是这样非常耗时，有没有什么办法可以发动一次kernel就获得结果呢，下面不妨看看我们的分析过程。

expand和where介绍当谈到 Torch 中的 expand 函数时，我们实际上是指 PyTorch（Torch 的 Python 接口）中的 expand 方法。下面是对 expand 方法和 where 函数的介绍，包括它们的输入和输出：expand 方法：torch.Tensor.expand() 是 PyTorch 中 Tensor 类的一个方法，用于扩展张量的维度。输入：input 是要扩展的张量，size 是一个元组，指定了要扩展的每个维度的大小。输出：返回一个新的张量，形状是

这里写自定义目录标题多重网格算法介绍问题描述——五点差分法求解二维泊松方程五点差分法Gauss迭代算法限制算子介绍提升算子二重网格算法多重网格算法多重网格cuda代码编写串行代码mg.c两重网格cuda并行代码jacobi迭代的cuda编程device_jacobiMakefilecuda_mg.cucuda_mg.hmain.c多重网格算法介绍问题描述——五点差分法求解二维泊松方程{−Δu=f, in Ω,u=g, on ∂Ω,\left\{\begin{al

{ut+uux−(0.01/π)uxx=0,x∈[−1,1],t∈[0,0.99],u(0,x)=−sin⁡(πx),u(t,−1)=u(t,1)=0.\left\{\begin{array}{l}u_t+u u_x-(0.01 / \pi) u_{x x}=0, \quad x \in[-1,1], \quad t \in[0,0.99], \\ u(0, x)=-\sin (\pi x), \\ u(t,-1)=u(t, 1)=0 .\end{array}\right.⎩⎨⎧​ut​+uux​−(0.0

−Δuu​finΩgon∂Ω​这里我们选取Ω012。

1：Error: CUDA driver version is insufficient for CUDA runtime version，这个问题一般是任务在登录节点运行出现的，如果放到计算节点（slurm提交）就会消失。2：__global__开头的函数没有被GPU调用，这种情况可能是slurm文件有问题，比如说没有指定分区gpu（许多人直接拿openmp的slurm用，但是openmp指定分区cpu）。3：所有代码都正确，但是结果有问题，经过考察发现这种问题一般是设置的blockdim=(NX,N

添加链接描述本人在这篇博客已经介绍过稀疏矩阵向量乘法以及相关的一些存储方式，这里重点开始介绍openmp的并行优化。下面这个图片就是采用CSR格式读取稀疏矩阵向量乘法的主体架构，全称为Compressed Sparse Row Matrix压缩稀疏矩阵行格式，该格式对矩阵进行行压缩，通过row_ptr,col_ind,values来确定矩阵。对于第i行来说，该行非零元素的列索引存储在col_ind[row_ptr[i]:row_ptr[i + 1]]中，该行非零元素值存储在values[row_ptr[i

−Δu=f,x∈Ω=(0,1)2−Δu=f,x∈Ω=(0,1)2u=g,x∈∂Ωu=g,x∈∂Ωu=exp⁡(−(x−a)2+(y−b)22l2),a,b,lu=exp(−2l2(x−a)2+(y−b)2​),a,b,l是参数，其中a∈[0,1],b∈[0,1],l∈。

稀疏矩阵主要指的是存在大量零元素的矩阵，常见的稀疏矩阵比如说社交媒体的邻接矩阵，微分方程数值解钟差分法形成的矩阵和有限元方法搭建的刚度矩阵。对于稀疏矩阵，如果对于矩阵每个元素都分配内存存储，将会造成大量的内存浪费，同时做矩阵运算的时候由于读取过程中反复读取零元素将会增大内存访问的时间。因此，为了降低存储空间同时提高访问效率，一种高效的稀疏矩阵存储方法势在必行，下面主要介绍两种稀疏矩阵乘法，我们以一个5×55 \times 55×5的稀疏矩阵存储来举例子。

⇒⇒CPU端定义变量⇒⇒将数据加载到GPU⇒⇒在GPU端执行计算⇒⇒将计算结果拷贝回CPU.CPU端的执行需要一个.c文件来编写代码.GPU端的执行需要一个.cu文件来编写代码.由于C++不能直接调用cuda kernel，还需要编写一个.h文件辅助.即对于一个一般的计算任务，CUDA编程至少应该包含.c文件，.cu文件和.h三个文件。CPU和GPU。

这里先放Makefile和run.slurmMakefilerun.slurm使用阻塞通信实现不同进程之间的消息传递，rank=0: send,recvrank=1: recv,send, send,recv…:recv,send,send,recvrank = size - 1:recv,send上面涉及到相邻进程的消息互换，所以为了方便起见，我们也可以使用send_recv函数来完成相邻进程之间的消息传递，但是另一方面我们发现上面这种写法其实需要对rank进行分类，中间的rank需要两次

直接使用anaconda，下面红色位置使用命令：pip install mpi4py即可，第一次写代码使用导入MPI库，第一次运行代码可能会报错ImportError: DLL load failed点击链接添加链接描述，安装msmpisetup.exe和MS-MPISDK。下面重点介绍MPI的python使用comm = MPI.COMM_WORLD是全局通信器comm.send(msg, dest=1)发送数据，comm.send(data,dest,tag)这是句法。s = comm

MPI,OpenMP的差别run.slurmMakefile九点差分法的矩阵形式omp_orphaning.comp_end.comp_simd_end.c编译这个omp_simd_end.c的时候，需要先module load gcc，这步是为了使用最新版本的gcc，然后需要在Makefile里面写入相关的向量化命令申明函数也需要利用#pragma omp declare simd //申明向量化函数#-----------------------------------

OpenMP的调度主要包括static,dynamic,guided,runtime,auto这5种，下面我们重点介绍前面4种调度方式。我们将针对九点差分法求解泊松方程来考虑这几种调度方式的使用。使用的从句是schedule(type,chunk)Makefilerun.slurmomp_schedule.c**这里注意schedule从句和ordered从句只能跟#pragma omp parallel for，不能跟#pragma omp parallel，**为了方便后面打印不同线程的分配任务

MPI的特点：每个进程都是独立的处理器，有自己私有的内存空间，不同进程之间只能通过消息发送接收来传递信息OpenMP的程序有一系列线程（thread)控制，其中一个进程可以开启多个线程，每个线程都有自己私有的变量，并且线程之间有共享内存，可以通过共享内存交互信息。我们以最简单的代码来对OpenMP做一个基本的介绍omp_hello.c运行代码首先需要设置环境变量，这里假设使用4个线程，-lm连接数学库，不过对于这个代码来说不需要使用。下面我们来看OpenMP代码的基本组成。1：加载头文件#incl

给定一个泊松方程如下,{−Δu=f in Ω,u=u0 on ∂Ω.\left\{\begin{aligned}-\Delta u &= f &&\text{ in } \Omega, \\ u &=u_0 &&\text { on } \partial \Omega.\end{aligned}\right.{−Δuu​=f=u0​​​ in Ω, on ∂Ω.​其中Ω=[0,1]2\Omega = [0,1]^2Ω=[0,1]2，对于该泊松方程，往往可以采取差分法求出计算区域离散点的近似解，具体做

这个代码我们来观察，如果rank=0，那么这个0号进程就会把token发送给1号进程，如果rank=1，那么1号进程就会接收来自0号进程的token，同时把token发送回去给0号进程，只有当rank=0的时候才会打印结果。这个代码很简单，但是如果0号进程把token发送出去，1号进程没有接收到token，此时1号进程是否会启动自己的发送命令呢？

#include <iostream>#include <cmath>#include <ctime>using namespace std;#define PI M_PIconst int M = 11;const int N = 11;const double area[2][2] = {{0,1},{0,1}};double UU(double x,double y,int order[2]){ if(order[0] == 0 &a

GPU训练前提：电脑里得有显卡图片最底下那个NVDIA就是我们需要的显卡。我们仍然以求解泊松方程为例，这里不多说，只是单纯给出在哪些地方需要使用cuda。pinn在损失函数Lossf,error这两个函数定义中需要把数据传输到设备cuda上。这里特别要注意torch.autograd.grad函数的使用。只放代码：import torchimport timeimport numpy as npimport matplotlib.pyplot as pltimport torch.n

问题描述完成按列求和规约程序的并行优化，主要要求：∙\bullet∙只能使用一块Tesla K40m GPU进行计算。∙\bullet∙测试附件slurm中几种典型规模的加速效果。原始代码解析：代码运行：进入文件夹以后，我们先加载cuda模块，这里特别注意，我们需要使用命令：module add cuda/9.0加载9.0版本的cuda。然后再使用make编译，如果事先直接使用make或者加载错误版本的cuda，然后编译，这样会生成一些中间文件，这个时候编译报错，此时重新加载9.0版本的cu

问题描述求解微分方程：φxx+φyy−αφ=F.\varphi_{xx} + \varphi_{yy} - \alpha \varphi = F.φxx​+φyy​−αφ=F.其中测试问题的时候选取φ=(1−x2)(1−y2)\varphi = (1 - x^2)(1 - y^2)φ=(1−x2)(1−y2)，计算区域是[−1,1]2[-1,1]^2[−1,1]2使用五点差分格式来求解问题:φi−1,j−2φi,j+φi+1,jh2+φi,j−1−2φi,j+φi,j+1h2−αφi,j=Fi,j.

问题描述给定一个数据集，里面存放着BATCH_{}SIZE个L1×L2L_1 \times L_2L1​×L2​大小的矩阵，在给出的代码中，BATCH_SIZE=128,L1=512,L2=512BATCH\_{}SIZE=128,L_1=512,L_2=512BATCH_SIZE=128,L1​=512,L2​=512。再给出一个3维张量kernel，其中kernel的形状为BATCH_SIZE×KSIZE×KSIZEBATCH\_{}SIZE\times KSIZE \times KSIZEBATCH

区域分解算法以下都以泊松方程来举例问题引入给定二维空间的泊松方程：−Δu=f,x∈Ω-\Delta u = f,x \in \Omega−Δu=f,x∈Ωu=g,x∈∂Ωu = g,x \in \partial \Omegau=g,x∈∂Ω我们把区域分成若干部分分别求解，然后相邻区域做一个通信，这样就可以得到一个并行算法。下面我们先把区域分成两部分做一个介绍：考虑区域Ω1,Ω2\Omega_1,\Omega_2Ω1​,Ω2​其中红色的是Ω1\Omega_1Ω1​，蓝色的是Ω2\Omega

方程引出y′=Ay+F=g(t,y),y(0)=[00]y^{'}=Ay+F=g(t,y),y(0)=\left[\begin{array}{cc}0\\0\end{array}\right]y′=Ay+F=g(t,y),y(0)=[00​]A=[01−10]A=\left[\begin{array}{cc}0&1\\-1&0\end{array}\right]A=[0−1​10​]F=[0t]F=\left[\begin{array}{cc}0\\t\end{array}\right

微分方程引入−u′′+ρ(u)=f,x∈(0,1)-u^{''}+\rho\sqrt(u)=f,x \in(0,1)−u′′+ρ(​u)=f,x∈(0,1)u(0)=α,u(1)=βu(0)=\alpha,u(1)=\betau(0)=α,u(1)=β我们假设剖分为N，那么就利用中心差分得到：h=1N−1h=\frac{1}{N-1}h=N−11​−ui−1+2ui−ui+1−(R(ui)+fi)h2=0,1≤i≤N−2-u_{i-1}+2u_i-u_{i+1}-(R(u_i)+f_i)h^2=0

非线性方程的求解引入非线性方程：F(X)=[2α(x0−1)+4x0(x02+x12−0.25)2α(x1−1)+4x1(x02+x12−0.25)]F(X)=\left[\begin{array}{c}2\alpha(x_0-1)+4x_0(x_0^2+x_1^2-0.25)\\2\alpha(x_1-1)+4x_1(x_0^2+x_1^2-0.25) \end{array}\right]F(X)=[2α(x0​−1)+4x0​(x02​+x12​−0.25)2α(x1​−1)+4x1​(x02​+x

第一种方法只需要利用前面提到的矩阵和向量搭建前面几篇博客提到了矩阵，向量的搭建，本人也介绍来求解任意大小的线性系统的方法，下面本人先介绍不调用Petsc任何有关的包的求解方法，为此先引入基本定义：−Δu=f,x∈Ω=(0,1)2-\Delta u = f,x \in \Omega=(0,1)^2−Δu=f,x∈Ω=(0,1)2u=u0(x,y),(x,y)∈∂Ωu=u_{0}(x,y),(x,y)\in \partial \Omegau=u0​(x,y),(x,y)∈∂Ω这里取f=4∗cos(2x)