涉及临时内存分配时避免 blas?

Avoid blas when involving temporary memory allocation?

我有一个程序可以重复计算矩阵乘积 x'Ay。更好的做法是通过调用 MKL 的 blas 来计算这个,即 cblas_dgemvcblas_ddot,这需要将内存分配给一个临时向量,或者最好简单地取 x_i * a_ij * y_j 的总和?换句话说,MKL的blas理论上有没有增加什么价值?

我为我的笔记本电脑进行了基准测试。除了 g++_no_blas 的性能是其他测试的两倍之外(为什么?),每个测试几乎没有区别。 O2、O3和Ofast之间也没有区别。

  1. g++_blas_static 57ms
  2. g++_blas_dynamic 58ms
  3. g++_no_blas 100 毫秒
  4. icpc_blas_static 57ms
  5. icpc_blas_dynamic58ms
  6. icpc_no_blas58ms

util.h

#ifndef UTIL_H
#define UTIL_H

#include <random>
#include <memory>
#include <iostream>

struct rng 
{
        rng() : unif(0.0, 1.0)
        {
        }

        std::default_random_engine re; 
        std::uniform_real_distribution<double> unif;

        double rand_double()
        {
                return unif(re);
        }

        std::unique_ptr<double[]> generate_square_matrix(const unsigned N)
        {
                std::unique_ptr<double[]> p (new double[N * N]);
                for (unsigned i = 0; i < N; ++i)
                {
                        for (unsigned j = 0; j < N; ++j)
                        {
                                p.get()[i*N + j] = rand_double();
                        }
                }
                return p;
        }

        std::unique_ptr<double[]> generate_vector(const unsigned N)
        {
                std::unique_ptr<double[]> p (new double[N]);
                for (unsigned i = 0; i < N; ++i)
                {
                        p.get()[i] = rand_double();
                }
                return p;
        }
};

#endif // UTIL_H

main.cpp

#include <iostream>
#include <iomanip>
#include <memory>
#include <chrono>
#include "util.h"
#include "mkl.h"

double vtmv_blas(double* x, double* A, double* y, const unsigned n)
{
        double temp[n];
        cblas_dgemv(CblasRowMajor, CblasNoTrans, n, n, 1.0, A, n, y, 1, 0.0, temp, 1); 
        return cblas_ddot(n, temp, 1, x, 1); 
}

double vtmv_non_blas(double* x, double* A, double* y, const unsigned n)
{
        double r = 0;
        for (unsigned i = 0; i < n; ++i)
        {
                for (unsigned j = 0; j < n; ++j)
                {
                        r += x[i] * A[i*n + j] * y[j];
                }
        }
        return r;
}

int main()
{
        std::cout << std::fixed;
        std::cout << std::setprecision(2);
        constexpr unsigned N = 10000;
        rng r;

        std::unique_ptr<double[]> A = r.generate_square_matrix(N);
        std::unique_ptr<double[]> x = r.generate_vector(N);
        std::unique_ptr<double[]> y = r.generate_vector(N);

        auto start = std::chrono::system_clock::now();
        const double prod = vtmv_blas(x.get(), A.get(), y.get(), N); 
        auto end = std::chrono::system_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(
                end - start);

        std::cout << "Result: " << prod << std::endl;
        std::cout << "Time (ms): " << duration.count() << std::endl;

GCC no blas 很差,因为它不使用矢量化 SMID 指令,而其他的都使用。 icpc 将自动矢量化你的循环。

您没有显示矩阵大小,但通常 gemv 受内存限制。由于矩阵比临时向量大得多,因此消除它可能无法大幅提高性能。