quadro - usar gpu para matlab

CUDA para resolver muchos sistemas lineales "pequeños/moderados"

Question

Alguna información de fondo sobre el problema que estoy tratando de acelerar usando CUDA:

Tengo una gran cantidad de sistemas lineales del mismo tamaño pequeño / moderado que necesito resolver de forma independiente. Cada sistema lineal es cuadrado, real, denso, invertible y no simétrico. En realidad, estos son sistemas matriciales, por lo que cada sistema se ve como AX = B, donde A, X y B son matrices (nxn).

En esta pregunta previa, solicito los tamaños de lotes y matrices de CUBLAS , donde aprendo las operaciones por lotes de CuBLAS, obtengo el mejor rendimiento para matrices de tamaño 100x100 o menores.

Todavía tengo un problema porque las matrices con las que estoy trabajando tienen 100 <n <700. Por lo tanto, las matrices son de tamaño moderado donde las operaciones por lotes cuBLAS no ofrecen el mejor rendimiento, y las BLAS regulares (cusolverDnDgetrf, cusolverDnDgetrs) tampoco son mejores. rendimiento que MATLAB (mira los tiempos a continuación).

Hice algunos tiempos en comparación con MATLAB, para resolver un solo sistema, y ​​encontré BLAS regular es mejor para matrices de tamaño (4096x4096) o más grande. Hago una matriz aleatoria de tamaño (nxn), para n = 64,256,512,1024,4096,16384, y solo el tiempo de la factorización y la resolución de retroceso / avance, no se transfieren a través de PCIE.

DOBLE PRECISIÓN CUDA (GTX 1080ti) vs MATLAB (barra invertida)

(GPU) 64: 0.001157 sec (MATLAB) 64: 0.000205 sec

(GPU) 256: 0.01161 seg (MATLAB) 256: 0.007762 sec

(GPU) 512: 0.026348 seg (MATLAB) 512: 0.008550 sec

(GPU) 1024: 0.064357 seg (MATLAB) 1024: 0.036280 sec

(GPU) 4096: 0.734908 seg (MATLAB) 4096: 1.174442 sec

(GPU) 16384: 32.962229 seg (MATLAB) 16384: 68.691236 seg

Estos tiempos me hacen concluir que iterar uno por uno sobre mis matrices llamando al método de inversión no por lotes será más lento que MATLAB. Además, para mis matrices de tamaño moderado, el método de inversión por lotes discontinuo cuBLAS no tendrá un buen rendimiento, de acuerdo con los tamaños de lotes y matrices de CUBLAS .

¿Hay algún otro enfoque que deba considerar para acelerar mi código con CUDA? ¿O estoy malentendiendo algo?

 /*  How to use
 *      ./cuSolverDn_LinearSolver                     // Default: cholesky
 *     ./cuSolverDn_LinearSolver -R=chol -filefile>   // cholesky factorization
 *     ./cuSolverDn_LinearSolver -R=lu -file<file>     // LU with partial pivoting
 *     ./cuSolverDn_LinearSolver -R=qr -file<file>     // QR factorization
 *
 *  Remark: the absolute error on solution x is meaningless without knowing condition number of A.
 *     The relative error on residual should be close to machine zero, i.e. 1.e-15.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>

#include <cuda_runtime.h>

#include "cublas_v2.h"
#include "cusolverDn.h"
#include "helper_cuda.h"

#include "helper_cusolver.h"
int linearSolverLU(
    cusolverDnHandle_t handle,
    int n,
    const double *Acopy,
    int lda,
    const double *b,
    double *x)
{
    int bufferSize = 0;
    int *info = NULL;
    double *buffer = NULL;
    double *A = NULL;
    int *ipiv = NULL; // pivoting sequence
    int h_info = 0;
    double start, stop;
    double time_solve;

    checkCudaErrors(cusolverDnDgetrf_bufferSize(handle, n, n, (double*)Acopy, lda, &bufferSize));

    checkCudaErrors(cudaMalloc(&info, sizeof(int)));
    checkCudaErrors(cudaMalloc(&buffer, sizeof(double)*bufferSize));
    checkCudaErrors(cudaMalloc(&A, sizeof(double)*lda*n));
    checkCudaErrors(cudaMalloc(&ipiv, sizeof(int)*n));


    // prepare a copy of A because getrf will overwrite A with L
    checkCudaErrors(cudaMemcpy(A, Acopy, sizeof(double)*lda*n, cudaMemcpyDeviceToDevice));
    checkCudaErrors(cudaMemset(info, 0, sizeof(int)));

    start = second();
    start = second();

    checkCudaErrors(cusolverDnDgetrf(handle, n, n, A, lda, buffer, ipiv, info));
    checkCudaErrors(cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));

    if ( 0 != h_info ){
        fprintf(stderr, "Error: LU factorization failed/n");
    }

    //checkCudaErrors(cudaMemcpy(x, b, sizeof(double)*n, cudaMemcpyDeviceToDevice));
    checkCudaErrors(cudaMemcpy(x, b, sizeof(double)*lda*n, cudaMemcpyDeviceToDevice));
    //checkCudaErrors(cusolverDnDgetrs(handle, CUBLAS_OP_N, n, 1, A, lda, ipiv, x, n, info));
    checkCudaErrors(cusolverDnDgetrs(handle, CUBLAS_OP_N, n, n, A, lda, ipiv, x, n, info));
    checkCudaErrors(cudaDeviceSynchronize());
    stop = second();

    time_solve = stop - start;
    fprintf (stdout, "timing: LU = %10.6f sec/n", time_solve);

    if (info  ) { checkCudaErrors(cudaFree(info  )); }
    if (buffer) { checkCudaErrors(cudaFree(buffer)); }
    if (A     ) { checkCudaErrors(cudaFree(A)); }
    if (ipiv  ) { checkCudaErrors(cudaFree(ipiv));}

    return 0;
}
void generate_random_dense_matrix(int M, int N, double **outA)
{
    int i, j;
    double rMax = (double)RAND_MAX;
    double *A   = (double *)malloc(sizeof(double) * M * N);

    // For each column
    for (j = 0; j < N; j++)
    {
        // For each row
        for (i = 0; i < M; i++)
        {
            double dr = (double)rand();
            A[j * M + i] = (dr / rMax) * 100.0;
        //printf("A[j * M + i] = %f /n",A[j * M + i]);
        }
    }

    *outA = A;
}

int main (int argc, char *argv[])
{
    struct testOpts opts;
    cusolverDnHandle_t handle = NULL;
    cublasHandle_t cublasHandle = NULL; // used in residual evaluation
    cudaStream_t stream = NULL;

    int rowsA = 0; // number of rows of A
    int colsA = 0; // number of columns of A
    int nnzA  = 0; // number of nonzeros of A
    int baseA = 0; // base index in CSR format
    int lda   = 0; // leading dimension in dense matrix

    // CSR(A) from I/O
    int *h_csrRowPtrA = NULL;
    int *h_csrColIndA = NULL;
    double *h_csrValA = NULL;

    double *h_A = NULL; // dense matrix from CSR(A)
    double *h_x = NULL; // a copy of d_x
    double *h_b = NULL; // b = ones(m,1)
    double *h_r = NULL; // r = b - A*x, a copy of d_r

    double *d_A = NULL; // a copy of h_A
    double *d_x = NULL; // x = A / b
    double *d_b = NULL; // a copy of h_b
    double *d_r = NULL; // r = b - A*x

    // the constants are used in residual evaluation, r = b - A*x
    const double minus_one = -1.0;
    const double one = 1.0;

    double x_inf = 0.0;
    double r_inf = 0.0;
    double A_inf = 0.0;
    int errors = 0;


    colsA = 660;
    rowsA = colsA;
    int NN = colsA;
    int MM = rowsA;
    lda = rowsA;

    // Generate inputs
    srand(9384);  
    generate_random_dense_matrix(MM, NN, &h_A);
    generate_random_dense_matrix(MM, NN, &h_b);

    parseCommandLineArguments(argc, argv, opts);

    if (NULL == opts.testFunc)
    {
        //opts.testFunc = "chol"; // By default running Cholesky as NO solver selected with -R option.
    opts.testFunc = "lu";  
    //opts.testFunc = "qr"; 
    }

    findCudaDevice(argc, (const char **)argv);

    /*
    printf("step 1: read matrix market format/n");

    if (opts.sparse_mat_filename == NULL)
    {
        opts.sparse_mat_filename =  sdkFindFilePath("gr_900_900_crg.mtx", argv[0]);
        if (opts.sparse_mat_filename != NULL)
            printf("Using default input file [%s]/n", opts.sparse_mat_filename);
        else
            printf("Could not find gr_900_900_crg.mtx/n");
    }
    else
    {
        printf("Using input file [%s]/n", opts.sparse_mat_filename);
    }

    if (opts.sparse_mat_filename == NULL)
    {
        fprintf(stderr, "Error: input matrix is not provided/n");
        return EXIT_FAILURE;
    }

    if (loadMMSparseMatrix<double>(opts.sparse_mat_filename, ''d'', true , &rowsA, &colsA,
               &nnzA, &h_csrValA, &h_csrRowPtrA, &h_csrColIndA, true))
    {
        exit(EXIT_FAILURE);
    }
    baseA = h_csrRowPtrA[0]; // baseA = {0,1}

    printf("sparse matrix A is %d x %d with %d nonzeros, base=%d/n", rowsA, colsA, nnzA, baseA);

    if ( rowsA != colsA )
    {
        fprintf(stderr, "Error: only support square matrix/n");
        exit(EXIT_FAILURE);
    }

    printf("step 2: convert CSR(A) to dense matrix/n");

    lda = opts.lda ? opts.lda : rowsA;
    if (lda < rowsA)
    {
        fprintf(stderr, "Error: lda must be greater or equal to dimension of A/n");
        exit(EXIT_FAILURE);
    }
    */
    //h_A = (double*)malloc(sizeof(double)*lda*colsA);
    h_x = (double*)malloc(sizeof(double)*lda*colsA);
    //h_b = (double*)malloc(sizeof(double)*rowsA);
    h_r = (double*)malloc(sizeof(double)*lda*rowsA);
    assert(NULL != h_A);
    assert(NULL != h_x);
    assert(NULL != h_b);
    assert(NULL != h_r);

    /*
    memset(h_A, 0, sizeof(double)*lda*colsA); 
    for(int row = 0 ; row < rowsA ; row++)
    {
        const int start = h_csrRowPtrA[row  ] - baseA;
        const int end   = h_csrRowPtrA[row+1] - baseA;
        for(int colidx = start ; colidx < end ; colidx++)
        {
            const int col = h_csrColIndA[colidx] - baseA;
            const double Areg = h_csrValA[colidx];
            h_A[row + col*lda] = Areg;
        }
    }

    printf("step 3: set right hand side vector (b) to 1/n");
    for(int row = 0 ; row < rowsA ; row++)
    {
        h_b[row] = 1.0;
    }
    */
    // verify if A is symmetric or not.
    if ( 0 == strcmp(opts.testFunc, "chol") )
    {
        int issym = 1;
        for(int j = 0 ; j < colsA ; j++)
        {
            for(int i = j ; i < rowsA ; i++)
            {
                double Aij = h_A[i + j*lda];
                double Aji = h_A[j + i*lda];
                if ( Aij != Aji )
                {
                    issym = 0;
                    break;
                }
            }
        }
        if (!issym)
        {
            printf("Error: A has no symmetric pattern, please use LU or QR /n");
            exit(EXIT_FAILURE);
        }
    }

    checkCudaErrors(cusolverDnCreate(&handle));
    checkCudaErrors(cublasCreate(&cublasHandle));
    checkCudaErrors(cudaStreamCreate(&stream));

    checkCudaErrors(cusolverDnSetStream(handle, stream));
    checkCudaErrors(cublasSetStream(cublasHandle, stream));


    checkCudaErrors(cudaMalloc((void **)&d_A, sizeof(double)*lda*colsA));
    checkCudaErrors(cudaMalloc((void **)&d_x, sizeof(double)*lda*colsA));
    checkCudaErrors(cudaMalloc((void **)&d_b, sizeof(double)*lda*rowsA));
    checkCudaErrors(cudaMalloc((void **)&d_r, sizeof(double)*lda*rowsA));

    printf("step 4: prepare data on device/n");
    checkCudaErrors(cudaMemcpy(d_A, h_A, sizeof(double)*lda*colsA, cudaMemcpyHostToDevice));
    checkCudaErrors(cudaMemcpy(d_b, h_b, sizeof(double)*lda*rowsA, cudaMemcpyHostToDevice));

    printf("step 5: solve A*x = b /n");
    // d_A and d_b are read-only
    if ( 0 == strcmp(opts.testFunc, "chol") )
    {
        linearSolverCHOL(handle, rowsA, d_A, lda, d_b, d_x);
    }
    else if ( 0 == strcmp(opts.testFunc, "lu") )
    {
        //printf("hi /n");
        linearSolverLU(handle, rowsA, d_A, lda, d_b, d_x);
    }
    else if ( 0 == strcmp(opts.testFunc, "qr") )
    {
        linearSolverQR(handle, rowsA, d_A, lda, d_b, d_x);
    }
    else
    {
        fprintf(stderr, "Error: %s is unknown function/n", opts.testFunc);
        exit(EXIT_FAILURE);
    }
    printf("step 6: evaluate residual/n");
    checkCudaErrors(cudaMemcpy(d_r, d_b, sizeof(double)*lda*rowsA, cudaMemcpyDeviceToDevice));

    // r = b - A*x
    checkCudaErrors(cublasDgemm_v2(
        cublasHandle,
        CUBLAS_OP_N,
        CUBLAS_OP_N,
        rowsA,
        colsA,
        colsA,
        &minus_one,
        d_A,
        lda,
        d_x,
        rowsA,
        &one,
        d_r,
        rowsA));

    checkCudaErrors(cudaMemcpy(h_x, d_x, sizeof(double)*lda*colsA, cudaMemcpyDeviceToHost));
    checkCudaErrors(cudaMemcpy(h_r, d_r, sizeof(double)*lda*rowsA, cudaMemcpyDeviceToHost));

    x_inf = vec_norminf(colsA, h_x);
    r_inf = vec_norminf(rowsA, h_r);
    A_inf = mat_norminf(rowsA, colsA, h_A, lda);


    printf("x[0] = %f/n", h_x[0]);
    printf("r[0] = %f/n", h_r[0]);
    printf("|b - A*x| = %E /n", r_inf);
    printf("|A| = %E /n", A_inf);
    printf("|x| = %E /n", x_inf);
    printf("|b - A*x|/(|A|*|x|) = %E /n", r_inf/(A_inf * x_inf));

    if (handle) { checkCudaErrors(cusolverDnDestroy(handle)); }
    if (cublasHandle) { checkCudaErrors(cublasDestroy(cublasHandle)); }
    if (stream) { checkCudaErrors(cudaStreamDestroy(stream)); }

    if (h_csrValA   ) { free(h_csrValA); }
    if (h_csrRowPtrA) { free(h_csrRowPtrA); }
    if (h_csrColIndA) { free(h_csrColIndA); }

    if (h_A) { free(h_A); }
    if (h_x) { free(h_x); }
    if (h_b) { free(h_b); }
    if (h_r) { free(h_r); }

    if (d_A) { checkCudaErrors(cudaFree(d_A)); }
    if (d_x) { checkCudaErrors(cudaFree(d_x)); }
    if (d_b) { checkCudaErrors(cudaFree(d_b)); }
    if (d_r) { checkCudaErrors(cudaFree(d_r)); }

    return 0;
}

Answer 1

Intente utilizar dos o más transmisiones en paralelo (con un sistema lineal cada una) en la GPU, posiblemente esto ayude a utilizar una parte más grande de la GPU.

Para medir el tiempo y la utilización del hardware, utilice el generador de perfiles visual en lugar de las mediciones de tiempo de la CPU.

Otro punto es que las GPU GTX (del consumidor) funcionan bastante mal en la doble presición. Si tienes la oportunidad, intenta usar una GPU Tesla en su lugar.

Answer 2

MATLAB proporciona una forma de llamar a la interfaz por lotes de cublas para matrices de GPU utilizando pagefun .