mundo - programacion cuda linux
Reduce las filas de la matriz con CUDA (3)
Windows 7, NVidia GeForce 425M.
Escribí un código CUDA simple que calcula las sumas de fila de una matriz. La matriz tiene una representación unidimensional (puntero a un flotador).
La versión en serie del código está debajo (tiene 2 bucles, como se esperaba):
void serial_rowSum (float* m, float* output, int nrow, int ncol) {
float sum;
for (int i = 0 ; i < nrow ; i++) {
sum = 0;
for (int j = 0 ; j < ncol ; j++)
sum += m[i*ncol+j];
output[i] = sum;
}
}
Dentro del código CUDA, llamo a la función kernel barriendo la matriz por filas. A continuación, el fragmento de llamada del kernel:
dim3 threadsPerBlock((unsigned int) nThreadsPerBlock); // has to be multiple of 32
dim3 blocksPerGrid((unsigned int) ceil(nrow/(float) nThreadsPerBlock));
kernel_rowSum<<<blocksPerGrid, threadsPerBlock>>>(d_m, d_output, nrow, ncol);
y la función kernel que realiza la suma paralela de las filas (aún tiene 1 bucle):
__global__ void kernel_rowSum(float *m, float *s, int nrow, int ncol) {
int rowIdx = threadIdx.x + blockIdx.x * blockDim.x;
if (rowIdx < nrow) {
float sum=0;
for (int k = 0 ; k < ncol ; k++)
sum+=m[rowIdx*ncol+k];
s[rowIdx] = sum;
}
}
Hasta aquí todo bien. Los resultados seriales y paralelos (CUDA) son iguales.
El punto es que la versión de CUDA tarda casi el doble de tiempo que la de serie para calcular, incluso si cambio el parámetro nThreadsPerBlock : nThreadsPerBlock de 32 a 1024 (número máximo de hilos por bloque permitido para mi tarjeta).
IMO, la dimensión de la matriz es lo suficientemente grande como para justificar la paralelización: 90,000 x 1,000 .
A continuación, nThreadsPerBlock el tiempo transcurrido para las versiones en serie y paralelas utilizando diferentes nThreadsPerBlock . Tiempo informado en msec sobre un promedio de 100 muestras:
Matriz: nrow = 90000 x ncol = 1000
Serie: Tiempo promedio transcurrido por muestra en mseg ( 100 muestras): 289.18 .
CUDA ( 32 ThreadsPerBlock): Tiempo promedio transcurrido por muestra en mseg ( 100 muestras): 497.11 .
CUDA ( 1024 ThreadsPerBlock): tiempo promedio transcurrido por muestra en mseg ( 100 muestras): 699.66 .
Por si acaso, la versión con nThreadsPerBlock es la más rápida / más lenta.
Entiendo que hay una especie de sobrecarga al copiar de Host a Dispositivo y viceversa, pero tal vez la lentitud se debe a que no estoy implementando el código más rápido.
Ya que estoy lejos de ser un experto en CUDA:
¿Estoy codificando la versión más rápida para esta tarea? ¿Cómo podría mejorar mi código? ¿Puedo deshacerme del bucle en la función kernel?
Cualquier pensamiento apreciado.
EDIT 1
Aunque describo un rowSum estándar, estoy interesado en la operación AND / OR de filas que tienen valores (0;1} , como rowAND / rowOR . Dicho esto, no me permite explotar el multiplicador cuBLAS por 1 ''s Truco de vectores COL columna, según lo sugerido por algunos comentaristas.
EDIT 2
Como sugieren los usuarios otros usuarios y aquí endosados:
Olvídese de intentar escribir sus propias funciones , use la biblioteca Thrust en su lugar y la magia viene.
Si este es el alcance (sumando las filas) de las operaciones que necesita hacer con estos datos, no esperaría un beneficio considerable de la GPU. Tiene exactamente una operación aritmética por elemento de datos, y para eso está pagando el costo de transferir ese elemento de datos a la GPU. Y más allá de un determinado tamaño de problema (lo que sea necesario para mantener la máquina ocupada) no se obtiene un beneficio adicional de los tamaños de problema más grandes, porque la intensidad aritmética es O (n).
Este no es un problema particularmente emocionante para resolver en la GPU.
Pero como lo ha indicado talonmies, usted tiene un problema coalescente en la forma en que lo ha diseñado, lo que ralentizará aún más las cosas. Echemos un vistazo a un pequeño ejemplo:
C1 C2 C3 C4
R1 11 12 13 14
R2 21 22 23 24
R3 31 32 33 34
R4 41 42 43 44
Arriba hay un simple ejemplo pictórico de una pequeña porción de su matriz. El almacenamiento de datos de la máquina es tal que los elementos (11), (12), (13) y (14) se almacenan en ubicaciones de memoria adyacentes.
Para el acceso combinado , queremos un patrón de acceso tal que las ubicaciones de memoria adyacentes se soliciten a partir de la misma instrucción, ejecutada a lo largo de la urdimbre.
Necesitamos pensar en la ejecución de su código desde el punto de vista de un warp , es decir, 32 hilos ejecutándose en lock-step. ¿Qué está haciendo tu código? ¿Qué elementos está recuperando (pidiendo) en cada paso / instrucción? Echemos un vistazo a esta línea de código:
sum+=m[rowIdx*ncol+k];
Los hilos adyacentes en la urdimbre tienen valores adyacentes (es decir, consecutivos) para rowIdx ya que ha creado esa variable. Entonces, cuando k = 0, ¿qué elemento de datos solicita cada subproceso cuando intentamos recuperar el valor m[rowIdx*ncol+k] ?
En el bloque 0, el subproceso 0 tiene un rowIdx de 0. El subproceso 1 tiene un rowIdx de 1, etc. Por lo tanto, los valores solicitados por cada subproceso en esta instrucción son:
Thread: Memory Location: Matrix Element:
0 m[0] (11)
1 m[ncol] (21)
2 m[2*ncol] (31)
3 m[3*ncol] (41)
¡Pero esto no es un acceso combinado! Los elementos (11), (21), etc. no están adyacentes en la memoria. Para el acceso combinado, nos gustaría que la fila del elemento de la matriz se lea así:
Thread: Memory Location: Matrix Element:
0 m[?] (11)
1 m[?] (12)
2 m[?] (13)
3 m[?] (14)
Si luego trabajas hacia atrás para determinar el valor de ? debería ser, se te ocurrirá una instrucción como esta:
sum+=m[k*ncol+rowIdx];
Esto le dará acceso combinado, pero no le dará la respuesta correcta, porque ahora estamos sumando columnas de matriz en lugar de filas de matriz. Podemos solucionar esto reorganizando su almacenamiento de datos para que esté en orden de columna principal en lugar de orden de fila mayor. (Deberías ser capaz de buscar ideas en google, ¿no?) Conceptualmente, esto es equivalente a transponer tu matriz m . Si esto es conveniente para usted o no, está fuera del alcance de su pregunta, como yo lo veo, y no es realmente un problema de CUDA. Puede ser una tarea sencilla para usted, ya que está creando la matriz en el host o transfiriendo la matriz de un host a otro. Pero, en resumen, no conozco una forma de sumar las filas de la matriz con un 100% de acceso combinado, si la matriz está almacenada en orden mayor de fila. (Podría recurrir a una secuencia de reducciones de filas, pero eso me parece doloroso).
No es raro, cuando estamos pensando en formas de acelerar el código en la GPU, considerar la posibilidad de reorganizar nuestro almacenamiento de datos para facilitar la GPU. Este es un ejemplo.
Y, sí, lo que estoy esbozando aquí todavía conserva un bucle en el kernel.
Como comentario adicional, sugeriría que se sincronicen las porciones de copia de datos y las porciones de kernel (cálculo) por separado. No puedo decir a partir de su pregunta si está sincronizando solo el kernel o la operación completa (GPU), incluidas las copias de datos. Si ajusta las copias de datos por separado, puede descubrir que solo el tiempo de copia de datos excede su tiempo de CPU. Cualquier esfuerzo puesto en la optimización de su código CUDA no afectará el tiempo de copia de datos. Este podría ser un punto de datos útil antes de dedicar mucho tiempo a esto.
Como mencionó, necesita un algoritmo de reducción general que no sea solo suma. Trataré de dar 3 enfoques aquí. El enfoque del núcleo puede tener el mayor rendimiento. El enfoque de empuje es más fácil de implementar. El enfoque cuBLAS funciona solo con suma y tiene un buen rendimiento.
Enfoque Kernel
Aquí hay un documento muy bueno que presenta cómo optimizar la reducción paralela estándar. La reducción estándar se puede dividir en 2 etapas.
- Los múltiples bloques de hilos reducen una parte de los datos;
- Un bloque de hilo se reduce del resultado de la etapa 1 al elemento final 1.
Para su problema de reducción múltiple (reducir filas de tapetes), solo la etapa 1 es suficiente. La idea es reducir 1 fila por bloque de hilos. Para otras consideraciones, como multirruta por bloque de subprocesos o 1 fila por múltiples bloques de subprocesos, puede consultar el documento provisto por @Novak . Esto puede mejorar más el rendimiento, especialmente para matrices con mala forma.
Enfoque de empuje
La multirreducción general se puede realizar mediante thrust::reduction_by_key en unos minutos. Puede encontrar algunas discusiones aquí Determinando el elemento mínimo y su posición en cada columna de matriz con CUDA Thrust .
Sin embargo, thrust::reduction_by_key no supone que cada fila tenga la misma longitud, por lo que obtendrá una penalización de rendimiento. Otra publicación ¿Cómo normalizar columnas de matriz en CUDA con un rendimiento máximo? da comparación de perfiles entre thrust::reduction_by_key y cuBLAS enfoque en la suma de las filas. Puede darle una comprensión básica sobre el rendimiento.
enfoque CuBLAS
La suma de filas / columnas de una matriz A se puede ver como una multiplicación matriz-vector donde los elementos del vector son todos uno. se puede representar mediante el siguiente código de matlab.
y = A * ones(size(A,2),1);
donde y es la suma de las filas de A.
cuBLAS libary proporciona una función de multiplicación de matriz-vector de alto rendimiento cublas<t>gemv() para esta operación.
El resultado del tiempo muestra que esta rutina es solo 10 ~ 50% más lenta que simplemente leer todos los elementos de A una vez, lo que puede verse como el límite superior teórico del rendimiento para esta operación.
La reducción de las filas de una matriz se puede resolver mediante el uso de CUDA Thrust de tres maneras (pueden no ser las únicas, pero abordar este punto está fuera del alcance). Como también es reconocido por el mismo OP, usar CUDA Thrust es preferible para este tipo de problema. Además, es posible un enfoque usando cuBLAS .
ENFOQUE # 1 - reduce_by_key
Este es el enfoque sugerido en esta página de ejemplo de Thrust . Incluye una variante usando make_discard_iterator .
ENFOQUE # 2 - transform
Este es el enfoque sugerido por Robert Crovella en CUDA Thrust: reduce_by_key en solo algunos valores en una matriz, basado en valores en una matriz "clave" .
ENFOQUE # 3 - inclusive_scan_by_key
Este es el enfoque sugerido por Eric en Cómo normalizar columnas de matriz en CUDA con un rendimiento máximo. .
ENFOQUE # 4 - cublas<t>gemv
Utiliza cuBLAS gemv para multiplicar la matriz relevante por una columna de 1 ''s.
EL CÓDIGO COMPLETO
Aquí está el código que condensa los dos enfoques. Los archivos Utilities.cu y Utilities.cuh se mantienen aquí y se omiten aquí. El TimingGPU.cu y TimingGPU.cuh se mantienen aquí y también se omiten.
#include <cublas_v2.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/reduce.h>
#include <thrust/functional.h>
#include <thrust/random.h>
#include <thrust/sequence.h>
#include <stdio.h>
#include <iostream>
#include "Utilities.cuh"
#include "TimingGPU.cuh"
// --- Required for approach #2
__device__ float *vals;
/**************************************************************/
/* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */
/**************************************************************/
template <typename T>
struct linear_index_to_row_index : public thrust::unary_function<T,T> {
T Ncols; // --- Number of columns
__host__ __device__ linear_index_to_row_index(T Ncols) : Ncols(Ncols) {}
__host__ __device__ T operator()(T i) { return i / Ncols; }
};
/******************************************/
/* ROW_REDUCTION - NEEDED FOR APPROACH #2 */
/******************************************/
struct row_reduction {
const int Ncols; // --- Number of columns
row_reduction(int _Ncols) : Ncols(_Ncols) {}
__device__ float operator()(float& x, int& y ) {
float temp = 0.f;
for (int i = 0; i<Ncols; i++)
temp += vals[i + (y*Ncols)];
return temp;
}
};
/**************************/
/* NEEDED FOR APPROACH #3 */
/**************************/
template<typename T>
struct MulC: public thrust::unary_function<T, T>
{
T C;
__host__ __device__ MulC(T c) : C(c) { }
__host__ __device__ T operator()(T x) { return x * C; }
};
/********/
/* MAIN */
/********/
int main()
{
const int Nrows = 5; // --- Number of rows
const int Ncols = 8; // --- Number of columns
// --- Random uniform integer distribution between 10 and 99
thrust::default_random_engine rng;
thrust::uniform_int_distribution<int> dist(10, 99);
// --- Matrix allocation and initialization
thrust::device_vector<float> d_matrix(Nrows * Ncols);
for (size_t i = 0; i < d_matrix.size(); i++) d_matrix[i] = (float)dist(rng);
TimingGPU timerGPU;
/***************/
/* APPROACH #1 */
/***************/
timerGPU.StartCounter();
// --- Allocate space for row sums and indices
thrust::device_vector<float> d_row_sums(Nrows);
thrust::device_vector<int> d_row_indices(Nrows);
// --- Compute row sums by summing values with equal row indices
//thrust::reduce_by_key(thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_row_index<int>(Ncols)),
// thrust::make_transform_iterator(thrust::counting_iterator<int>(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
// d_matrix.begin(),
// d_row_indices.begin(),
// d_row_sums.begin(),
// thrust::equal_to<int>(),
// thrust::plus<float>());
thrust::reduce_by_key(
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)),
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
d_matrix.begin(),
thrust::make_discard_iterator(),
d_row_sums.begin());
printf("Timing for approach #1 = %f/n", timerGPU.GetCounter());
// --- Print result
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums[i] << "/n";
}
/***************/
/* APPROACH #2 */
/***************/
timerGPU.StartCounter();
thrust::device_vector<float> d_row_sums_2(Nrows, 0);
float *s_vals = thrust::raw_pointer_cast(&d_matrix[0]);
gpuErrchk(cudaMemcpyToSymbol(vals, &s_vals, sizeof(float *)));
thrust::transform(d_row_sums_2.begin(), d_row_sums_2.end(), thrust::counting_iterator<int>(0), d_row_sums_2.begin(), row_reduction(Ncols));
printf("Timing for approach #2 = %f/n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_2[i] << "/n";
}
/***************/
/* APPROACH #3 */
/***************/
timerGPU.StartCounter();
thrust::device_vector<float> d_row_sums_3(Nrows, 0);
thrust::device_vector<float> d_temp(Nrows * Ncols);
thrust::inclusive_scan_by_key(
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)),
thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols)) + (Nrows*Ncols),
d_matrix.begin(),
d_temp.begin());
thrust::copy(
thrust::make_permutation_iterator(
d_temp.begin() + Ncols - 1,
thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(Ncols))),
thrust::make_permutation_iterator(
d_temp.begin() + Ncols - 1,
thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC<int>(Ncols))) + Nrows,
d_row_sums_3.begin());
printf("Timing for approach #3 = %f/n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_3[i] << "/n";
}
/***************/
/* APPROACH #4 */
/***************/
cublasHandle_t handle;
timerGPU.StartCounter();
cublasSafeCall(cublasCreate(&handle));
thrust::device_vector<float> d_row_sums_4(Nrows);
thrust::device_vector<float> d_ones(Ncols, 1.f);
float alpha = 1.f;
float beta = 0.f;
cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_T, Ncols, Nrows, &alpha, thrust::raw_pointer_cast(d_matrix.data()), Ncols,
thrust::raw_pointer_cast(d_ones.data()), 1, &beta, thrust::raw_pointer_cast(d_row_sums_4.data()), 1));
printf("Timing for approach #4 = %f/n", timerGPU.GetCounter());
for(int i = 0; i < Nrows; i++) {
std::cout << "[ ";
for(int j = 0; j < Ncols; j++)
std::cout << d_matrix[i * Ncols + j] << " ";
std::cout << "] = " << d_row_sums_4[i] << "/n";
}
return 0;
}
RESULTADOS DE TIEMPO (probado en un Kepler K20c)
Matrix size #1 #1-v2 #2 #3 #4 #4 (no plan)
100 x 100 0.63 1.00 0.10 0.18 139.4 0.098
1000 x 1000 1.25 1.12 3.25 1.04 101.3 0.12
5000 x 5000 8.38 15.3 16.05 13.8 111.3 1.14
100 x 5000 1.25 1.52 2.92 1.75 101.2 0.40
5000 x 100 1.35 1.99 0.37 1.74 139.2 0.14
Parece que los enfoques n. ° 1 y n. ° 3 superan al enfoque n. ° 2, excepto en el caso de pequeñas cantidades de columnas. El mejor enfoque, sin embargo, es el enfoque n.º 4, que es significativamente más conveniente que los demás, siempre que el tiempo necesario para crear el plan pueda amortizarse durante el cálculo.