【CUDA并行编程之五】计算向量的欧式距离

时间：2014-12-25 13:09:42 阅读：259 评论：0 收藏：0 [点我收藏+]

本文将介绍如何用cuda来计算两个向量之间的欧式距离，其中涉及到了如果将二维矩阵传入到核函数进行计算的问题，并且介绍两个内存分配和拷贝的API：cudaMallocPitch以及cudaMemcpy2D。

一、需求分析

现在我们要解决这么一个问题：计算一个D维的向量A[D]到二维矩阵B[N][D]的每一行的欧式距离，并且将每一组距离保存在一个向量dis[N]中并返回。我们还是通过串行和并行两种方式来进行实现。

二、串行实现

实现方法就是用一个二重循环进行相乘，然后将结果保存。上代码：

dis_cal_sequence.cc:

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

using namespace std;

const int N = 100;
const int D = 8;
const int MAX = 10;

void cal_dis(int **train_data, int *test_data, int *dis)
{
	for(int k=0;k<N;k++)
	{
		int sum = 0;
		int temp = 0;
		for(int i=0;i<D;i++)
		{
			temp =  *(*(train_data+k)+i) - test_data[i];
			sum += temp * temp;
		}
		dis[k] = sum;
	}
}

void print(int **data)
{
	cout<<"training data:"<<endl;
 	for(int i=0;i<N;i++)
	{
		for(int j=0;j<D;j++)
		{
			cout<<*(*(data+i)+j)<<" ";		
		}
		cout<<endl;
	}
}

void print(int *data,int n)
{
	for(int i=0;i<n;i++)
	{
		cout<<data[i]<<" ";
	}
	cout<<endl;
}

int main()
{
	int **h_train_data , *h_test_data , *distance;
	
	//allocate space in heap for the variable
	h_train_data = new int*[N];
	for(int i=0;i<N;i++)
	{
		h_train_data[i] = new int[D];
	}
	h_test_data = new int[D];
	distance = new int[N];

	//initialize training data
	srand( (unsigned)time(NULL) );
	for( int i=0;i<N;i++ )
	{
		for( int j=0;j<D;j++)
		{
			h_train_data[i][j] = rand()%MAX;
		}
	}
	print(h_train_data);

	//initialize testing data
	for( int j=0;j<D;j++ )
	{
	  	h_test_data[j] = rand() % MAX;
	}
	cout<<"testing data:"<<endl;
	print(h_test_data,D);

	//calculate the distance
	cal_dis( h_train_data,h_test_data,distance );

	cout<<"distance data:"<<endl;
	print(distance,N);

	return 0;
}

结果：

代码没有什么太多值得将的，主要讲一下C/C++中二维矩阵的传参问题。还是做一个总结吧。

对于静态数组：

void fun( int array[][D])// OK
void fun( int (*array)[D])//OK
void fun( int **array) //error

int main()
{
	int A[N][D];
	fun(A);
}

对于动态数组：

void fun( int array[][D])// Error
void fun( int (*array)[D])//Error
void fun( int **array) //OK

int main()
{
	int **A;
	A = new int*[N];
	for(int i=0;i<N;i++)
		A[i] = new int[D];
	fun(A);
}

当然在函数内部两种访问数组的方式都可以：array[i][j]和*( *(array+i) + j)。

空间复杂度O(N*D)，时间复杂度O(N*D)。

三、向CUDA内核中传递二位向量

在介绍代码之前首先将一个如何将一个二维向量传入到核函数中进行计算的问题，被这个问题困扰很久。我们一般都是将二维矩阵转换成为一维向量进行计算，那么如果我偏偏想用二维向量计算该怎么办呢？（PS：在此吐个槽，国内网站上的技术问答有时候很难找到合适的答案，搜来搜去浪费了很多时间，但是在stackoverflow上搜很快就搜到了~）

在CUDA上分配二维数组还是有点让人困惑的，可能会犯如下两种错误：

错误一：

int rowCount = 10;  
float** d_array=0; // array on device  
cudaMalloc(d_array, rowCount*sizeof(float*));  
for(int i = 0 ; i < rowCount ; i++)  
{  
	//this results in error "Access violation writing location"  
	cudaMalloc((void **)&d_array[i], (i + 1) * sizeof(float)); // column length increases with "i" here  
}

这里有点类似于二维向量的动态分配过程。这样做的问题是，cudaMalloc在内核里面一次性的获取到内存空间，一旦主线程分配完一维空间，那么就会造成当分配第二维空间的时候会抛出"Access violation writing location"的异常。

错误二：

int rowCount = 10;  
float** d_array =(float**)malloc(rowCount*sizeof(float*)); //malloc 1st dimension  
for(int i = 0 ; i < rowCount ; i++)  
{  
cudaMalloc(&d_array[i], (i + 1) * sizeof(float)); // cuda malloc 2nd dimension  
}

这个问题在于一维数组是在主机内存上分配的（使用了malloc），而第二维空间是在设备上分配的（使用了cudaMalloc）。这会导致内存获取失败。

那么如何来解决这个问题？我们使用两个API：cudaMallocPitch&cudaMemcpy2D。

3.1：cudaMallocPitch()

cudaMallocPitch( void**devPtr, size_t* pitch ,size_t widthInBytes , size_t height)

devPtr：指向矩阵的内存空间的头指针。

pitch：分配存储器的宽度，以字节为单位

width：分配矩阵的列数

height：分配矩阵的行数

在设备上分配widthInBytes * height字节的线性内存，并返回分配内存的指针*devPtr。函数将确保在任何给出的行中对应的指针是连续的。pitch返回的指针*pitch是分配的宽度。Pitch作为内存分配的一个分开的参数，用来计算2D数组中的地址。一个给定行和列的类型T的数组元素，地址等于：

T* pElement = (T*)((char*)BaseAddress + Row * pitch) + Column;

对于2D数组的分配，建议使用cudaMallocPitch()分配内存。由于pitch对列限制受限于硬件，特别是当应用程序从设备内存的不同区域执行一个2D的内存拷贝。

例如：要在GPU上开辟一个行数为height，列数为width的float型矩阵空间。

int width = 64, height = 64;
float* devPtr;
size_t pitch;
cudaMallocPitch((void**)&devPtr,&pitch, width * sizeof(float), height)

假设GPU中global memory被划分为128 Byte的段（0-127，128-255, 256-383,……），你要为内存分配float型矩阵数据。假设矩阵的行数为N，列数为M，在C和C++中数据是按行优先存储的。当N*4是128倍数的时候(float型数据占4个Byte)，那么用cudaMalloc分配出来的内存空间也是对齐的，也就是说你一行分配的字节数刚好是128 Bytes的倍数的时候，cudaMalloc()也是对齐的；

再看另一种情况，N=33，此时33*4=132，不是128的倍数，当warp从第二行开头（从132开始读，第一行是0到131）读取global memory的时候，首地址并非global memory划分的对齐段的首地址，那么这样的访问就是非合并的，cudaMallocPitch()就是为了解决每行首地址是否是global memory对齐段的问题，如果用cudaMallocPitch()来分配N=33(即列数为33)的矩阵时，每一行大小会变成256个Bytes(0-131为我们需要使用的空间，132-255未使用)，而不是cudaMalloc中的132个Bytes，这样分配以后，每行的首地址将会是与globla memory分段地址对齐的（都是128的整数倍），warp在访问的时候就可以对齐了！

3.2：cudaMemcpy2D

cudaMemcpy2D(void* dst,size_t dpitch,const void* src,size_t spitch,size_t width,size_t height,enum cudaMemcpyKind kind);
dst: 目的矩阵内存头指针
dpitch: dst指向的2D数组中的内存宽度，以字节为单位，是cuda为了读取方便，
对齐过的内存宽度，可能大于一行元素占据的实际内存。
src：源矩阵内存头指针
spitch: src指向的2D数组中的内存宽度，以字节为单位
width: src指向的2D数组中一行元素占据的实际宽度。以字节为单位，等于width*sizeof(type)
height: src指向的2D数组的行数
kind：拷贝数据的方向

从src指向的内存区域拷贝数据到dst指向的内存区域。kind表示拷贝方向：cudaMemcpyHostToHost, cudaMemcpyDeviceToHost, cudaMemcpyHostToDevice,或者cudaMemcpyDeviceToDevice。

四、并行实现

方法分析：并行的情况下让GPU的每一个线程去处理两个一维向量之间的欧式距离，那么当N个线程并行的情况下就能将所有结果计算出来，这样的时间复杂度会降到O(D)。

上代码，dis_cal_parallel.cu：

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

using namespace std;

const int N = 10;
const int D = 8;
const int MAX = 10;

__global__ void cal_dis(int *train_data, int *test_data, int *dis,int pitch)
{
	int tid = blockIdx.x;
	if(tid<N)
	{
		int temp = 0;
		int sum = 0;
		for(int i=0;i<D;i++)
		{
			temp = *((int*)((char*)train_data + tid * pitch) + i) - test_data[i];
			sum += temp * temp;
		}
		dis[tid] = sum;
	}
}

void print(int data[][D])
{
	cout<<"training data:"<<endl;
 	for(int i=0;i<N;i++)
	{
		for(int j=0;j<D;j++)
		{
			cout<<*(*(data+i)+j)<<" ";		
		}
		cout<<endl;
	}
}

void print(int *data,int n)
{
	for(int i=0;i<n;i++)
	{
		cout<<data[i]<<" ";
	}
	cout<<endl;
}

int main()
{
	int h_train_data[N][D] , h_test_data[D] , distance[N];

	int *d_train_data , *d_test_data , *d_dis;

	size_t pitch_d;
	size_t pitch_h = D * sizeof(int) ;

	//allocate memory on GPU 
	cudaMallocPitch( &d_train_data , &pitch_d , D * sizeof(float) , N ); 
	cudaMalloc( (void**)&d_test_data ,  D*sizeof(int) );
	cudaMalloc( (void**)&d_dis , 		  N*sizeof(int) );

	//initialize training data
	srand( (unsigned)time(NULL) );
	for( int i=0;i<N;i++ )
	{
		for( int j=0;j<D;j++)
		{
			h_train_data[i][j] = rand()%MAX;
		}
	}
	print(h_train_data);

	//initialize testing data
	for( int j=0;j<D;j++ )
	{
	  	h_test_data[j] = rand() % MAX;
	}
	cout<<"testing data:"<<endl;
	print(h_test_data,D);

	//copy training and testing data from host to device
	cudaMemcpy2D( d_train_data , pitch_d , h_train_data , pitch_h , D * sizeof(int) , N , cudaMemcpyHostToDevice );
	cudaMemcpy( d_test_data,  h_test_data ,  D*sizeof(int), cudaMemcpyHostToDevice);

	//calculate the distance
	cal_dis<<<N,1>>>( d_train_data,d_test_data,d_dis,pitch_d );

	//copy distance data from device to host
	cudaMemcpy( distance , d_dis  , N*sizeof(int) , cudaMemcpyDeviceToHost);

	cout<<"distance:"<<endl;;
	print(distance , N);

	cudaFree(d_train_data);
	cudaFree(d_test_data);
	cudaFree(d_dis);

	return 0;
}

结果：