用cublas实现行优先矩阵乘法和列优先矩阵乘法

最新推荐文章于 2024-03-12 18:10:56 发布

lalalala256

最新推荐文章于 2024-03-12 18:10:56 发布

阅读量2.7k

点赞数 4

分类专栏： GPU

本文链接：https://blog.csdn.net/qq_25147897/article/details/70806487

版权

GPU 专栏收录该内容

15 篇文章 1 订阅

订阅专栏

引用：http://www.cnblogs.com/scut-fm/p/3756242.html
cublas库是CUDA标准的线代库，但没有专门针对稀疏矩阵的操作。
其中cublasSgemm实现C=α*A*B+β*C功能
函数原型为

/* GEMM */
CUBLASAPI cublasStatus_t CUBLASWINAPI cublasSgemm_v2 (cublasHandle_t handle, 
                                                      cublasOperation_t transa,
                                                      cublasOperation_t transb, 
                                                      int m,
                                                      int n,
                                                      int k,
                                                      const float *alpha, /* host or device pointer */  
                                                      const float *A, 
                                                      int lda,
                                                      const float *B,
                                                      int ldb, 
                                                      const float *beta, /* host or device pointer */  
                                                      float *C,
                                                      int ldc);

函数中handle为blas库对象
transa 为矩阵A属性参数
transb 为矩阵B属性参数
M A，C行数，即C的高
N B，C列数，即C的宽
K A的列数和B的行数
alpha 运算式的α的地址
A A在显存中的地址
lda A的主维度，所表达的矩阵的高
B B在显存中的地址
ldb B的主维度，所表达的矩阵的高
beta 运算式的β的地址
C C在显存中的地址
ldc C的主维度，所表达的矩阵的高

若列优先计算矩阵乘法
|1 0| * |2 0 3| = |2 0 3|
|1 1 0|
A * B = C
则*alpha=1.0，*β=0.0
列存储形式：
A={1,0}表示矩阵，矩阵宽2，高1
1 0
B={2,1,0,1,3,0}，矩阵宽3，高2
2 0 3
1 1 0
C={2,0,3}，矩阵宽3，高1
2 0 3
ABC均按列存储，所以AB属性设置为CUBLAS_OP_N，得到的C矩阵也是按列存储
MNK分别为1,3,2
lda，ldb，ldc分别为三个矩阵的高，即1,2,1
全部代码如下

// CUDA runtime 库 + CUBLAS 库
#include "cuda_runtime.h"
#include "cuda_runtime_api.h"
#include "cublas_v2.h"

#include <time.h>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "MyTimer.h"
//本代码实现了cublas的列存储矩阵乘法
//#include <common_functions.h>
/*
|1 0| * |2 0 3| = |2 0 3|
        |1 1 0|
A     * B       = C
cublasSgemm 实现的功能为C = alpha*A*B + beta*C
ABC均为列存储
*/
using namespace std;

// 定义测试矩阵的维度
int const aW = 2;
int const aH = 1;
int const bW = 3;
int const bH = 2;
int const cW = bW;
int const cH = aH;

int main()
{
    // 定义状态变量
    cublasStatus_t status;

    // 在 内存 中为将要计算的矩阵开辟空间
    float *h_A = (float*)malloc (aW*aH*sizeof(float));
    float *h_B = (float*)malloc (bW*bH*sizeof(float));

    // 在 内存 中为将要存放运算结果的矩阵开辟空间
    float *h_C = (float*)malloc (cW*cH*sizeof(float));

    // 为待运算矩阵的元素赋值
    h_A[0]=1; h_A[1]=0;
    h_B[0]=2; h_B[1]=1; h_B[2]=0; h_B[3]=1; h_B[4]=3; h_B[5]=0;

    // 打印待测试的矩阵
    cout << "矩阵 A :" << endl;
    for (int i=0; i<aW*aH; i++){
        cout << h_A[i] << " ";
        if ((i+1)%aW == 0) cout << endl;
    }
    cout << endl;
    cout << "矩阵 B :" << endl;//按照行优先打印的，不正确
    for (int i=0; i<bW*bH; i++){
        cout << h_B[i] << " ";
        if ((i+1)%bW == 0) cout << endl;
    }
    cout << endl;

    /*
    ** GPU 计算矩阵相乘
    */

    // 创建并初始化 CUBLAS 库对象
    cublasHandle_t handle;
    status = cublasCreate(&handle);

    if (status != CUBLAS_STATUS_SUCCESS)
    {
        if (status == CUBLAS_STATUS_NOT_INITIALIZED) {
            cout << "CUBLAS 对象实例化出错" << endl;
        }
        getchar ();
        return EXIT_FAILURE;
    }

    float *d_A, *d_B, *d_C;
    // 在 显存 中为将要计算的矩阵开辟空间
    cudaMalloc (
        (void**)&d_A,    // 指向开辟的空间的指针
        aW*aH * sizeof(float)    //　需要开辟空间的字节数
    );
    cudaMalloc (
        (void**)&d_B,
        bW*bH * sizeof(float)
    );

    // 在 显存 中为将要存放运算结果的矩阵开辟空间
    cudaMalloc (
        (void**)&d_C,
        cW*cH * sizeof(float)
    );

    // 将矩阵数据传递进 显存 中已经开辟好了的空间
    cublasSetVector (
        aW*aH,    // 要存入显存的元素个数
        sizeof(float),    // 每个元素大小
        h_A,    // 主机端起始地址
        1,    // 连续元素之间的存储间隔
        d_A,    // GPU 端起始地址
        1    // 连续元素之间的存储间隔
    );
    cublasSetVector (
        bW*bH,
        sizeof(float),
        h_B,
        1,
        d_B,
        1
    );

    // 同步函数
    cudaThreadSynchronize();

    Timer myTimer;
    myTimer.start();
    // 传递进矩阵相乘函数中的参数，具体含义请参考函数手册。
    float a=1; float b=0;
    // 矩阵相乘。该函数必然将数组解析成列优先数组
    cublasSgemm (
        handle,    // blas 库对象
        CUBLAS_OP_N,    // 矩阵 A 属性参数
        CUBLAS_OP_N,    // 矩阵 B 属性参数
        cH,    // A, C 的行数
        cW,    // B, C 的列数
        aW,    // A 的列数和 B 的行数
        &a,    // 运算式的 α 值
        d_A,    // A 在显存中的地址
        aH,    // lda
        d_B,    // B 在显存中的地址
        bH,    // ldb
        &b,    // 运算式的 β 值
        d_C,    // C 在显存中的地址(结果矩阵)
        cH    // ldc
    );
    // 同步函数
    cudaThreadSynchronize();
    myTimer.stop();
    printf("use: %lf\n", myTimer.getElapsedTime());

    // 从 显存 中取出运算结果至 内存中去
    cublasGetVector (
        cW*cH,    //  要取出元素的个数
        sizeof(float),    // 每个元素大小
        d_C,    // GPU 端起始地址
        1,    // 连续元素之间的存储间隔
        h_C,    // 主机端起始地址
        1    // 连续元素之间的存储间隔
    );

    // 打印运算结果
    cout << "计算结果的转置 ( (A*B)的转置 )：" << endl;

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

    // 清理掉使用过的内存
    free (h_A);
    free (h_B);
    free (h_C);
    cudaFree (d_A);
    cudaFree (d_B);
    cudaFree (d_C);

    // 释放 CUBLAS 库对象
    cublasDestroy (handle);

    getchar();

    return 0;
}

C语言中行优先是更加习惯的使用，因此，在介绍下使用cublas实现行优先矩阵乘法的方法
C=A*B则C'=B'*A'
行存储形式：
A={1,0}表示矩阵，矩阵宽2，高1
1 0
B={2,0,3,1,1,0}，矩阵宽3，高2
2 0 3
1 1 0
C={2,0,3}，矩阵宽3，高1
2 0 3
ABC均按行存储
我们设置AB属性为CUBLAS_OP_N，cublas将AB按照列优先解析，并将A按宽1高2，B按宽2高3，C按宽1高3的方式设置。
则cublas认为B是
2 1
0 1
3 0
即B的转置，AC同理
那么计算C=A*B时，我们按照上面设置将参数传入cublasSgemm函数，得到的即为C'=B'*A'即实现了行优先存储乘法（行优先和列优先在计算时，输入输出矩阵数据分别按照行优先和列优先给出，赋值时宽高不必做转换，该是什么就是什么，但
cublasSgemm函数中的宽高需要改成其转置，A和B也要互换，具体看代码）
行优先存储全部代码如下

/*
 * cublasLearn.cpp
 *
 *  Created on: 2017年4月6日
 *      Author: wjt
 */

// CUDA runtime 库 + CUBLAS 库
#include "cuda_runtime.h"
#include "cuda_runtime_api.h"
#include "cublas_v2.h"

#include <time.h>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "MyTimer.h"
//本代码实现了cublas的行存储矩阵乘法
//#include <common_functions.h>
/*
|1 0| * |2 0 3| = |2 0 3|
        |1 1 0|
A     * B       = C
cublasSgemm 实现的功能为C = alpha*A*B + beta*C
ABC均为列存储
*/
using namespace std;

// 定义测试矩阵的维度
int const aW = 2;
int const aH = 1;
int const bW = 3;
int const bH = 2;
int const cW = bW;
int const cH = aH;

int main()
{
    // 定义状态变量
    cublasStatus_t status;

    // 在 内存 中为将要计算的矩阵开辟空间
    float *h_A = (float*)malloc (aW*aH*sizeof(float));
    float *h_B = (float*)malloc (bW*bH*sizeof(float));

    // 在 内存 中为将要存放运算结果的矩阵开辟空间
    float *h_C = (float*)malloc (cW*cH*sizeof(float));

    // 为待运算矩阵的元素赋值
    h_A[0]=1; h_A[1]=0;
    h_B[0]=2; h_B[1]=0; h_B[2]=3; h_B[3]=1; h_B[4]=1; h_B[5]=0;

    // 打印待测试的矩阵
    cout << "矩阵 A :" << endl;
    for (int i=0; i<aW*aH; i++){
        cout << h_A[i] << " ";
        if ((i+1)%aW == 0) cout << endl;
    }
    cout << endl;
    cout << "矩阵 B :" << endl;
    for (int i=0; i<bW*bH; i++){
        cout << h_B[i] << " ";
        if ((i+1)%bW == 0) cout << endl;
    }
    cout << endl;

    /*
    ** GPU 计算矩阵相乘
    */

    // 创建并初始化 CUBLAS 库对象
    cublasHandle_t handle;
    status = cublasCreate(&handle);

    if (status != CUBLAS_STATUS_SUCCESS)
    {
        if (status == CUBLAS_STATUS_NOT_INITIALIZED) {
            cout << "CUBLAS 对象实例化出错" << endl;
        }
        getchar ();
        return EXIT_FAILURE;
    }

    float *d_A, *d_B, *d_C;
    // 在 显存 中为将要计算的矩阵开辟空间
    cudaMalloc (
        (void**)&d_A,    // 指向开辟的空间的指针
        aW*aH * sizeof(float)    //　需要开辟空间的字节数
    );
    cudaMalloc (
        (void**)&d_B,
        bW*bH * sizeof(float)
    );

    // 在 显存 中为将要存放运算结果的矩阵开辟空间
    cudaMalloc (
        (void**)&d_C,
        cW*cH * sizeof(float)
    );

    // 将矩阵数据传递进 显存 中已经开辟好了的空间
    cublasSetVector (
        aW*aH,    // 要存入显存的元素个数
        sizeof(float),    // 每个元素大小
        h_A,    // 主机端起始地址
        1,    // 连续元素之间的存储间隔
        d_A,    // GPU 端起始地址
        1    // 连续元素之间的存储间隔
    );
    cublasSetVector (
        bW*bH,
        sizeof(float),
        h_B,
        1,
        d_B,
        1
    );

    // 同步函数
    cudaThreadSynchronize();

    Timer myTimer;
    myTimer.start();
    // 传递进矩阵相乘函数中的参数，具体含义请参考函数手册。
    float a=1; float b=0;
    // 矩阵相乘。该函数必然将数组解析成列优先数组
    cublasSgemm (
        handle,    // blas 库对象
        CUBLAS_OP_N,    // 矩阵 A 属性参数
        CUBLAS_OP_N,    // 矩阵 B 属性参数
        bW,    // A, C 的行数
        aH,    // B, C 的列数
        bH,    // A 的列数和 B 的行数
        &a,    // 运算式的 α 值
        d_B,    // A 在显存中的地址
        bW,    // lda
        d_A,    // B 在显存中的地址
        aW,    // ldb
        &b,    // 运算式的 β 值
        d_C,    // C 在显存中的地址(结果矩阵)
        cW    // ldc
    );
    // 同步函数
    cudaThreadSynchronize();
    myTimer.stop();
    printf("use: %lf\n", myTimer.getElapsedTime());

    // 从 显存 中取出运算结果至 内存中去
    cublasGetVector (
        cW*cH,    //  要取出元素的个数
        sizeof(float),    // 每个元素大小
        d_C,    // GPU 端起始地址
        1,    // 连续元素之间的存储间隔
        h_C,    // 主机端起始地址
        1    // 连续元素之间的存储间隔
    );

    // 打印运算结果
    cout << "计算结果的转置 ( (A*B)的转置 )：" << endl;

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

    // 清理掉使用过的内存
    free (h_A);
    free (h_B);
    free (h_C);
    cudaFree (d_A);
    cudaFree (d_B);
    cudaFree (d_C);

    // 释放 CUBLAS 库对象
    cublasDestroy (handle);

    getchar();

    return 0;
}