打开Blas矩阵矩阵不快

如何解决打开Blas矩阵矩阵不快

我正在MacBook上玩Open BLAS。

我使用自己的实现和cblas_dgemm函数来实现基本矩阵矩阵乘法

我的Matrix.hpp文件如下所示：

template <class T>
class Matrix
{
public:
    int rows = -1;
    int cols = -1;

    T *values = nullptr;

    /* constructor */
    Matrix(int rows,int cols);

    /* set matrix */
    void setMatrix(T *&&new_values);

    /* multiples two matrices together,with BLAS optimization */
    Matrix<T> *matMatMult(Matrix<T> &rhs);

    /* multiples two matrices together,without BLAS optimization */
    Matrix<T> *matMatMult_unopt(Matrix<T> &rhs);

    /* prints the matrix */
    void print();

    /* destructor */
    ~Matrix();

protected:
    int size_of_values = -1;
};

值以行主格式存储。

我的Matrix.cpp文件如下所示：

#include "Matrix.hpp"
#include <iostream>
#include "cblas.h"

/*  
constructor
creates a matrix. All values initialized to 0.
*/
template <class T>
Matrix<T>::Matrix(int nrows,int ncols) : rows(nrows),cols(ncols),size_of_values(nrows * ncols)
{
    this->values = new T[this->size_of_values];

    for (int i = 0; i < size_of_values; i++)
        this->values[i] = 0;
}

/*  
destructor 
cleans up the Matrix instance.
*/
template <class T>
Matrix<T>::~Matrix()
{
    if (this->values != nullptr)
        delete[] this->values;
}

/*  
prints the matrix with rows and columns.
*/
template <class T>
void Matrix<T>::Matrix::print()
{
    for (int i = 0; i < this->rows; i++)
    {
        for (int j = 0; j < this->cols; j++)
            std::cout << this->values[i * this->cols + j] << " ";
        std::cout << std::endl;
    }
}

/* multiplies two matrices together. */
template <class T>
Matrix<T> *Matrix<T>::matMatMult(Matrix<T> &rhs)
{
    /* check dimensions make sense for matrix multiplication */
    if (this->rows != rhs.cols)
        throw std::invalid_argument("The right hand side matrix has incorrect column dimensions.");

    /* This holds the result */
    auto result = new Matrix<T>(this->rows,rhs.cols);

    cblas_dgemm(CblasRowMajor,CblasNoTrans,this->rows,rhs.cols,this->cols,1,this->values,rhs.values,result->values,rhs.cols);

    return result;
}

template <class T>
Matrix<T> *Matrix<T>::matMatMult_unopt(Matrix &rhs)
{
    /* check dimensions make sense return without doing any multiplication */
    if (this->cols != rhs.rows)
        throw std::invalid_argument("The right hand side matrix has incorrect column dimensions.");

    /* create an output matrix that will hold our values */
    auto result = new Matrix(this->rows,this->preallocated);

    /*======================================
     *  matrix multiplication is O(n^3).
     *  Although this is loop ordering takes advantage of caching,it
     *  does not take advantage of BLAS routines (for row by row access).
     *======================================*/
    for (int i = 0; i < this->rows; i++)
        for (int k = 0; k < this->cols; k++)
            for (int j = 0; j < rhs.cols; j++)
                result->values[i * result->cols + j] += this->values[i * this->cols + k] * rhs.values[k * rhs.cols + j];

    return result;
}

/* Transposes a matrix in place */
template <class T>
void Matrix<T>::Matrix::transpose()
{
    T *transpose_values = new T[this->size_of_values];

    for (int i = 0; i < this->rows; i++)
        for (int j = 0; j < this->cols; j++)
            transpose_values[i * this->cols + j] = this->values[j * this->cols + i];

    delete[] this->values;

    this->values = transpose_values;

    int temp;
    temp = this->rows;
    this->rows = this->cols;
    this->cols = temp;
}

/* sets the values inside a matrix. */
template <class T>
void Matrix<T>::setMatrix(T *&&new_values)
{
    auto tmp = new_values;
    new_values = nullptr;
    delete[] this->values;
    this->values = tmp;
}

这是我的计时测试：

TEST(UnitTests,SpeedTest)
{
    /* flag to check if we pass or fail the test */
    bool test = true;

    int rows = 10000;
    int cols = 10000;

    // the values we want to set in our matrices
    double *A_values = new double[rows * cols]{10};

    double *B_values = new double[rows * cols]{20};

    std::unique_ptr<Matrix<double>> A(new Matrix<double>(rows,cols));

    std::unique_ptr<Matrix<double>> B(new Matrix<double>(rows,cols));

    A->setMatrix(std::move(A_values));

    B->setMatrix(std::move(B_values));

    /* timing tests for the BLAS optimized matrix multiplication */
    auto start_blas = system_clock::now();

    A->matMatMult(*B);

    auto end_blas = system_clock::now() - start_blas;

    auto time_blas = duration<double>(end_blas).count();

    std::cout << "optimized,time in seconds: " << time_blas << std::endl;

    // the values we want to set in our matrices
    double *C_values = new double[rows * cols]{10};

    double *D_values = new double[rows * cols]{20};

    std::unique_ptr<Matrix<double>> C(new Matrix<double>(rows,cols));

    std::unique_ptr<Matrix<double>> D(new Matrix<double>(rows,cols));

    C->setMatrix(std::move(C_values));

    D->setMatrix(std::move(D_values));

    /* timing tests for the BLAS optimized matrix multiplication */
    auto start_blas2 = system_clock::now();

    C->matMatMult(*D);

    auto end_blas2 = system_clock::now() - start_blas2;

    auto time_blas2 = duration<double>(end_blas2).count();

    std::cout << "unoptimized,time in seconds: " << time_blas2 << std::endl;

    EXPECT_TRUE(test);
}

现在的问题是，使用BLAS的矩阵矩阵乘法与不使用BLAS的朴素矩阵在时间上没有显着差异。

有人可以解释为什么会这样吗，或者我做错了什么以免加速？