cuda - Impacto de la dispersión de matrices en cblas sgemm en Ubuntu 14.04

ubuntu-14.04 cublas

Question

Recientemente descubrí que el rendimiento de una llamada a cblas_sgemm para la multiplicación de matrices mejora drásticamente si las matrices tienen un "gran" número de ceros en ellas. Mejora hasta el punto en que supera a su primo cubila alrededor de 100 veces. Esto podría atribuirse muy probablemente a alguna detección automática de dispersión y conversión de formato adecuada mediante la función cblas_sgemm.

Desafortunadamente, tal comportamiento no se exhibe por su contraparte de cuda, es decir, cublasSgemm.

Entonces, la pregunta es, ¿cómo puedo obtener el mismo tipo de optimización en cublasSgemm para matrices que pueden tener una gran cantidad de ceros.

¿Y qué técnica usa cblas_sgemm para ajustarse automáticamente a matrices dispersas?

Por favor, no recomiende cuSparse / CUSP, etc. porque

1. No estoy seguro acerca de la escasez de matrices de entrada de antemano
2. Estoy trabajando en un algoritmo iterativo donde, para las primeras iteraciones iniciales, las matrices pueden ser dispersas pero gradualmente se vuelven densas a medida que pasa el tiempo.

Gracias por adelantado

Editado para incluir código para reproducir el escenario anterior

#include <iostream>
#include <stdio.h>
#include <time.h>
#include <cblas.h>
#include <cublas_v2.h>
using namespace std;

int main()
{
const int m = 5000; 

timespec blas_start, blas_end, cublas_start, cublas_end;
long totalnsec; //total nano sec
double totalsec, totaltime;
int i, j;
float *A = new float[m]; // 1 x m
float *B = new float[m*m]; // m x m
float *C = new float[m]; // 1 x m

// input martix A: every 32nd element is non-zero
for(i = 0; i < m; i++)
{       
    A[i] = 0;
    if( i % 32 == 0)    //adjust for sparsity
        A[i] = i;
}

// input matrix B: identity matrix
// col major = row major
for(i = 0; i < m; i++)
    for(j = 0; j < m; j++)
    {
        if (i==j)           
            B[j*m + i] = 1;
        else
            B[j*m + i] = 0;     
    }

clock_gettime(CLOCK_REALTIME, &blas_start);
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1, A, m, B, m, 0, C, m);
clock_gettime(CLOCK_REALTIME, &blas_end);
    /*
for(i = 0; i < 12; i++)
    printf("%f ", C[i]);
    */

//cublas section

cudaError_t cudaStat;   
cublasHandle_t handle;
cublasCreate(&handle);
//Declaring Device Variables
float *A_d, *B_d, *C_d;

//Allocating Memory for Device Variables
cudaStat = cudaMalloc(&A_d, sizeof(float)*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for A_d/n");

cudaStat = cudaMalloc(&B_d, sizeof(float)*m*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for B_d/n");

cudaStat = cudaMalloc(&C_d, sizeof(float)*m);
if(cudaStat != cudaSuccess) printf("Error Allocating Memory for C_d/n");

// Moving values of A, B onto Device variables
cublasSetVector(m, sizeof(float), A, 1, A_d, 1);
cublasSetMatrix(m, m, sizeof(float), B, m, B_d, m); 

// Do the actual multiplication
float alpha = 1.0f, beta = 0.0f;
cudaDeviceSynchronize();
clock_gettime(CLOCK_REALTIME, &cublas_start);

cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, 1, m, m, &alpha, A_d, 1, B_d, m, &beta, C_d, 1);

cudaDeviceSynchronize();
clock_gettime(CLOCK_REALTIME, &cublas_end);

cublasGetVector(m, sizeof(float), C, 1, C_d, 1);
    /*
for(i = 0; i < 12; i++)
    printf("%f ", C[i]);
    */

// Print times
// blas time
totalsec = (double)blas_end.tv_sec - (double)blas_start.tv_sec;
totalnsec = blas_end.tv_nsec - blas_start.tv_nsec;
if(totalnsec < 0)
{
    totalnsec += 1e9;
    totalsec -= 1;
}
totaltime = totalsec + (double)totalnsec*1e-9;
cout<<"BLAS Time = "<< totaltime << "/n";

//cublas
totalsec = (double)cublas_end.tv_sec - (double)cublas_start.tv_sec;
totalnsec = cublas_end.tv_nsec - cublas_start.tv_nsec;
if(totalnsec < 0)
{
    totalnsec += 1e9;
    totalsec -= 1;
}
totaltime = totalsec + (double)totalnsec*1e-9;
cout<<"CUBLAS Time = "<< totaltime << "/n";

return 0;
}

Lo ejecutó para obtener los siguientes resultados

malang@ubuntu:~/uas/stackoverflow$ nvcc -arch=sm_12 blascomp.cu -o blascomp.o -lblas -lcublas
malang@ubuntu:~/uas/stackoverflow$ ./blascomp.o
BLAS Time = 0.000964504
CUBLAS Time = 0.0365322

EDITAR

Editado después de la respuesta de @Eric

El uso de cublasSgemv ha mejorado enormemente el rendimiento en la GPU. Pero, todavía tengo este problema de que cblas_sgemm sea mucho más eficiente para matrices dispersas en la CPU. ¿Cuáles podrían ser las posibles razones?

EDITAR Ejecuta los siguientes comandos sobre la sugerencia de @Eric @osgx @Robert Crovella

erisp@ubuntu:~/uas/stackoverflow$ ldd ./gemmcomp.o
    linux-gate.so.1 =>  (0xb76f6000)
    libblas.so.3 => /usr/lib/libblas.so.3 (0xb765e000)
    libstdc++.so.6 => /usr/lib/i386-linux-gnu/libstdc++.so.6 (0xb7576000)
    libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb73c7000)
    libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7381000)
/lib/ld-linux.so.2 (0xb76f7000)
    libgcc_s.so.1 => /lib/i386-linux-gnu/libgcc_s.so.1 (0xb7364000)

erisp@ubuntu:~/uas/stackoverflow$ ll -d /usr/lib/libblas* /etc/alternatives/libblas.*
lrwxrwxrwx 1 root root   26 مارچ  13  2015 /etc/alternatives/libblas.a -> /usr/lib/libblas/libblas.a
lrwxrwxrwx 1 root root   27 مارچ  13  2015 /etc/alternatives/libblas.so -> /usr/lib/libblas/libblas.so
lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3 -> /usr/lib/libblas/libblas.so.3
lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3gf -> /usr/lib/libblas/libblas.so.3
drwxr-xr-x 2 root root 4096 مارچ  13  2015 /usr/lib/libblas/
lrwxrwxrwx 1 root root   27 مارچ  13  2015 /usr/lib/libblas.a -> /etc/alternatives/libblas.a
lrwxrwxrwx 1 root root   28 مارچ  13  2015 /usr/lib/libblas.so -> /etc/alternatives/libblas.so
lrwxrwxrwx 1 root root   30 مارچ  13  2015 /usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3
lrwxrwxrwx 1 root root   32 مارچ  13  2015 /usr/lib/libblas.so.3gf -> /etc/alternatives/libblas.so.3gf

Answer 1

Tu código tiene un problema: estás usando la API BLAS incorrecta. gemm() rutina de multiplicación matriz-matriz gemm() para hacer una operación de multiplicación de matriz vectorial.

Para vec-mat-mul o mat-vec-mul debes usar gemv() . Por supuesto, gemm() puede dar el resultado correcto con una matriz que tiene solo 1 fila. Pero este es un caso de esquina inesperado que gemv() debería manejar, por lo que es posible que no obtenga el máximo rendimiento en GPU y / o CPU.

Podrías cambiar a gemv() y benchmark de nuevo.

EDITAR

Aquí está mi resultado de referencia con un solo hilo MKL. Los valores de A y B son los mismos que en su código. No puedo reproducir el resultado de ''0.000964504s'' en la CPU. Puede verificar la corrección de su vector de resultados. Existe la posibilidad de que tu biblioteca cblas tenga un error.

Usando gemm()

BLAS Time = 0.0169784
CUBLAS Time = 0.00356155

Usando gemv()

BLAS Time = 0.0167557
CUBLAS Time = 0.0013809

EDIT2

Ahora puedo reproducir el resultado de FAST en unbuntu 14.04 con el paquete libblas-dev .

El motivo se responde en la siguiente pregunta.

cblas gemm tiempo depende de los valores de la matriz de entrada - Ubuntu 14.04

En la versión particular de BLAS, hay código para verificar cero elemento. El costo de comprobación es O (n ^ 2), por lo que vale la pena hacerlo en la multiplicación matriz-matriz cuyo costo es O (n ^ 3).

Para GPU gemm (), como el orden de la computación es diferente (bloque por bloque en lugar de línea por línea), dicha optimización puede no ser factible. Pero es factible para GPU gemv (), donde podría guardarse el tiempo dedicado a cargar la matriz desde la memoria global.