cuda中用cublas庫做矩陣乘法

這里矩陣C=A*B，原始文檔給的公式是C=alpha*A*B+beta*C，所以這里alpha=1，beta=0。

主要使用cublasSgemm這個函數，這個函數的第二個參數有三種類型，這里CUBLAS_OP_N求出來矩陣結果是按行排列，所以不需要轉置了。

如果用CUBLAS_OP_T參數求得的結果是按列排列，做成C形式的矩陣應該還需要轉置一下，並且后面跟的參數也不太一樣，這個參數我就沒再嘗試了。

代碼如下：

#include "cuda_runtime.h"
#include "cublas_v2.h"

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <ctime>

using namespace std;

int main()
{
    srand(time(0));
    int M = 2;            //矩陣A的行，矩陣C的行
    int N = 3;            //矩陣A的列，矩陣B的行
    int K = 4;            //矩陣B的列，矩陣C的列

    float *h_A = (float*)malloc(sizeof(float)*M*N);
    float *h_B = (float*)malloc(sizeof(float)*N*K);
    float *h_C = (float*)malloc(sizeof(float)*M*K);

    for (int i = 0; i < M*N; i++)
    {
        h_A[i] = rand() % 10;
        cout << h_A[i] << "  ";
        if ((i + 1) % N == 0)
            cout << endl;        
    }
    cout << endl;

    for (int i = 0; i < N*K; i++)
    {
        h_B[i] = rand() % 10;
        cout << h_B[i] << "  ";
        if ((i + 1) % K == 0)
            cout << endl;
    }
    cout << endl;

    float *d_A, *d_B, *d_C,*d_CT;
    cudaMalloc((void**)&d_A, sizeof(float)*M*N);
    cudaMalloc((void**)&d_B, sizeof(float)*N*K);
    cudaMalloc((void**)&d_C, sizeof(float)*M*K);

    cudaMemcpy(d_A, h_A, M*N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, N*K * sizeof(float), cudaMemcpyHostToDevice);

    float alpha = 1;
    float beta = 0;

    //C=A*B
    cublasHandle_t handle;
    cublasCreate(&handle);
    cublasSgemm(handle,
        CUBLAS_OP_N,  
        CUBLAS_OP_N,   
        K,                    //矩陣B的列數
        M,                    //矩陣A的行數
        N,                    //矩陣A的列數
        &alpha,           
        d_B,            
        K,                    
        d_A,         
        N,         
        &beta,          
        d_C,           
        K);

    cudaMemcpy(h_C, d_C, M*K * sizeof(float), cudaMemcpyDeviceToHost);

    for (int i = 0; i < M*K; i++)
    {
        cout << h_C[i] << "  ";
        if ((i+1)%K==0)
            cout << endl;
    }

    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    free(h_A);
    free(h_B);
    free(h_C);
    return 0;
}

結果：