贴一个cublas库的矩阵乘法封装

最新推荐文章于 2023-03-15 11:17:41 发布
最新推荐文章于 2023-03-15 11:17:41 发布
阅读量161 收藏
点赞数
分类专栏： 笔记 文章标签： 线性代数 算法 c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/m0_54206076/article/details/120373815
版权

底层直接用vector模板，对于double和float可以直接开cublas，也可以用普通的三层循环，小规模矩阵乘法，向量乘法都可以。当然，功能可以说是比较简陋，够用就成

对于cublas的调优，基本就是他了。大矩阵能用，小矩阵灵活，80008000的不惧，44的也可。

#ifndef MATRIX
#define MATRIX

#include <cuda/cuda_runtime.h>
#include <cuda/cublas_v2.h>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <typeinfo>

namespace Lnr
{
    template <typename T>
    struct Vec;

    template <typename T>
    struct Mat
    {
        Mat() = delete;

        Mat(unsigned r = 1, unsigned c = 1, T d = 0)
            : mat(std::vector<T>(r * c, d)), row(r), col(c) {}

        Mat(const std::vector<T> &rhs, unsigned r, unsigned c)
            : mat(rhs), row(r), col(c) {}

        Mat(std::vector<T> &&rhs, unsigned r, unsigned c)
            : mat(rhs), row(r), col(c) {}

        Mat(const Mat &rhs)
            : mat(rhs.mat), row(rhs.row), col(rhs.col) {}

        Mat(Mat &&rhs)
            : mat(rhs.mat), row(rhs.row), col(rhs.col) {}

        Mat &operator=(const Mat &rhs)
        {
            mat = rhs.mat;
            row = rhs.row;
            col = rhs.col;
            return *this;
        }

        Mat &operator=(Mat &&rhs)
        {
            if (this != &rhs)
            {
                mat = std::move(rhs.mat);
                row = rhs.row;
                col = rhs.col;
            }
            return *this;
        }

        ~Mat() {}

        const T *data() const
        {
            return mat.data();
        }

        T *data()
        {
            return mat.data();
        }

        const std::vector<T> getdata() const { return mat; }

        //矩阵转置
        void transpose()
        {
            if (row * col != 0)
            {
                Mat<T> result(col, row);
#pragma omp parallel for num_threads(8) schedule(dynamic)
                for (size_t i = 0; i != row; ++i)
                {
                    for (size_t j = 0; j != col; ++j)
                    {
                        std::swap(result(j, i), (*this)(i, j));
                    }
                }
                *this = std::move(result);
            }
        }

        inline const T &operator()(unsigned r, unsigned c) const
        {
            if (r < row && c < col)
            {
                return mat[r * col + c];
            }
            else
            {
                std::cout << "  out of row or col at " << __LINE__ << std::endl;
            }
            return mat[0];
        }

        inline T &operator()(unsigned r, unsigned c)
        {
            if (r < row && c < col)
            {
                return mat[r * col + c];
            }
            else
            {
                std::cout << "  out of row or col  at " << __LINE__ << std::endl;
            }
            return mat[0];
        }

        //矩阵向量乘积
        Vec<T> operator*(const Vec<T> &rhs) const
        {
            if (col == rhs.row)
            {
                Vec<T> result(row);

#pragma omp parallel for num_threads(8) schedule(dynamic)
                for (size_t i = 0; i != row; ++i)
                {
                    for (size_t j = 0; j != col; ++j)
                    {
                        result(i) += (*this)(i, j) * rhs(j);
                    }
                }

                return result;
            }
            return Vec<T>(1);
        }

        //矩阵相乘
        Mat<T> operator*(const Mat<T> &rhs) const
        {
            if (col == rhs.row)
            {
                Mat<T> result(row, rhs.col);
                Mat<T> buff(rhs);
                buff.transpose();

#pragma omp parallel for num_threads(8) schedule(dynamic)
                for (size_t i = 0; i != row; ++i)
                {
                    for (size_t j = 0; j != buff.row; ++j)
                    {
                        for (size_t k = 0; k != col; ++k)
                        {
                            result(i, j) += (*this)(i, k) * buff(j, k);
                        }
                    }
                }

                return result;
            }

            return Mat<T>(1, 1);
        }

        //CUDA初始化
        void cuInit(int &devID)
        {
            devID = 0;
            cudaError_t error;
            cudaDeviceProp deviceProp;

            error = cudaSetDevice(devID);
            err(error);

            error = cudaGetDevice(&devID);
            err(error);

            error = cudaGetDeviceProperties(&deviceProp, devID);
            err(error);
        }

        //cublas矩阵相乘
        Mat<T> cuMult(const Mat<T> &rhs, const unsigned &rowLmt) const
        {
            std::vector<T> result;
            result.reserve(row * rhs.col);

            cudaError_t error;
            cublasStatus_t status;
            cublasHandle_t handle;
            cudaDeviceProp deviceProp;
            int devID = 0;
            error = cudaGetDeviceProperties(&deviceProp, devID);
            err(error);

            int block_size = (deviceProp.major < 2) ? 8 : 16; //原来为16：32，过8000崩溃，设置小一倍，增加稳定性。

            dim3 threads(block_size, block_size);
            dim3 grid(rhs.col / threads.x, (row > rowLmt ? rowLmt : row) / threads.y);
            cublasCreate(&handle);

            unsigned int mem_size_A = sizeof(T) * col * (row > rowLmt ? rowLmt : row);
            unsigned int mem_size_B = sizeof(T) * rhs.col * rhs.row;
            unsigned int mem_size_C = sizeof(T) * rhs.col * (row > rowLmt ? rowLmt : row);

            T *d_A, *d_B, *d_C;

            error = cudaMalloc((void **)&d_A, mem_size_A);
            err(error);

            error = cudaMalloc((void **)&d_B, mem_size_B);
            err(error);

            error = cudaMalloc((void **)&d_C, mem_size_C);
            err(error);

            T *h_CUBLAS = (T *)malloc(mem_size_C);

            const T *h_B = rhs.data();

            //B不变所以提到循环之前。
            error = cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice);
            err(error);

            int loopNum = 0;
            for (int i = row; i > 0; i -= rowLmt)
            {
                int rowLoop = i > rowLmt ? rowLmt : i;

                const T *h_A = (*this).data() + loopNum * rowLmt * col; //易出错地方，每循环一次，
                                                                        //指针移动 “循环次数 * 限定行数 * 左矩阵列数。

                error = cudaMemcpy(d_A, h_A, col * rowLoop * sizeof(T), cudaMemcpyHostToDevice);
                err(error);

                {
                    cublasCore(handle, error, status, rhs.col, rowLoop, d_B, d_A, d_C, T(0));

                    error = cudaMemcpy(h_CUBLAS, d_C, mem_size_C, cudaMemcpyDeviceToHost);
                    err(error);

                    result.insert(result.end(), h_CUBLAS, h_CUBLAS + rowLoop * rhs.col); //输出尺寸：计算行 X 后矩阵列数
                }

                ++loopNum;
            }

            {
                status = cublasDestroy_v2(handle);
                stat(status);

                error = cudaFree(d_A);
                err(error);

                error = cudaFree(d_B);
                err(error);

                error = cudaFree(d_C);
                err(error);

                free(h_CUBLAS);
            }
            return Mat<T>(std::move(result), row, rhs.col);
        }

    private:
        friend struct Vec<T>;

        void err(cudaError_t &error) const
        {
            if (error != cudaSuccess)
            {
                printf("cuda returned error code %d, line(%d)\n", error, __LINE__);
                exit(EXIT_FAILURE);
            }
        }

        void stat(cublasStatus_t &status) const
        {
            if (status != CUBLAS_STATUS_SUCCESS)
            {
                printf(" cublas returned error code %d, line(%d)\n", status, __LINE__);
                exit(EXIT_FAILURE);
            }
        }

        void cublasCore(cublasHandle_t &handle, cudaError_t &error, cublasStatus_t &status, unsigned rhscol, unsigned rowLoop, T *d_B, T *d_A, T *d_C, float f) const
        {
            const T alpha = 1.0;
            const T beta = 0.0;
            status = cublasSgemm_v2(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                                    rhscol, rowLoop, col,
                                    &alpha,
                                    d_B, rhscol,
                                    d_A, col,
                                    &beta,
                                    d_C, rhscol);
            stat(status);
            error = cudaDeviceSynchronize();
            err(error);
        }

        void cublasCore(cublasHandle_t &handle, cudaError_t &error, cublasStatus_t &status, unsigned rhscol, unsigned rowLoop, T *d_B, T *d_A, T *d_C, double d) const
        {
            const T alpha = 1.0;
            const T beta = 0.0;
            status = cublasDgemm_v2(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                                    rhscol, rowLoop, col,
                                    &alpha,
                                    d_B, rhscol,
                                    d_A, col,
                                    &beta,
                                    d_C, rhscol);
            stat(status);
            error = cudaDeviceSynchronize();
            err(error);
        }

        std::vector<T> mat;
        unsigned row;
        unsigned col;
    };

    template <typename T>
    struct Vec
    {
        Vec() = delete;

        Vec(unsigned r = 1, T d = 0)
            : vec(std::vector<T>(r, d)), row(r) {}

        Vec(const Vec &rhs)
            : vec(rhs.vec), row(rhs.row) {}

        Vec(Vec &&rhs)
            : vec(rhs.vec), row(rhs.row) {}

        Vec &operator=(const Vec &rhs)
        {
            vec = rhs.vec;
            row = rhs.row;
            return *this;
        }

        Vec &operator=(Vec &&rhs)
        {
            if (this != &rhs)
            {
                vec = std::move(rhs.vec);
                row = rhs.row;
            }
            return *this;
        }

        ~Vec() {}

        //向量取值
        T &operator()(unsigned num)
        {
            if (num < row)
            {
                return vec[num];
            }
        }

        //const向量取值
        const T &operator()(unsigned num) const
        {
            if (num < row)
            {
                return vec[num];
            }
        }

        //向量取首元素指针
        inline T *data()
        {
            return vec.data();
        }

        //向量取首元素指针
        inline const T *data() const
        {
            return vec.data();
        }

        //向量点乘
        T operator*(const Vec &x) const
        {
            T matdot = 0;
            if (row == x.row)
            {
                for (size_t i = 0; i != row; ++i)
                {
                    matdot += x(i) * (*this)(i);
                }
            }
            return matdot;
        }

        //向量矩阵乘积
        Vec<T> operator=(const Mat<T> &rhs) const
        {
            if (row = rhs.row)
            {
                Mat<T> result(rhs);
                result.transpose();
                return result * (*this);
            }
        }

    private:
        friend struct Mat<T>;

        std::vector<T> vec;
        unsigned row;
    };

} // namespace Lnr

#endif
Keras中 用Lambda层封装的矩阵乘法
sxh371302的博客
01-05 2099
Keras中 用Lambda层封装的矩阵乘法 keras中对于层的封装绝对是一个大坑，今天分享一些自己在写代码中遇到的bug。如果是对层进行矩阵乘法，假设输入层分别为A:（batch_size, 300, 64),B:(batch_size,300,64)，我们需要对A进行转置再与B进行矩阵乘法，但是Keras中无法直接对tensor进行矩阵乘法，会报错，必须要对其进行封装。 # A:(batch_size, 300, 64),B:(batch_size,300,64) A = Permute((2, 1)
cuda矩阵相乘_CUDA / CUBLAS矩阵向量乘法
weixin_39629947的博客
12-20 401
I previously posted a question regarding matrix-vector multiplication in CUDA and about writing my own kernel. After doing this, I decided to implement my problem using CUBLAS as suggested by some use...
{模板}打包矩阵乘法
FarmerJohn
05-26 415
Code...
解矩阵方程+矩阵乘法，已封装，可直接调用
03-12
解矩阵方程+矩阵乘法，已封装，可直接调用， // private static double[][] matrixC ; /** * * 解矩阵方程 * * @param A * 矩阵 * @param b * 列阵 * @return matrixC 列阵 */ /** * * 矩阵相乘 * * @param A * 矩阵 * @param b 列阵 * * @return c 列阵 */ public static double[][] MatrixMultiply(double[][] A, double[][] b) {
【矩阵运算c++实现】矩阵封装实现Matrix类
安安爸Chris的专栏
11-24 2161
Matrix类封装了矩阵运算里的常用几种函数 Matrix.h #include <string> #include <sstream> typedef struct MatrixShape { int row; int col; int size() const { return row*col; } std::string ToStr() const { std::stringstream ss;.
CUDA 10.2版本CUBLAS库文档：GPU线性代数运算指南
对于使用CUDA进行编程的初学者而言，该文档是一个非常有价值的入门资源，它不仅介绍了cuBLAS库的基本使用，还包含了性能调优和故障排除的指导，从而帮助开发者更有效地利用GPU的计算能力。 cuBLAS库中的函数可以...
python 矩阵乘法 点乘_numpy中矩阵乘法，星乘(*)和点乘(.dot)的区别
weixin_42388898的博客
12-23 3411
import numpya = numpy.array([[,],[,]])b = numpy.array([[,],[,]])星乘表示矩阵内各对应位置相乘，矩阵a*b下标(0,0)=矩阵a下标(0,0) x 矩阵b下标(0,0)；点乘表示求矩阵内积，二维数组称为矩阵积(mastrixproduct)。数学上的概念 不一样1、乘积用于矩阵相乘，表示为C=A*B，A的列数与B的行数必须相同，C也是...
CUDA Samples 之 matrixMul 矩阵乘法详解
qq_43232373的博客
10-08 1026
CUDA实现矩阵乘法详解
C++矩阵乘法计算 || GPU && CPU 实现
春江水暖鸭先知
09-03 5241
废话 使用英伟达加速矩阵运算一般使用cublas库进行操作，在cublas库中，使用cublasSgemv()和cublasSgemm()分别进行矩阵向量间的乘法运算与矩阵矩阵间的乘法运算。本文将具体的解释上述两个函数的参数以及具体的应用例子。参照官方解释 矩阵乘法函数解释 仅仅使用单精度(float)进行解释，若想尝试双精度(double),可以进入官网查看相应函数. ...
cublas64-10.rar
最新发布
09-09
本文将深入探讨cublas64_10.dll，它是CUDA BLAS库的一个重要组成部分。 cublas64_10.dll是64位架构下的CUDA BLAS库的一个动态链接库文件，版本号为10。在Windows操作系统中，动态链接库文件（DLL）用于封装特定的...
矩阵乘法及矩阵快速幂实现及思想
qq_44654498的博客
05-27 343
记录矩阵乘法及矩阵快速幂的实现
Cuda 编程 矩阵转置
ECNU_sei_mfj的博客
07-10 936
将矩阵分块进行转置 #include <stdio.h> #include<stdlib.h> #include <cuda_runtime.h> #include<device_launch_parameters.h> #include<device_functions.h> #define imin(a,b) (a<b?a:b) //const int N=33*1024; const int N = 16; const int
cublas库实现矩阵乘法（任意维数）
qq632544991p的专栏
11-17 2477
话不多说，直接进入主题 cublas是CUDA上矩阵运算的库，可以在gpu上实现很高的效率。然而关于它的使用，并没有详细的中文资料 笔者，经过多次调试成功，分享一点儿心得 #include #include #include // Multiply the arrays A and B on GPU and save the result in C // C(m,n) = A(
cuda矩阵转置
lumanman_的博客
09-05 1691
一直以为经典的cuda矩阵转置只能用于矩阵的宽高都能被线程块大小整除的情况，也是很奇怪，不知道怎么形成这个概念的，然后这次又要用到，本想着大干一番，把宽高不能被线程块整除的矩阵转置攻克了，可是没想到一测试，人本来就可以实现，这就尴尬了, 所以在此记录下来，纠正自己的这个错误。 代码： #define BLOCK_DIM 16__global__ void myTranspose
cuBLAS使用2
qq_44632658的博客
03-15 2428
无
cublas中的矩阵运算
热门推荐
l691899397的博客
08-17 1万+
Cublas是一个可以与cuda一同使用的函数库，它提供了多种矩阵运算的API，但是它列主序的存储方式却让人十分疑惑，今天我就以cublas中的矩阵乘法运算简单说一下我的理解。 Cublas中的矩阵乘法运算函数有5个，分别是cublasSgemm、cublasDgemm、cublasCgemm、cublasZgemm、cublasHgemm，分别包括了不同数据类型的计算，比如单精度浮点、双精度浮
基于深度学习的人脸识别系统系列——使用CUBLAS加速计算人脸向量的余弦距离
Allyli0022的专栏
07-04 3363
前言 基于深度学习的人脸识别系统，一共用到了5个开源库：OpenCV（计算机视觉库）、Caffe（深度学习库）、Dlib（机器学习库）、libfacedetection（人脸检测库）、cudnn（gpu加速库）。  用到了一个开源的深度学习模型：VGG model。  最终的效果是很赞的，识别一张人脸的速度是0.039秒，而且最重要的是：精度高啊！！！  CPU：intel i5-45
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

不停感叹的老林_<C 语言编程核心突破>

不打赏的人, 看完也学不会.

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值