CPU-GPU异构并行化APSP算法

news2024/11/25 18:45:10

一、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

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

[ai笔记7] google浏览器ai学习提效定制优化+常用插件推荐

[ai笔记7] google浏览器ai学习提效定制优化+常用插件推荐

欢迎来到文思源想的ai空间&#xff0c;这是技术老兵重学ai以及成长思考的第7篇分享&#xff01; 工欲善其事必先利其器&#xff0c;为了ai学习的效能提升&#xff0c;放假期间对google浏览器做了一次系统整改&#xff0c;添加了一些配置和插件&#xff0c;这里既有一些显示、主…
阅读更多...
opencv通道分离与合并

opencv通道分离与合并

void QuickDemo::channels_demo(Mat & image) {std::vector<Mat>mv;//通道分离合并split(image,mv);//原图 指针(Mat)imshow("蓝色", mv[0]);imshow("绿色", mv[1]);imshow("红色", mv[2]); } split(image,mv);//原图 指针(Mat) 这里…
阅读更多...
AI：128-基于机器学习的建筑物能源消耗预测

AI：128-基于机器学习的建筑物能源消耗预测

🚀点击这里跳转到本专栏，可查阅专栏顶置最新的指南宝典~ 🎉🎊🎉 你的技术旅程将在这里启航！ 从基础到实践，深入学习。无论你是初学者还是经验丰富的老手，对于本专栏案例和项目实践都有参考学习意义。 ✨✨✨ 每一个案例都附带有在本地跑过的关键代码，详细讲解供…
阅读更多...
OpenAI Sora 初体验

OpenAI Sora 初体验

OpenAI Sora 初体验 就在刚刚&#xff0c;OpenAI 再次投下一枚重磅炸弹——Sora&#xff0c;一个文本到视频生成模型。 我第一时间体验了 Sora。看过 Sora 的能力后&#xff0c;我真的印象深刻。对细节的关注、无缝的角色刻画以及生成视频的绝对质量真正将可能性提升到了一个新…
阅读更多...
C++数据结构与算法——字符串

C++数据结构与算法——字符串

C第二阶段——数据结构和算法&#xff0c;之前学过一点点数据结构&#xff0c;当时是基于Python来学习的&#xff0c;现在基于C查漏补缺&#xff0c;尤其是树的部分。这一部分计划一个月&#xff0c;主要利用代码随想录来学习&#xff0c;刷题使用力扣网站&#xff0c;不定时更…
阅读更多...
STM32 寄存器操作 GPIO 与下降沿中断

STM32 寄存器操作 GPIO 与下降沿中断

一、如何使用stm32寄存器点灯&#xff1f; 1.1 寄存器映射表 寄存器本质就是一个开关&#xff0c;当我们把芯片寄存器配置指定的状态时即可使用芯片的硬件能力。 寄存器映射表则是开关的地址说明。对于我们希望点亮 GPIO_B 的一个灯来说&#xff0c;需要关注以下的两个寄存器…
阅读更多...
PLC_博图系列☞FBD

PLC_博图系列☞FBD

PLC_博图系列☞FBD 文章目录 PLC_博图系列☞FBD背景介绍FBD优势局限性 FBD 元素 关键字&#xff1a; PLC、 西门子、 博图、 Siemens 、 FBD 背景介绍 这是一篇关于PLC编程的文章&#xff0c;特别是关于西门子的博图软件。我并不是专业的PLC编程人员&#xff0c;也不懂电路…
阅读更多...
1036 跟奥巴马一起编程 (15)

1036 跟奥巴马一起编程 (15)

美国总统奥巴马不仅呼吁所有人都学习编程&#xff0c;甚至以身作则编写代码&#xff0c;成为美国历史上首位编写计算机代码的总统。2014 年底&#xff0c;为庆祝“计算机科学教育周”正式启动&#xff0c;奥巴马编写了很简单的计算机代码&#xff1a;在屏幕上画一个正方形。现在…
阅读更多...
【JavaEE】_文件与IO

【JavaEE】_文件与IO

目录 1.文件概述 1.1 文件的概念 1.2 文件的存储 1.3 文件的分类 1.4 目录结构 1.5 文件操作 1.5.1 文件系统操作 1.5.2 文件内容操作 2. Java文件系统操作 2.1 File类所处的包 2.2 构造方法 2.3 方法 2.3.1 与文件路径、文件名有关的方法 2.3.2 文件是否存在与普…
阅读更多...
c++之function和bind详解-SurfaceFlinger学习必备语法基础

c++之function和bind详解-SurfaceFlinger学习必备语法基础

背景 C中的function和bind是为了更方便地进行函数对象的封装和调用而设计的&#xff0c;在SurfaceFlinger源码中也是有很多使用部分。 比如分析Vsync相关源码时候有相关回调时候 可以看到这里的mRegistration就有个参数是 std::bind,怎么这里就可以进行回调呢&#xff1f; 所…
阅读更多...
代码随想录 Leetcode860. 柠檬水找零

代码随想录 Leetcode860. 柠檬水找零

题目&#xff1a; 代码(首刷自解 2024年2月15日&#xff09;&#xff1a; class Solution { public:bool lemonadeChange(vector<int>& bills) {vector<int> leftchange(2,0);//leftchange[0]代表5元数量&#xff0c;1代表10for (int i 0; i < bills.size…
阅读更多...
最小生成树（Kruskal算法及相关例题）

最小生成树（Kruskal算法及相关例题）

1.Kruskal算法概念以及基本思路 &#xff08;1&#xff09;概念&#xff1a; 克鲁斯卡尔算法是求连通网的最小生成树的另一种方法。它的时间复杂度为O&#xff08;ElogE&#xff09;(E是图G的边的总数)&#xff0c;适合于求边稀疏的网的最小生成树 。 其基本思想是&#xff…
阅读更多...
OS文件管理

OS文件管理

文件管理 文件的属性 文件所包含的属性&#xff1a; 文件名&#xff1a;由创建文件的用户决定文件名&#xff0c;主要为了方便用户找到文件&#xff0c;同一目录下不允许有重名文件。标识符&#xff1a;一个系统内的各文件标识符唯一&#xff0c;对用户来说毫无可读性&#…
阅读更多...
2.12：C语言测试题

2.12：C语言测试题

1.段错误&#xff1a;申请堆区内存未返回&#xff0c;str指向NULL 2.段错误&#xff1a;局部变量&#xff0c;本函数结束&#xff0c;p也释放 3.越界访问&#xff0c;可能正常输出hello&#xff0c;可能报错 4.可能段错误&#xff0c;释放后&#xff0c;str未指向NULL&#x…
阅读更多...
CentOS7.9+Kubernetes1.29.2+Docker25.0.3高可用集群二进制部署

CentOS7.9+Kubernetes1.29.2+Docker25.0.3高可用集群二进制部署

CentOS7.9Kubernetes1.29.2Docker25.0.3高可用集群二进制部署 Kubernetes高可用集群&#xff08;Kubernetes1.29.2Docker25.0.3&#xff09;二进制部署二进制软件部署flannel v0.22.3网络&#xff0c;使用的etcd是版本3&#xff0c;与之前使用版本2不同。查看官方文档进行了解…
阅读更多...
The method toList() is undefined for the type Stream

The method toList() is undefined for the type Stream

The method toList() is undefined for the type Stream &#xff08;JDK16&#xff09; default List<T> toList() { return (List<T>) Collections.unmodifiableList(new ArrayList<>(Arrays.asList(this.toArray()))); }
阅读更多...
C语言strstr函数

C语言strstr函数

简介 strstr(str1,str2) 函数用于判断字符串str2是否是str1的子串。如果是&#xff0c;则该函数返回 str1字符串从 str2第一次出现的位置开始到 str1结尾的字符串&#xff1b;否则&#xff0c;返回NULL。 实验 #include "stdio.h" #include "string.h"c…
阅读更多...
2.14：二维数组、非函数实现strcat、strcmp、strcpy、strlen

2.14：二维数组、非函数实现strcat、strcmp、strcpy、strlen

1.编程实现二维数组的杨辉三角 程序代码&#xff1a; 1 #include<stdio.h>2 #include<string.h>3 #include<stdlib.h>4 int main(int argc, const char *argv[])5 {6 int n;7 printf("please enter n:");8 scanf("%d",&…
阅读更多...
2024/02/13

2024/02/13

21 、C 22 、D 23、B 如果5先出栈那么1&#xff0c;2&#xff0c;3&#xff0c;4就已经入栈了&#xff0c;5出后4出&#xff0c;1要出栈必须先让3&#xff0c;2出栈&#xff0c;所以 不可能输出B 24、10&#xff0c;12&#xff0c;120 25、2&#xff0c;5 26、段错…
阅读更多...
2024.02.12作业

2024.02.12作业

1. 段错误 2. 段错误 3. hello 4. world 5. int a; int* a; int **a; int a[10]; int* a[10]; int(* a)[10]; int* a(int); int (*a[10])(int); 6. 6&#xff1b; 2&#xff1b; 2 7. 2 8. 2 9. b 10. a 11. a 12. c 13. b 14. c 15. a 16. c 17. b 18. a 19…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023