首页 文章

CUDA / CUBLAS矩阵 - 向量乘法

提问于
浏览
1

我之前发布了一个关于CUDA中矩阵向量乘法的问题以及关于编写自己的内核的问题 . 在这之后,我决定使用CUBLAS按照一些用户(感谢@Robert Crovella)的建议来实现我的问题,希望获得更高的性能(我的项目是性能驱动的) .

只是为了澄清:我想将NxN矩阵与1xN向量相乘 .

我一直在看下面粘贴的代码几天,我无法弄清楚为什么乘法会给我一个不正确的结果 . 我担心使用<vector>数组导致问题(这是使用这些数据类型的更大系统的一部分) . 我不是故意将此线程用作调试工具,但我认为这也会对其他试图实现此目的的用户有所帮助,因为我没有在互联网上遇到特别全面的问题来解决我的问题(以及对于我的问题) v2 API) . 提前致谢!

#include <cuda.h>
#include <vector>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include <cublas_v2.h>
#include <time.h>

//#include "timenow.cu"

// error check macros
#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

// for CUBLAS V2 API
#define cublasCheckErrors(fn) \
    do { \
        cublasStatus_t __err = fn; \
        if (__err != CUBLAS_STATUS_SUCCESS) { \
            fprintf(stderr, "Fatal cublas error: %d (at %s:%d)\n", \
                (int)(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)

// random data filler
void fillvector(float *data, int N){
    for(int i=0; i<N; i++){
        data[i] = float(rand() % 10);
    }
}

//printer
void printer(bool printOut, float *data, int N){
    if(printOut == true){
    for(int i=0; i<N; i++){
        printf("%2.1f ", data[i]);
    }
    printf("\n");
    }
}

/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////

int main(){

bool printOut = true;

int N;
std::cout << "Enter N: " ;
std::cin >> N;

std::vector<float> x0;
x0.resize(N);

std::vector<float> p;
p.resize(N*N);

// matrix A
std::vector<float> A[N];
for(int i=0;i<N;i++){
        A[i].resize(N);
    fillvector(A[i].data(), N);
    printer(printOut, A[i].data(), N);
}
printf("\n");
fillvector(x0.data(), N);
printer(printOut, x0.data(), N);

printf("\nStarting CUDA computation...");
///double startTime = timenow();

// device pointers
float *d_A, *d_p, *d_b, *d_x0, *d_v, *d_temp;

cudaMalloc((void**)&d_A, N*N*sizeof(float));
cudaMalloc((void**)&d_temp, N*sizeof(float));
cudaMalloc((void**)&d_x0, N*sizeof(float));
cudaCheckErrors("cuda malloc fail");

// might need to flatten A...
cublasSetVector(N, sizeof(float), &x0, 1, d_x0, 1);
//daMemcpy(d_x0, &x0, N*sizeof(float), cudaMemcpyHostToDevice);
cublasSetMatrix(N, N, sizeof(float), &A, N, d_A, N);
cudaCheckErrors("cuda memcpy of A or x0 fail");

float *temp;
temp = (float *)malloc(N*sizeof(temp));

cublasHandle_t handle;
cublasCheckErrors(cublasCreate(&handle));

float alpha = 1.0f;
float beta = 0.0f;
cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_N, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1));

cublasGetVector(N, sizeof(float), &temp, 1, d_temp, 1);
//cudaMemcpy(temp, d_temp, N*sizeof(float), cudaMemcpyDeviceToHost);
cudaCheckErrors("returning to host failed");

printf("\n");
printer(printOut, temp, N);

/*alpha = -1.0;
cublasSaxpy(handle, N, &alpha, d_temp, 1, d_v, 1);
cublasGetVector(N, sizeof(float) * N, d_v, 1, &v, 1);
printf("\n");
for(int i=0; i<N; i++){
    printf("%2.1f ",v[i]);
}*/

printf("\nFinished CUDA computations...");
//double endTime = timenow();

//double timeDiff = endTime - startTime;
//printf("\nRuntime: %2.3f seconds \n", timeDiff);

cudaFree(d_temp);
cudaFree(d_A);
cudaFree(d_p);
cudaFree(d_x0);

return 0;
}

1 回答

  • 2
    • 我们不以这种方式引用向量的第一个元素:
    cublasSetVector(N, sizeof(float), &x0, 1, d_x0,
    

    相反,你应该这样做:

    cublasSetVector(N, sizeof(float), &(x0[0]), 1, d_x0, 1);
    

    同样适合您的 SetMatrix 呼叫引用 A

    cublasSetMatrix(N, N, sizeof(float), &(A[0]), N, d_A, N);
    
    • 您的 GetVector 电话有2个错误:
    cublasGetVector(N, sizeof(float), &temp, 1, d_temp, 1);
    

    您有 tempd_temp 参数reversed(您正在从设备复制到主机)并且您不应该使用 temp 的地址:它已经是指针 . 这样做:

    cublasGetVector(N, sizeof(float), d_temp, 1, temp, 1);
    
    • 您没有对所有cublas调用进行正确的错误检查,例如get / set matrix / vector调用 . 使用与其他Cublas调用相同的方法 .

    • 您正在创建 A 作为向量数组 . 这不适用于 cublasSetMatrix . 相反,我们需要创建 A 作为平面向量,具有足够的大小(N * N)来存储整个矩阵 .

    • 最后,cublas希望它使用的矩阵按列主要顺序存储 . 如果以行主顺序传递C样式数组,则应在 cublasSgemv 中使用该矩阵的转置:

    cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_T, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1));
    

    以下代码修复了以下各种问题:

    $ cat t235.cu
    #include <cuda.h>
    #include <vector>
    #include <iostream>
    #include <fstream>
    #include <stdio.h>
    #include <stdlib.h>
    #include <cmath>
    #include <cublas_v2.h>
    #include <time.h>
    
    //#include "timenow.cu"
    
    // error check macros
    #define cudaCheckErrors(msg) \
        do { \
            cudaError_t __err = cudaGetLastError(); \
            if (__err != cudaSuccess) { \
                fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                    msg, cudaGetErrorString(__err), \
                    __FILE__, __LINE__); \
                fprintf(stderr, "*** FAILED - ABORTING\n"); \
                exit(1); \
            } \
        } while (0)
    
    // for CUBLAS V2 API
    #define cublasCheckErrors(fn) \
        do { \
            cublasStatus_t __err = fn; \
            if (__err != CUBLAS_STATUS_SUCCESS) { \
                fprintf(stderr, "Fatal cublas error: %d (at %s:%d)\n", \
                    (int)(__err), \
                    __FILE__, __LINE__); \
                fprintf(stderr, "*** FAILED - ABORTING\n"); \
                exit(1); \
            } \
        } while (0)
    
    // random data filler
    void fillvector(float *data, int N){
        for(int i=0; i<N; i++){
            data[i] = float(rand() % 10);
        }
    }
    
    //printer
    void printer(bool printOut, float *data, int N){
        if(printOut == true){
        for(int i=0; i<N; i++){
            printf("%2.1f ", data[i]);
        }
        printf("\n");
        }
    }
    
    /////////////////////////////////////////////////////////////////////
    /////////////////////////////////////////////////////////////////////
    
    int main(){
    
    bool printOut = true;
    
    int N;
    std::cout << "Enter N: " ;
    std::cin >> N;
    
    std::vector<float> x0;
    x0.resize(N);
    
    std::vector<float> p;
    p.resize(N*N);
    
    // matrix A
    std::vector<float> A(N*N);
    fillvector(A.data(), N*N);
    for (int i=0; i< N; i++){
      printer(printOut, &(A[(i*N)]), N);
      printf("\n");}
    fillvector(x0.data(), N);
    printer(printOut, x0.data(), N);
    
    printf("\nStarting CUDA computation...");
    ///double startTime = timenow();
    
    // device pointers
    float *d_A, *d_x0, *d_temp;
    
    cudaMalloc((void**)&d_A, N*N*sizeof(float));
    cudaMalloc((void**)&d_temp, N*sizeof(float));
    cudaMalloc((void**)&d_x0, N*sizeof(float));
    cudaCheckErrors("cuda malloc fail");
    
    // might need to flatten A...
    cublasCheckErrors(cublasSetVector(N, sizeof(float), &(x0[0]), 1, d_x0, 1));
    //daMemcpy(d_x0, &x0, N*sizeof(float), cudaMemcpyHostToDevice);
    cublasCheckErrors(cublasSetMatrix(N, N, sizeof(float), &(A[0]), N, d_A, N));
    //cudaCheckErrors("cuda memcpy of A or x0 fail");
    
    float *temp;
    temp = (float *)malloc(N*sizeof(temp));
    
    cublasHandle_t handle;
    cublasCheckErrors(cublasCreate(&handle));
    
    float alpha = 1.0f;
    float beta = 0.0f;
    cublasCheckErrors(cublasSgemv(handle, CUBLAS_OP_T, N, N, &alpha, d_A, N, d_x0, 1, &beta, d_temp, 1));
    
    cublasCheckErrors(cublasGetVector(N, sizeof(float), d_temp, 1, temp, 1));
    //cudaMemcpy(temp, d_temp, N*sizeof(float), cudaMemcpyDeviceToHost);
    //cudaCheckErrors("returning to host failed");
    
    printf("\n");
    printer(printOut, temp, N);
    
    /*alpha = -1.0;
    cublasSaxpy(handle, N, &alpha, d_temp, 1, d_v, 1);
    cublasGetVector(N, sizeof(float) * N, d_v, 1, &v, 1);
    printf("\n");
    for(int i=0; i<N; i++){
        printf("%2.1f ",v[i]);
    }*/
    
    printf("\nFinished CUDA computations...\n");
    //double endTime = timenow();
    
    //double timeDiff = endTime - startTime;
    //printf("\nRuntime: %2.3f seconds \n", timeDiff);
    
    cudaFree(d_temp);
    cudaFree(d_A);
    //cudaFree(d_p);
    cudaFree(d_x0);
    
    return 0;
    }
    $ nvcc -arch=sm_20 -O3 -o t235 t235.cu -lcublas
    $ ./t235
    Enter N: 5
    3.0 6.0 7.0 5.0 3.0
    
    5.0 6.0 2.0 9.0 1.0
    
    2.0 7.0 0.0 9.0 3.0
    
    6.0 0.0 6.0 2.0 6.0
    
    1.0 8.0 7.0 9.0 2.0
    
    0.0 2.0 3.0 7.0 5.0
    
    Starting CUDA computation...
    83.0 86.0 92.0 62.0 110.0
    
    Finished CUDA computations...
    $
    

相关问题