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

Question

我有一个反复计算矩阵乘积x'Ay的程序。通过调用MKL的blas(即cblas_dgemv和cblas_ddot)来计算这一点是否更好，这需要将内存分配给临时向量，或者更好地采用x_i * a_ij * y_j的总和?换句话说，从理论上讲，MKL的blas有任何价值吗?

我对我的笔记本电脑进行了基准测试。除了g++_no_blas的性能比其他测试差两倍(为什么?)之外，每个测试实际上没有什么区别。O2、O3和Ofast之间也无差异。

1. g++ _blas_static 57女士
2. g++ _blas_dynamic 58女士
3. g++ _no_blas 100 ms
4. icpc_blas_static 57女士
5. icpc_blas_dynamic 58女士
6. icpc_no_blas 58女士

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;

Answer 1

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

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