Часть кода, над которым я работаю, требует как можно более быстрого выполнения векторного умножения матриц, то есть с использованием оптимизированной сторонней библиотеки, такой как cublas (хотя тот же принцип применим к любому процессору).
Проблема в том, что между элементами в векторе есть своего рода шаг:
Матрица хранится в виде трехмерного массива с плавающей запятой 3Nx3N.
Вектор сохраняется как массив N 1D float4s, но должны использоваться только первые три элемента каждого float4, четвертый должен игнорироваться.
N составляет порядка миллионов.
Если бы вектор хранился как float3 вместо float4, я мог бы просто привести указатель к float, как в этом рабочем примере:
//Compile with nvcc test.cu -O3 -lcublas -o test
/*
Multiply a 3Nx3N float matrix, M, by a vector, X, of N float3 elements
The result, Y, is a 3N float vector
-----------------------
What if X is a vector of N float4?
How can I tell cublas to skip the forth element?
*/
#include<iostream>
#include<thrust/device_vector.h>
#include<cuda_runtime.h>
#include<cublas_v2.h>
using namespace std;
int main(){
int N = 3;
thrust::device_vector<float3> X(N);
thrust::device_vector<float> Y(3*N);
for(int i=0; i<N; i++)
X[i] = make_float3(1,1,1); //make_float4(1,1,1,0); //in the case of float4 i.e., The result should be the same
thrust::device_vector<float> M(3*N*3*N, 1);cublasHandle_t handle;
cublasCreate(&handle);
float beta = 0.0f;
float alpha = 1.0f;
cublasSgemv(handle, CUBLAS_OP_T,
3*N, 3*N,
&alpha,
thrust::raw_pointer_cast(&M[0]), 3*N,
(float*) thrust::raw_pointer_cast(&X[0]), 1,
&beta,
thrust::raw_pointer_cast(&Y[0]), 1);
cout<<"Performed Y = M·X\n\tX = ";
for(int i=0; i<N; i++){
float3 Xi = X[i];
cout<<Xi.x<<" "<<Xi.y<<" "<<Xi.z<<" ";
}
cout<<"\n\tY = ";
for(int i=0; i<3*N; i++){
cout<<Y[i]<<" ";
}
cout<<endl;
return 0;
}
Но как я могу выполнить эту операцию, если вектор X хранится как float4 s?
Учитывая, что float4 * можно интерпретировать как float * с в 4 раза большим количеством элементов, вопрос может быть более общим (хотя меня интересует только случай float4);
Если между 3 «полезными» элементами есть шаг. Я хочу сказать кублам, что массив не слипается в памяти. Но что-то вроде: есть 3 элемента в начале, следующие три — это элементы «шага» после этого и т. Д.
Подобно тому, что вы можете сделать в OpenGL с объектами массива вершин.
РЕДАКТИРОВАТЬ:
Ответы показали, что наиболее жизнеспособный метод — просто скопировать выделенный массив во временный преобразованный массив float3, понятный cublas..
На данный момент есть два варианта:
1. Use cudaMemcpy2D
2. Use a thrust transformation
3. Use a custom copy kernel
Я написал этот код для проверки трех случаев:
//Compile with Compile with: nvcc test.cu -O3 -lcublas -o test
#include<iostream>
#include<thrust/device_vector.h>
#include<cuda.h>
#include<cuda_runtime.h>
#include<cublas_v2.h>
using namespace std;struct Timer{
cudaEvent_t start, stop;
float time;
void tic(){
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
}
float toc(){
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&time, start, stop);
cudaEventDestroy(start);
cudaEventDestroy(stop);
return time;
}
};struct copy_functor{
copy_functor(){}
__device__ float3 operator() (const float4& X4){
return make_float3(X4.x, X4.y, X4.z);
}
};__global__ void copy_kernel(const float4* __restrict__ X4, float3* __restrict__ X3, int N){
int id = blockIdx.x*blockDim.x + threadIdx.x;
if(id < N){
float4 x4 = X4[id];
X3[id] = make_float3(x4.x, x4.y, x4.z);
}
}
int main(){
int N = 1000000;
int Ntest = 1000;
Timer t;
thrust::device_vector<float3> X3(N, make_float3(0,0,0));
thrust::device_vector<float4> X4(N, make_float4(1,1,1,10));/*************************CUDAMEMCPY2D*******************/
t.tic();
for(int i= 0; i<Ntest; i++){
cudaMemcpy2DAsync(thrust::raw_pointer_cast(&X3[0]),
3*sizeof(float),
thrust::raw_pointer_cast(&X4[0]),
4*sizeof(float),
3*sizeof(float),
N,
cudaMemcpyDeviceToDevice);
cudaDeviceSynchronize();
}
printf ("Time for cudaMemcpy2DAsync: %f ms\n", t.toc()/(float)Ntest);/************************THRUST***********************/
t.tic();
for(int i= 0; i<Ntest; i++){
transform(X4.begin(), X4.end(), X3.begin(), copy_functor());
cudaDeviceSynchronize();
}
printf ("Time for thrust transformation: %f ms\n", t.toc()/(float)Ntest);
/*********************COPY KERNEL*****************************/
t.tic();
for(int i= 0; i<Ntest; i++){
copy_kernel<<< N/128 + 1, 128 >>>(thrust::raw_pointer_cast(&X4[0]),
thrust::raw_pointer_cast(&X3[0]), N);
cudaDeviceSynchronize();
}
printf ("Time for copy kernel: %f ms\n", t.toc()/(float)Ntest);return 0;
}
Обратите внимание, что я выполняю в среднем 1000 копий.
Вывод этого кода в GTX 980 следующий:
Time for cudaMemcpy2DAsync: 1.465522 ms
Time for thrust transformation: 0.178745 ms
Time for copy kernel: 0.168507 ms
cudaMemcpy2D на порядок медленнее остальных.
тяги и копирование ядра очень похожи и самый быстрый способ
Такое поведение остается с любым количеством элементов.
EDIT2:
Другие ответы предполагают, что GEMM может быть использован для передачи шага. Без необходимости временного массива.
Интерпретация матричного вектора мул. как матрица матрица мул. будет сделано так:
cublasSgemm(handle, CUBLAS_OP_T, CUBLAS_OP_T,
3*N, 1 /*m*/, 3*N,
&alpha,
thrust::raw_pointer_cast(&M[0]), 3*N,
(float*) thrust::raw_pointer_cast(&X3[0]), 1 /*ldb*/,
&beta,
thrust::raw_pointer_cast(&Y[0]), 3*N);
Однако, на данный момент, я не знаю, как пройти X4 вместо X3. Похоже, что решение находится в параметрах m и ldb.
Вы можете рассматривать ваш 1-D вектор float4 как матрицу Nx3 2-D с плавающей строкой 4 и использовать cudaMemcpy2DAsync
изменить шаг с 4 до 3 с
cudaMemcpy2DAsync(dst,
3*sizeof(float),
src,
4*sizeof(float),
3*sizeof(float),
N,
cudaMemcpyDeviceToDevice);
Тогда dst
может рассматриваться как трехмерный вектор с плавающей запятой и передаваться в gemv()
непосредственно.
Учитывая масштаб вашего N
время копирования не заметно по сравнению с gemv()
,
РЕДАКТИРОВАТЬ
Результат теста @Apo показывает, что лучше использовать ядро копии вместо cudaMemcpy2DAsync
, Я был слишком ожидаем cudaMemcpy2DAsync
и думал, что это будет хорошо оптимизировано и будет иметь лучшую производительность для всех случаев.
Других решений пока нет …