CPU-GPU异构并行化APSP算法

一、Floyd-Warshall算法

介绍

Floyd-Warshall算法（英语：Floyd-Warshall algorithm），中文亦称弗洛伊德算法或佛洛依德算法，是解决任意两点间的最短路径的一种算法，可以正确处理有向图或负权（但不可存在负权回路）的最短路径问题，同时也被用于计算有向图的闭包传递。

原理

其本质为动态规划，给定有向图图 G = ( V , E ) G = (V, E) G=(V,E)，其中 V ( v e r t i c e s ) V(vertices) V(vertices)为顶点数， E ( e d g e s ) E(edges) E(edges)为边数，并给出初始权重矩阵 w [ i ] [ j ] w[i][j] w[i][j]，表示顶点 i → j i \rightarrow j i→j的权重，其表达式为：
w i , j = { weight of edge ( i , j ) if ( i , j ) ∈ E ; ∞ if ( i , j ) ∉ E . \left.w_{i,j}=\left\{\begin{array}{ll}\text{weight of edge}\left(i,j\right)&\text{if}\left(i,j\right)\in E;\\\infty&\text{if}\left(i,j\right)\notin E.\end{array}\right.\right. wi,j​={weight of edge(i,j)∞​if(i,j)∈E;if(i,j)∈/E.​
即，对于 i → j i \rightarrow j i→j未连通的边通常设置为一个无穷大的数 I N F INF INF；对于动态规划算法需要定义状态 D i , j , k D_{i,j,k} Di,j,k​:从 i i i到 j j j只以( 1.. k 1..k 1..k)集合中的节点为中间节点的最短路径的长度；则可分为以下2种情况讨论：

  1. 如果最短路经过点 k k k： D i , j , k = D i , k , k − 1 + D k , j , k − 1 . D_{i,j,k}=D_{i,k,k-1}+D_{k,j,k-1}. Di,j,k​=Di,k,k−1​+Dk,j,k−1​.

  2. 若最短路径不经过点 k k k: D i , j , k = D i , j , k − 1  。  D_{i,j,k}=D_{i,j,k-1\text{ 。 }} Di,j,k​=Di,j,k−1 。 ​

若不能理解 k − 1 k - 1 k−1的含义，则可理解为下一层 k k k的状态需要上一层 k − 1 k - 1 k−1推导出(因为要逐个枚举中间节点，例如 D 1 , 3 = D 1 , 2 + D 2 , 3 D_{1,3} = D_{1,2} + D_{2,3} D1,3​=D1,2​+D2,3​，那么需要保证 D 1 , 2 , D 2 , 3 D_{1,2},D_{2,3} D1,2​,D2,3​是对应的最短距离，才能导致 D 1 , 3 D_{1,3} D1,3​是1号节点到3号节点的最短距离)即第 k k k层状态依赖于第 k − 1 k-1 k−1层状态，故不可对 k k k层循环做并行化处理；最后可以得到状态转移方程:
D i , j , k = min ⁡ ( D i , j , k − 1 , D i , k , k − 1 + D k , j , k − 1 ) D_{i,j,k}=\min(D_{i,j,k-1},D_{i,k,k-1}+D_{k,j,k-1}) Di,j,k​=min(Di,j,k−1​,Di,k,k−1​+Dk,j,k−1​)
在实际算法中，为了节约空间，可以直接在原来空间上进行迭代，这样空间可降至二维。

分析

  1. 时间复杂度： O ( V 3 ) O(V^3) O(V3)，其中 V V V是点集，对于 i , j i,j i,j两层for循环可使用OpenMP优化到线性
  2. 空间复杂度： O ( V 2 ) O(V^2) O(V2)

二、CPU-GPU并行化Floyd-APSP算法

为了求到全部的最短路径，不仅需要计算最短路径距离矩阵 D D D，还需要计算最短路径构造矩阵 C C C。其中 C C C矩阵的定义为：如果在顶点 i i i和顶点 j j j之间至少存在一条最短路径，则 C i , j C_{i,j} Ci,j​表示最短路径上编号最高的中间顶点，否则为undefined (NULL)。构造矩阵的初值都是未定义的，用数学表示如下：
c i , j ( k ) = { NULL i f   k = 0 ; k i f   k ≥ 1   a n d   d i , j ( k − 1 ) > d i , k ( k − 1 ) + d k , j ( k − 1 ) ; c i , j ( k − 1 ) otherwise. , \left.c_{i,j}^{(k)}=\left\{\begin{array}{ll}\text{NULL}&\mathrm{if~}k=0;\\k&\mathrm{if~}k\geq1\mathrm{~and~}d_{i,j}^{(k-1)}>d_{i,k}^{(k-1)}+d_{k,j}^{(k-1)};\\c_{i,j}^{(k-1)}&\text{otherwise.}\end{array}\right.\right., ci,j(k)​=⎩ ⎨ ⎧​NULLkci,j(k−1)​​if k=0;if k≥1 and di,j(k−1)​>di,k(k−1)​+dk,j(k−1)​;otherwise.​,
其中 C i , j k − 1 C_{i,j}^{k-1} Ci,jk−1​与上相同，由于下一层受到上一层的制约，为递推关系。

Algorithm1: Floyd-Warshall

  • Floyd-Warshall算法用于计算最短路径距离矩阵 D i , j D_{i,j} Di,j​和最短路径构造矩阵 C i , j C_{i,j} Ci,j​

在这里插入图片描述

Algorithm2:

  • 输出一对顶点 ( i , j ) (i, j) (i,j)之间最短路径的中间顶点的递归过程

在这里插入图片描述

可以想象为二叉树，一边是往左子树遍历，一边是往右子树遍历，即左根右的中序遍历。

分块联合算法

  • 该算法是为在CPU-GPU混合系统上实现高GPU利用率的快速APSP解决方案而设计的。

在分块联合算法中，将 n × n n \times n n×n的距离矩阵 D i , j D_{i,j} Di,j​和构造矩阵 C i , j C_{i,j} Ci,j​划分为 b × b b \times b b×b的子矩阵的分块，其中 b b b为分块因子，为以下问题讨论方便，假设 n % b = = 0 n \% b ==0 n%b==0，即 n n n能整除 b b b，并在每个块内有定义 A I , J = a [ i , j ] A_{I, J} = a[i, j] AI,J​=a[i,j]来标识块索引为 ( I , J ) (I,J) (I,J)的子矩阵，用数学符号表示为：
1 ≤ I , J ≤ [ n b ] , 1 ≤ i , j ≤ b 1 \leq I, J \leq [\frac{n}{b}] , \\ 1 \leq i,j \leq b 1≤I,J≤[bn​],1≤i,j≤b
如下图所示，展现了 n = 12 n = 12 n=12矩阵的示例，其中 b = 3 b=3 b=3

在这里插入图片描述

Algorithm3
  • 针对APSP问题的分块联合算法

在这里插入图片描述

将该算法划分4个阶段为：

  1. 首先将 n × n n \times n n×n的矩阵分解为长度为 [ n b ] × [ n b ] [\frac{n}{b}] \times [\frac{n}{b}] [bn​]×[bn​]的以 b × b b \times b b×b的矩阵，并外层枚举节点 ( K , K ) (K, K) (K,K)，其中 1 ≤ K ≤ [ n b ] 1 \leq K \leq [\frac{n}{b}] 1≤K≤[bn​]，并在子矩阵 b × b b \times b b×b中使用Floyd-WarShall方法，求解 D K , K , C K , K D_{K, K}, C_{K,K} DK,K​,CK,K​。
  2. 对节点 ( K , K ) (K, K) (K,K)所在的第 K K K列进行MMA即矩阵乘法加法操作
  3. 对节点 ( K , K ) (K, K) (K,K)所在的第 K K K行进行MMA即矩阵乘法加法操作
  4. 对于除以上涉及到的剩余节点
Algorithm4
  • APSP子问题的阻塞联合算法
  • 区别在于：算法4的4-16行运行在GPU中，算法4的合并操作17-20运行在CPU中。

在这里插入图片描述

在这里插入图片描述

Algorithm5

  • 子分块联合APSP的矩阵-矩阵""乘-加"算法

  • 代数中的MMA算法可以扩展为同时计算路径矩阵 D i , j D_{i,j} Di,j​和构造矩阵 C i , j C_{i,j} Ci,j​

z i , j ← min ⁡ ( z i , j , ∑ k = 1 b x i , k + y k , j ) z_{i,j} \leftarrow \min(z_{i,j}, \sum_{k=1}^b x_{i,k}+y_{k,j}) zi,j​←min(zi,j​,k=1∑b​xi,k​+yk,j​)

其中， z i , j z_{i,j} zi,j​为 b × b b \times b b×b的子矩阵。

在这里插入图片描述

Algorithm6

  • 阶段2-阶段4可以使用矩阵乘法更新，在本问题中，就是极小加代数。极小加代数的乘法和加法是分离执行的，极小加操作(MINPLUS)是运行在GPU中，矩阵加(MMA)运行在CPU中。

  • 这个操作减少了 Z ， C Z，C Z，C从CPU到GPU的数据传输，也就允许了CPU和GPU之间的高速通信。

在这里插入图片描述

Code

  • 以下代码均来自：https://github.com/EricLu1218/Parallel_Programming

模拟GPU的串行

#include <stdio.h>
#include <stdlib.h>

const int INF = ((1 << 30) - 1);
const int V = 50010;
void input(char *inFileName);
void output(char *outFileName);

void block_FW(int B);
int ceil(int a, int b);
void cal(int B, int Round, int block_start_x, int block_start_y, int block_width, int block_height);

int n, m;
static int Dist[V][V];

int main(int argc, char *argv[])
{
    input(argv[1]);
    int B = 512;
    block_FW(B);
    output(argv[2]);
    return 0;
}

void input(char *infile)
{
    FILE *file = fopen(infile, "rb");
    fread(&n, sizeof(int), 1, file);
    fread(&m, sizeof(int), 1, file);

    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            if (i == j)
            {
                Dist[i][j] = 0;
            }
            else
            {
                Dist[i][j] = INF;
            }
        }
    }

    int pair[3];
    for (int i = 0; i < m; ++i)
    {
        fread(pair, sizeof(int), 3, file);
        Dist[pair[0]][pair[1]] = pair[2];
    }
    fclose(file);
}

void output(char *outFileName)
{
    FILE *outfile = fopen(outFileName, "w");
    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < n; ++j)
        {
            if (Dist[i][j] >= INF)
                Dist[i][j] = INF;
        }
        fwrite(Dist[i], sizeof(int), n, outfile);
    }
    fclose(outfile);
}

int ceil(int a, int b) { return (a + b - 1) / b; }

void block_FW(int B)
{
    int round = ceil(n, B);
    for (int r = 0; r < round; ++r)
    {
        printf("%d %d\n", r, round);
        fflush(stdout);
        /* Phase 1 */
        cal(B, r, r, r, 1, 1);

        /* Phase 2 */
        cal(B, r, r, 0, r, 1);
        cal(B, r, r, r + 1, round - r - 1, 1);
        cal(B, r, 0, r, 1, r);
        cal(B, r, r + 1, r, 1, round - r - 1);

        /* Phase 3 */
        cal(B, r, 0, 0, r, r);
        cal(B, r, 0, r + 1, round - r - 1, r);
        cal(B, r, r + 1, 0, r, round - r - 1);
        cal(B, r, r + 1, r + 1, round - r - 1, round - r - 1);
    }
}

void cal(
    int B, int Round, int block_start_x, int block_start_y, int block_width, int block_height)
{
    int block_end_x = block_start_x + block_height;
    int block_end_y = block_start_y + block_width;

    for (int b_i = block_start_x; b_i < block_end_x; ++b_i)
    {
        for (int b_j = block_start_y; b_j < block_end_y; ++b_j)
        {
            // To calculate B*B elements in the block (b_i, b_j)
            // For each block, it need to compute B times
            for (int k = Round * B; k < (Round + 1) * B && k < n; ++k)
            {
                // To calculate original index of elements in the block (b_i, b_j)
                // For instance, original index of (0,0) in block (1,2) is (2,5) for V=6,B=2
                int block_internal_start_x = b_i * B;
                int block_internal_end_x = (b_i + 1) * B;
                int block_internal_start_y = b_j * B;
                int block_internal_end_y = (b_j + 1) * B;

                if (block_internal_end_x > n)
                    block_internal_end_x = n;
                if (block_internal_end_y > n)
                    block_internal_end_y = n;

                for (int i = block_internal_start_x; i < block_internal_end_x; ++i)
                {
                    for (int j = block_internal_start_y; j < block_internal_end_y; ++j)
                    {
                        if (Dist[i][k] + Dist[k][j] < Dist[i][j])
                        {
                            Dist[i][j] = Dist[i][k] + Dist[k][j];
                        }
                    }
                }
            }
        }
    }
}

单GPU的CUDA代码

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

const int INF = (1 << 30) - 1;
int vertex_num, edge_num, matrix_size;
int *dist;

double cal_time(struct timespec start, struct timespec end)
{
	struct timespec temp;
	if ((end.tv_nsec - start.tv_nsec) < 0)
	{
		temp.tv_sec = end.tv_sec - start.tv_sec - 1;
		temp.tv_nsec = 1000000000 + end.tv_nsec - start.tv_nsec;
	}
	else
	{
		temp.tv_sec = end.tv_sec - start.tv_sec;
		temp.tv_nsec = end.tv_nsec - start.tv_nsec;
	}
	return temp.tv_sec + (double)temp.tv_nsec / 1000000000.0;
}

__device__ __host__ size_t index_convert(int i, int j, int row_size)
{
	return i * row_size + j;
}

void input(char *input_file_path, int &block_factor)
{
	FILE *input_file = fopen(input_file_path, "rb");
	fread(&vertex_num, sizeof(int), 1, input_file);
	fread(&edge_num, sizeof(int), 1, input_file);

	matrix_size = ceil((double)vertex_num / (double)block_factor) * block_factor;
	cudaMallocHost((void **)&dist, matrix_size * matrix_size * sizeof(int));
	for (int i = 0; i < matrix_size; ++i)
	{
		for (int j = 0; j < matrix_size; ++j)
		{
			if (i != j)
				dist[index_convert(i, j, matrix_size)] = INF;
			else if (i < vertex_num)
				dist[index_convert(i, j, matrix_size)] = 0;
			else
				dist[index_convert(i, j, matrix_size)] = INF;
		}
	}

	int data[3];
	for (int i = 0; i < edge_num; ++i)
	{
		fread(data, sizeof(int), 3, input_file);
		dist[index_convert(data[0], data[1], matrix_size)] = data[2];
	}
	fclose(input_file);
}

void output(char *output_file_path)
{
	FILE *output_file = fopen(output_file_path, "w");
	for (int i = 0; i < vertex_num; ++i)
	{
		fwrite(&dist[index_convert(i, 0, matrix_size)], sizeof(int), vertex_num, output_file);
	}
	fclose(output_file);
}

__constant__ int size[3]; //matrix size, block_factor, grid_size

__global__ void phase1(int *d_dist, int round)
{
	__shared__ int share[4 * 1024];
	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset = size[1] * round;
	int j_offset = size[1] * round;

	share[index_convert(j, i, size[1])] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
#pragma unroll 32
	for (int k = 0; k < size[1]; ++k)
	{
		__syncthreads();
		if (share[index_convert(j, i, size[1])] > share[index_convert(j, k, size[1])] + share[index_convert(k, i, size[1])])
			share[index_convert(j, i, size[1])] = share[index_convert(j, k, size[1])] + share[index_convert(k, i, size[1])];
	}
	d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = share[index_convert(j, i, size[1])];
}

__global__ void phase2(int *d_dist, int round)
{
	__shared__ int share[3 * 4 * 1024];
	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset, j_offset;
	if (blockIdx.x == 0)
	{
		i_offset = size[1] * ((round + blockIdx.y + 1) % size[2]);
		j_offset = size[1] * round;
		share[index_convert(i, j, size[1])] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
		share[index_convert(i + size[1], j, size[1])] = share[index_convert(i, j, size[1])];
		share[index_convert(i + 2 * size[1], j, size[1])] = d_dist[index_convert(j_offset + i, j_offset + j, size[0])];
	}
	else
	{
		i_offset = size[1] * round;
		j_offset = size[1] * ((round + blockIdx.y + 1) % size[2]);
		share[index_convert(i, j, size[1])] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
		share[index_convert(i + size[1], j, size[1])] = d_dist[index_convert(i_offset + i, i_offset + j, size[0])];
		share[index_convert(i + 2 * size[1], j, size[1])] = share[index_convert(i, j, size[1])];
	}

#pragma unroll 32
	for (int k = 0; k < size[1]; ++k)
	{
		__syncthreads();
		if (share[index_convert(i, j, size[1])] >
			share[index_convert(i + size[1], k, size[1])] + share[index_convert(k + 2 * size[1], j, size[1])])
			share[index_convert(i, j, size[1])] =
				share[index_convert(i + size[1], k, size[1])] + share[index_convert(k + 2 * size[1], j, size[1])];
	}
	d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = share[index_convert(i, j, size[1])];
}

__global__ void phase3(int *d_dist, int round)
{
	__shared__ int share[3 * 4 * 1024];
	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset = size[1] * ((round + blockIdx.y + 1) % size[2]);
	int j_offset = size[1] * ((round + blockIdx.x + 1) % size[2]);
	int r_offset = size[1] * round;

	share[index_convert(i, j, size[1])] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
	share[index_convert(i + size[1], j, size[1])] = d_dist[index_convert(i_offset + i, r_offset + j, size[0])];
	share[index_convert(i + 2 * size[1], j, size[1])] = d_dist[index_convert(r_offset + i, j_offset + j, size[0])];
#pragma unroll 32
	for (int k = 0; k < size[1]; ++k)
	{
		__syncthreads();
		if (share[index_convert(i, j, size[1])] >
			share[index_convert(i + size[1], k, size[1])] + share[index_convert(k + 2 * size[1], j, size[1])])
			share[index_convert(i, j, size[1])] =
				share[index_convert(i + size[1], k, size[1])] + share[index_convert(k + 2 * size[1], j, size[1])];
	}
	d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = share[index_convert(i, j, size[1])];
}

int main(int argc, char **argv)
{
	double total_time, bfd_time;
	timespec total_time1, total_time2, bfd_time1, bfd_time2;

	clock_gettime(CLOCK_MONOTONIC, &total_time1);
	cudaSetDevice(0); // 设置运行的为第0块GPU
	int block_factor = 32;
	if (argc == 4)
		block_factor = atoi(argv[3]);
	input(argv[1], block_factor); // 读取数据并初始化dist
	int grid_size = matrix_size / block_factor; // 划分后的网格大小N = [n / b]

	int size_info[3] = {matrix_size, block_factor, grid_size}; // n, b, N = [n / b]
	cudaMemcpyToSymbol(size, size_info, 3 * sizeof(int)); //  将矩阵大小、块大小和网格大小的信息传递给CUDA设备

	int *d_dist;
	clock_gettime(CLOCK_MONOTONIC, &bfd_time1);
	cudaMalloc(&d_dist, (size_t)sizeof(int) * matrix_size * matrix_size); // 在GPU上分配内存
	// 在GPU上分配和复制内存，将距离矩阵dist从主机（CPU）内存拷贝到设备（GPU）内存
	cudaMemcpy(d_dist, dist, (size_t)sizeof(int) * matrix_size * matrix_size, cudaMemcpyHostToDevice);
	// 定义了CUDA的线程块和网格的维度
	dim3 block(block_factor, block_factor); // (b, b)
	dim3 grid2(2, grid_size - 1); // (2, N - 1)
	dim3 grid3(grid_size - 1, grid_size - 1); // (N - 1, N - 1)
	for (int r = 0; r < grid_size; ++r)
	{
		phase1<<<1, block>>>(d_dist, r);
		phase2<<<grid2, block>>>(d_dist, r);
		phase3<<<grid3, block>>>(d_dist, r);
	}
	cudaMemcpy(dist, d_dist, (size_t)sizeof(int) * matrix_size * matrix_size, cudaMemcpyDeviceToHost);
	clock_gettime(CLOCK_MONOTONIC, &bfd_time2);

	output(argv[2]);
	cudaFree(d_dist);
	cudaFree(dist);

	clock_gettime(CLOCK_MONOTONIC, &total_time2);
	bfd_time = cal_time(bfd_time1, bfd_time2);
	total_time = cal_time(total_time1, total_time2);
	printf(" vertex:   %d\n", vertex_num);
	printf(" I/O time: %.5f\n", total_time - bfd_time);
	printf(" cal time: %.5f\n", bfd_time);
	printf(" runtime:  %.5f\n", total_time);
	return 0;
}

2个GPU代码

#include <math.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

const int INF = (1 << 30) - 1;
int vertex_num, edge_num, matrix_size;
int *dist;

double cal_time(struct timespec start, struct timespec end)
{
	struct timespec temp;
	if ((end.tv_nsec - start.tv_nsec) < 0)
	{
		temp.tv_sec = end.tv_sec - start.tv_sec - 1;
		temp.tv_nsec = 1000000000 + end.tv_nsec - start.tv_nsec;
	}
	else
	{
		temp.tv_sec = end.tv_sec - start.tv_sec;
		temp.tv_nsec = end.tv_nsec - start.tv_nsec;
	}
	return temp.tv_sec + (double)temp.tv_nsec / 1000000000.0;
}

__device__ __host__ size_t index_convert(int i, int j, int row_size)
{
	return i * row_size + j;
}

void input(char *input_file_path, int block_factor)
{
	FILE *input_file = fopen(input_file_path, "rb");
	fread(&vertex_num, sizeof(int), 1, input_file);
	fread(&edge_num, sizeof(int), 1, input_file);

	matrix_size = ceil((double)vertex_num / (double)block_factor) * block_factor;
	cudaMallocHost((void **)&dist, matrix_size * matrix_size * sizeof(int));
	for (int i = 0; i < matrix_size; ++i)
	{
		for (int j = 0; j < matrix_size; ++j)
		{
			if (i != j)
				dist[index_convert(i, j, matrix_size)] = INF;
			else if (i < vertex_num)
				dist[index_convert(i, j, matrix_size)] = 0;
			else
				dist[index_convert(i, j, matrix_size)] = INF;
		}
	}

	int data[3];
	for (int i = 0; i < edge_num; ++i)
	{
		fread(data, sizeof(int), 3, input_file);
		dist[index_convert(data[0], data[1], matrix_size)] = data[2];
	}
	fclose(input_file);
}

void output(char *output_file_path)
{
	FILE *output_file = fopen(output_file_path, "w");
	for (int i = 0; i < vertex_num; ++i)
	{
		fwrite(&dist[index_convert(i, 0, matrix_size)], sizeof(int), vertex_num, output_file);
	}
	fclose(output_file);
}

__constant__ int size[3]; //matrix size, block_factor, grid_size

__global__ void phase1(int *d_dist, int round)
{
	__shared__ int pivot[1024];
	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset = 32 * round;
	int j_offset = 32 * round;

	pivot[index_convert(i, j, 32)] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
#pragma unroll 32
	for (int k = 0; k < 32; ++k)
	{
		__syncthreads();
		if (pivot[index_convert(i, j, 32)] > pivot[index_convert(i, k, 32)] + pivot[index_convert(k, j, 32)])
			pivot[index_convert(i, j, 32)] = pivot[index_convert(i, k, 32)] + pivot[index_convert(k, j, 32)];
	}
	d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = pivot[index_convert(i, j, 32)];
}

__global__ void phase2(int *d_dist, int round)
{
	__shared__ int self[1024], pivot[1024];
	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset, j_offset;
	if (blockIdx.x == 0 && blockIdx.y != round)
	{
		i_offset = 32 * blockIdx.y;
		j_offset = 32 * round;

		self[index_convert(i, j, 32)] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
		pivot[index_convert(i, j, 32)] = d_dist[index_convert(j_offset + i, j_offset + j, size[0])];
#pragma unroll 32
		for (int k = 0; k < 32; ++k)
		{
			__syncthreads();
			if (self[index_convert(i, j, 32)] > self[index_convert(i, k, 32)] + pivot[index_convert(k, j, 32)])
				self[index_convert(i, j, 32)] = self[index_convert(i, k, 32)] + pivot[index_convert(k, j, 32)];
		}
		d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = self[index_convert(i, j, 32)];
	}
	else if (blockIdx.y != round)
	{
		i_offset = 32 * round;
		j_offset = 32 * blockIdx.y;

		self[index_convert(i, j, 32)] = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
		pivot[index_convert(i, j, 32)] = d_dist[index_convert(i_offset + i, i_offset + j, size[0])];
#pragma unroll 32
		for (int k = 0; k < 32; ++k)
		{
			__syncthreads();
			if (self[index_convert(i, j, 32)] > pivot[index_convert(i, k, 32)] + self[index_convert(k, j, 32)])
				self[index_convert(i, j, 32)] = pivot[index_convert(i, k, 32)] + self[index_convert(k, j, 32)];
		}
		d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = self[index_convert(i, j, 32)];
	}
}

__global__ void phase3(int *d_dist, int round, int grid_offset)
{
	__shared__ int col[1024], row[1024];
	int self;

	int block_i = grid_offset + blockIdx.y;
	int block_j = blockIdx.x;
	if (block_i == round || block_j == round)
		return;

	int i = threadIdx.y;
	int j = threadIdx.x;

	int i_offset = 32 * block_i;
	int j_offset = 32 * block_j;
	int r_offset = 32 * round;

	self = d_dist[index_convert(i_offset + i, j_offset + j, size[0])];
	col[index_convert(i, j, 32)] = d_dist[index_convert(i_offset + i, r_offset + j, size[0])];
	row[index_convert(i, j, 32)] = d_dist[index_convert(r_offset + i, j_offset + j, size[0])];

#pragma unroll 32
	for (int k = 0; k < 32; ++k)
	{
		__syncthreads();
		if (self > col[index_convert(i, k, 32)] + row[index_convert(k, j, 32)])
			self = col[index_convert(i, k, 32)] + row[index_convert(k, j, 32)];
	}
	d_dist[index_convert(i_offset + i, j_offset + j, size[0])] = self;
}

int main(int argc, char **argv)
{
	const int block_factor = 32, device_num = 2;
	input(argv[1], block_factor);
	int grid_size = matrix_size / block_factor;

	int *d_dist[2];
#pragma omp parallel num_threads(device_num)
	{
		int device_id = omp_get_thread_num();
		cudaSetDevice(device_id);

		int size_info[3] = {matrix_size, block_factor, grid_size};
		cudaMemcpyToSymbol(size, size_info, 3 * sizeof(int));

		int grid_partition = grid_size / device_num;
		int grid_offset = device_id * grid_partition;
		int grid_count = grid_partition;
		if (device_id == device_num - 1)
			grid_count += grid_size % device_num;
		size_t grid_start = grid_offset * block_factor * matrix_size;

		cudaMalloc(&(d_dist[device_id]), (size_t)sizeof(int) * matrix_size * matrix_size);
#pragma omp barrier
		cudaMemcpy(&(d_dist[device_id][grid_start]), &(dist[grid_start]),
				   (size_t)sizeof(int) * block_factor * grid_count * matrix_size, cudaMemcpyHostToDevice);
		dim3 block(block_factor, block_factor);
		dim3 grid2(2, grid_size);
		dim3 grid3(grid_size, grid_count);
		for (int r = 0; r < grid_size; ++r)
		{
			if (grid_offset <= r && r < grid_offset + grid_count)
			{
				size_t copy_start = r * block_factor * matrix_size;
				if (device_id == 0)
					cudaMemcpy(&(d_dist[1][copy_start]), &(d_dist[0][copy_start]),
							   (size_t)sizeof(int) * block_factor * matrix_size, cudaMemcpyDeviceToDevice);
				else
					cudaMemcpy(&(d_dist[0][copy_start]), &(d_dist[1][copy_start]),
							   (size_t)sizeof(int) * block_factor * matrix_size, cudaMemcpyDeviceToDevice);
			}
#pragma omp barrier
			phase1<<<1, block>>>(d_dist[device_id], r);
			phase2<<<grid2, block>>>(d_dist[device_id], r);
			phase3<<<grid3, block>>>(d_dist[device_id], r, grid_offset);
		}
		cudaMemcpy(&(dist[grid_start]), &(d_dist[device_id][grid_start]),
				   (size_t)sizeof(int) * block_factor * grid_count * matrix_size, cudaMemcpyDeviceToHost);
		cudaFree(d_dist[omp_get_thread_num()]);
#pragma omp barrier
	}

	output(argv[2]);
	cudaFree(dist);
	return 0;
}

Reference

  1. https://zh.wikipedia.org/wiki/Floyd-Warshall%E7%AE%97%E6%B3%95
  2. Blocked United Algorithm for the All-Pairs Shortest Paths Problem on Hybrid CPU-GPU Systems
  3. https://github.com/EricLu1218/Parallel_Programming

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/388545.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

Github用人工智能（AI）帮你的代码修正安全漏洞

Github用人工智能（AI）帮你的代码修正安全漏洞

每周跟踪AI热点新闻动向和震撼发展 想要探索生成式人工智能的前沿进展吗&#xff1f;订阅我们的简报&#xff0c;深入解析最新的技术突破、实际应用案例和未来的趋势。与全球数同行一同&#xff0c;从行业内部的深度分析和实用指南中受益。不要错过这个机会&#xff0c;成为AI领…
阅读更多...
2024.02.13作业

2024.02.13作业

21. c 22. b 23. b 5先出栈意味着1234都在栈内&#xff0c;此时1不能比2&#xff0c;3先出栈 24. b, c, d: 10, 12, 120 25. 2, 5 26. 数组越界&#xff0c;可能出现段错误 27. 0, 41 28. 1, 320 29. *a *b; *b *a - *b; *a - *b; 30. 0x801005&#xff1b;0x8…
阅读更多...
GPT-4带来的思想火花

GPT-4带来的思想火花

GPT-4能够以其强大的生成能力和广泛的知识储备激发出众多思想火花。它能够在不同的情境下生成新颖的观点、独特的见解和富有创意的解决方案&#xff0c;这不仅有助于用户突破思维定势&#xff0c;还能促进知识与信息在不同领域的交叉融合。 对于研究者而言&#xff0c;GPT-4可能…
阅读更多...
STM32—DHT11温湿度传感器

STM32—DHT11温湿度传感器

文章目录 一.温湿度原理1.1 时序图 二.代码 一.温湿度原理 1.1 时序图 (1).下图一是DHT11总的时序图。 (2).图二对应图一的左边黑色部分&#xff0c;图三对应图一的绿色部分&#xff0c;图四的左部分图对应图一的红色部分&#xff0c;图四的右部分对应图一的黄色部分。 (3)…
阅读更多...
鸿蒙（HarmonyOS）项目方舟框架（ArkUI）之Navigation组件

鸿蒙（HarmonyOS）项目方舟框架（ArkUI）之Navigation组件

鸿蒙&#xff08;HarmonyOS&#xff09;项目方舟框架&#xff08;ArkUI&#xff09;之Navigation组件 一、操作环境 操作系统: Windows 10 专业版、IDE:DevEco Studio 3.1、SDK:HarmonyOS 3.1 二、Navigation组件 鸿蒙&#xff08;HarmonyOS&#xff09;项目方舟框架&#…
阅读更多...
嵌入式CAN通信协议原理(下)

嵌入式CAN通信协议原理(下)

本篇文章结合实际CAN控制器继续介绍协议相关的内容&#xff0c;还有示例讲解。 好了&#xff0c;继续吧&#xff01; 二. STM32 CAN 控制器介绍 STM32 的芯片中具有 bxCAN 控制器 (Basic Extended CAN)&#xff0c;它支持 CAN 协议 2.0A 和 2.0B 标准。 该 CAN 控制器支持最…
阅读更多...
【C语言】socketpair 的系统调用

【C语言】socketpair 的系统调用

一、 Linux 内核 4.19socketpair 的系统调用 SYSCALL_DEFINE4(socketpair, int, family, int, type, int, protocol,int __user *, usockvec) {return __sys_socketpair(family, type, protocol, usockvec); } 这段代码定义了一个名为 socketpair 的系统调用。系统调用是操作…
阅读更多...
JavaScript基础第三天

JavaScript基础第三天

JavaScript 基础第三天 今天我们学习for循环、while循环、终止循环和无限循环。 1. for 循环 1.1. 语法 // 1. 语法格式 // for(起始值; 结束条件; 累加器) { // // 要重复执行的代码 // }1.2. 示例代码 let sum 0; for (let i 0; i < 100; i) {sum i; } alert(&q…
阅读更多...
阿里云服务器4核16G配置报价和CPU内存性能参数表

阿里云服务器4核16G配置报价和CPU内存性能参数表

阿里云4核16G服务器优惠价格ECS云服务器经济型e实例26元1个月、149元半年、79元3个月&#xff0c;4核16G通用算力u1服务器、通用型g7、通用型g8i、AMD通用型g8a、性能增强通用型g8ae、高主频通用型hfg8i、AMD通用型g7a、内存型r7p等均提供4核16G配置。阿里云服务器网aliyunfuwu…
阅读更多...
C#通过重写虚方法实现加、减、乘、除运算 通过多态确定人类的说话行为

C#通过重写虚方法实现加、减、乘、除运算 通过多态确定人类的说话行为

目录 一、涉及到的知识点1 1.虚方法 2.重写方法 3.重写方法与重载方法的区别 4.通过KeyPressEventArgs.KeyChar限制键盘输入的内容 5.if-else if-else嵌套转换为switch case 二、 涉及到的知识点2 1.多态性 2.使用多态性的注意事项 3. 使用虚方法实现多态性 三、实…
阅读更多...
红日靶场2学习

红日靶场2学习

靶场下载来自&#xff1a; http://vulnstack.qiyuanxuetang.net/vuln/detail/3/ 靶场统一登录密码&#xff1a;1qazWSX 按大佬的说法是 环境需要模拟内网和外网两个网段&#xff0c;PC端虚拟机相当于网关服务器&#xff0c;所以需要两张网卡&#xff0c;一个用来向外网提供web…
阅读更多...
FT2232调试记录（3）

FT2232调试记录（3）

FT2232调试记录&#xff08;1&#xff09;: FT2232调试记录&#xff08;2&#xff09;: FT2232调试记录&#xff08;3&#xff09;: FT2232 SPI读写函数: 参照SPI提供的文档&#xff1a; 工程&#xff1a; SPI 写函数&#xff1a; FT_STATUS write_byte(FT_HANDLE handle…
阅读更多...
剑指offer——数值的整数次方

剑指offer——数值的整数次方

目录 1. 题目描述2. 一般思路2.1 有问题的思路2.2 全面但不高效的思路2.3 面试小提示 3. 全面又高效的思路 1. 题目描述 题目:实现函数 double Power(double base,int exponent)&#xff0c;求base 的exponent 次方。不得使用库函数&#xff0c;同时不需要考虑大数问题 2. 一般…
阅读更多...
4核16g云服务器多少钱？

4核16g云服务器多少钱？

4核16G服务器租用优惠价格26元1个月&#xff0c;腾讯云轻量4核16G12M服务器32元1个月、96元3个月、156元6个月、312元一年&#xff0c;阿腾云atengyun.com分享4核16服务器租用费用价格表&#xff0c;阿里云和腾讯云详细配置报价和性能参数表&#xff1a; 腾讯云4核16G服务器价…
阅读更多...
VSCode汉化以及跳转到函数定义设置

VSCode汉化以及跳转到函数定义设置

目录 VSCode汉化设置跳转到函数定义 VSCode汉化 一、目前默认是英文状态 二、按下ShiftCtrlP弹出这个框框 三、点击语言显示配置 四、点击下载语言 五、点击下载中文即可 六、下载完后&#xff0c;再按ShiftCtrlP&#xff0c;选择语言配置 七、选择第二个&#xff0c;中…
阅读更多...
Swift Combine 级联多个 UI 更新，包括网络请求 从入门到精通十六

Swift Combine 级联多个 UI 更新，包括网络请求 从入门到精通十六

Combine 系列 Swift Combine 从入门到精通一Swift Combine 发布者订阅者操作者 从入门到精通二Swift Combine 管道 从入门到精通三Swift Combine 发布者publisher的生命周期 从入门到精通四Swift Combine 操作符operations和Subjects发布者的生命周期 从入门到精通五Swift Com…
阅读更多...
点餐|外卖订餐小程序|基于微信小程序的外卖订餐系统设计与实现(源码+数据库+文档)

点餐|外卖订餐小程序|基于微信小程序的外卖订餐系统设计与实现(源码+数据库+文档)

点餐|外卖订餐小程序目录 目录 基于微信小程序的外卖订餐系统设计与实现 一、前言 二、系统功能设计 三、系统实现 1、用户微信端功能模块 2、管理员服务端功能模块 3、商家务端功能模块 四、数据库设计 1、实体ER图 五、核心代码 六、论文参考 七、最新计算机毕设…
阅读更多...
Spring Boot整合 Cache 以Redis服务 处理数据缓存

Spring Boot整合 Cache 以Redis服务 处理数据缓存

目录 一、SpringBoot框架 二、什么是cache 三、redis介绍 四、高速缓存 一、SpringBoot框架 Spring Boot是一个基于Java的开源框架&#xff0c;用于快速构建独立的、可运行的、生产级的Spring应用程序。它简化了Spring应用程序的配置和部署过程&#xff0c;并提供了许多开…
阅读更多...
git stash 正确用法

git stash 正确用法

目录 一、背景 二、使用 2.1 使用之前&#xff0c;先简单了解下 git stash 干了什么&#xff1a; 2.2 git stash 相关命令 2.3 使用流程 1. 执行 git stash 2. 查看刚才保存的工作进度 git stash list 3. 这时候在看分支已经是干净无修改的(改动都有暂存到 stash) 4. 现在…
阅读更多...
山脉的个数/攀登者

山脉的个数/攀登者

题目描述 攀登者喜欢寻找各种地图&#xff0c;并且尝试攀登到最高的山峰。 地图表示为一维数组&#xff0c;数组的索引代表水平位置&#xff0c;数组的元素代表相对海拔高度。其中数组元素0代表地面。 例如&#xff1a;[0,1,2,4,3,1,0,0,1,2,3,1,2,1,0]&#xff0c;代表如下…
阅读更多...
最新文章

Copyright @ 2022~2023