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种情况讨论：

如果最短路经过点 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.
若最短路径不经过点 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)

在实际算法中，为了节约空间，可以直接在原来空间上进行迭代，这样空间可降至二维。

分析

时间复杂度： O ( V 3 ) O(V^3) O(V3)，其中 V V V是点集，对于 i , j i,j i,j两层for循环可使用OpenMP优化到线性
空间复杂度： 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个阶段为：

首先将 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。
对节点 ( K , K ) (K, K) (K,K)所在的第 K K K列进行MMA即矩阵乘法加法操作
对节点 ( K , K ) (K, K) (K,K)所在的第 K K K行进行MMA即矩阵乘法加法操作
对于除以上涉及到的剩余节点

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∑bxi,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 #include 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 #include #include #include 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<<>>(d_dist, r);
		phase3<<>>(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 #include #include #include #include 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<<>>(d_dist[device_id], r);
			phase3<<>>(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

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

分类:热门推荐日期:2024-03-11浏览:1评论:0