c++矩陣計算效能對比:Eigen和GPU解讀

2022-12-17 14:00:51

生成隨機矩陣

生成隨機矩陣有多種方式，直接了當的方式是使用顯式迴圈的方式為矩陣的每個元素賦隨機值。

#include <iostream>
#include <random>

using namespace std;

// 生成亂數
double GenerateRandomRealValue()
{
    std::random_device rd;
    std::default_random_engine eng(rd());
    std::uniform_real_distribution<double> distr(1, 10);
    return distr(eng);
}

int main()
{
        // 3d矩陣
    double a[3][3];
    for (int i = 0; i < 3; ++i) {
        for (int j = 0;  j < 3; ++j) {
            a[i][j] = GenerateRandomRealValue();
        }
    }

    return 0;
}

另一種方式是使用Eigen庫，它提供了矩陣運算的庫。

生成隨機矩陣：

#include "Eigen/Dense"
#include <functional>

using namespace std;
using namespace Eigen;

MatrixXd Generate2DMatrixByEigen()
{
        // 直接使用內建的Random，產生均勻分佈隨機矩陣
    MatrixXd m = MatrixXd::Random(3,3);
    
    // 也可以呼叫自定義的亂數生成函數填充資料
    // MatrixXd m = MatrixXd::Zero(3,3).unaryExpr(std::bind(GenerateRandomRealValue));
    return m;
}

計算矩陣點積

使用顯式迴圈計算

直接上程式碼：

void CalcMatrixDotForLoop(const vector<vector<double>>& a, const vector<vector<double>>& b)
{
    std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
    if (a[0].size() != b.size()) {
        cout << "error:" << a.size() << "," << b[0].size() << endl;
        return;
    }

    vector<vector<double>> c;
    vector<double> c_row(b[0].size());
    for (int i = 0; i < a.size(); ++i) {
        for (int j = 0; j < b[0].size(); ++j) {
            for (int k = 0; k < b.size(); ++k) {
                c_row[j] += a[i][k] * b[k][j];
            }
        }
        c.emplace_back(c_row);
    }
    std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double, std::milli> time_span = t2 - t1;
    std::cout << "Loop takes " << time_span.count() << " msn";

    // cout << "matrix c:n";
    // for (int i = 0; i < c.size(); ++i) {
    //     for (int j = 0; j < c[0].size(); ++j) {
    //         cout << c[i][j] << ",";
    //     }
    //     cout << endl;
    // }
}

使用Eigen庫

程式碼：

void ModeEigen(const int a_row, const int a_col, const int b_row, const int b_col)
{
    std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
    auto c = a * b;
    std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
    std::chrono::duration<double, std::milli> time_span = t2 - t1;
    std::cout << "Eigen takes " << time_span.count() << " msn";
    // cout << "matrix c:n" << c << endl;
}

使用GPU

程式碼片斷：

auto t_begin = std::chrono::high_resolution_clock::now();

t1 = std::chrono::high_resolution_clock::now();
cudaMalloc((void**)&da,size);
cudaMalloc((void**)&db,size);
cudaMalloc((void**)&dc,size);
t2 = std::chrono::high_resolution_clock::now();
time_span = t2 - t1;
std::cout << "GPU malloc takes " << time_span.count() << " msn";

t1 = std::chrono::high_resolution_clock::now();
cudaMemcpy(da,a,size,cudaMemcpyHostToDevice);
cudaMemcpy(db,b,size,cudaMemcpyHostToDevice);
t2 = std::chrono::high_resolution_clock::now();
time_span = t2 - t1;
std::cout << "cudaMemcpy takes " << time_span.count() << " msn";

t1 = std::chrono::high_resolution_clock::now();
dim3 dg(32,32);
dim3 dbs((n+dg.x-1)/dg.x,(n+dg.y-1)/dg.y);
mextix<<<dbs,dg>>>(da,db,dc,n);
t2 = std::chrono::high_resolution_clock::now();
time_span = t2 - t1;
std::cout << "gpu takes " << time_span.count() << " msn";

t1 = std::chrono::high_resolution_clock::now();
cudaMemcpy(c,dc,size,cudaMemcpyDeviceToHost);
t2 = std::chrono::high_resolution_clock::now();
time_span = t2 - t1;
std::cout << "cudaMemcpy back takes " << time_span.count() << " msn";

cudaFree(da);
cudaFree(db);
cudaFree(dc);

auto t_end = std::chrono::high_resolution_clock::now();
time_span = t_end - t_begin;
std::cout << "GPU total takes " << time_span.count() << " msn";