BLAS中的元素向量向量乘法？答案

【问题标题】：Element-wise vector-vector multiplication in BLAS?BLAS中的元素向量向量乘法？
【发布时间】：2011-10-01 16:32:17
【问题描述】：

有没有办法使用 BLAS、GSL 或任何其他高性能库进行逐元素向量乘法？

【问题讨论】：

标签： c++ blas

【解决方案1】：

（按字面意思理解问题的标题......）

是的，它可以单独使用 BLAS 完成（尽管它可能不是最有效的方式。）

诀窍是将输入向量之一视为对角矩阵：

⎡a    ⎤ ⎡x⎤    ⎡ax⎤
⎢  b  ⎥ ⎢y⎥ =  ⎢by⎥
⎣    c⎦ ⎣z⎦    ⎣cz⎦

然后您可以使用其中一个矩阵向量乘法函数，该函数可以将对角矩阵作为输入而无需填充，例如SBMV

例子：

void ebeMultiply(const int n, const double *a, const double *x, double *y)
{
    extern void dsbmv_(const char *uplo,
                       const int *n,
                       const int *k,
                       const double *alpha,
                       const double *a,
                       const int *lda,
                       const double *x,
                       const int *incx,
                       const double *beta,
                       double *y,
                       const int *incy);

    static const int k = 0; // Just the diagonal; 0 super-diagonal bands
    static const double alpha = 1.0;
    static const int lda = 1;
    static const int incx = 1;
    static const double beta = 0.0;
    static const int incy = 1;

    dsbmv_("L", &n, &k, &alpha, a, &lda, x, &incx, &beta, y, &incy);
}

// Test
#define N 3
static const double a[N] = {1,3,5};
static const double b[N] = {1,10,100};
static double c[N];

int main(int argc, char **argv)
{
    ebeMultiply(N, a, b, c);
    printf("Result: [%f %f %f]\n", c[0], c[1], c[2]);
    return 0;
}

Result: [1.000000 30.000000 500.000000]

【讨论】：

我知道这已经很晚了，但我只想说，虽然 finnw 的这个答案是有效的，但我可能不建议使用它。在我的实际案例中，自己编写循环要快得多（2-3 次）。我不知道我的编译器优化了多少，但通常切换到 blas 会产生很好的加速（例如，在另一个方向上提高 2-3 倍）而不是放慢速度。当然，这取决于几个因素，但只是作为计算时间的警告。
我可以确认@oli 的发现。对我来说，N = 300 的 ssbmv 比两个嵌套的 for 循环慢约 20 倍。我正在使用带有 -O2 -fPIC -fstack-protector-strong 选项的 Intel Xeon X7560、OpenBLAS 和 GCC 8.3.0。我的猜测是 ?sbmv 太笼统了，无法充分利用矢量化指令。

【解决方案2】：

我发现 MKL 在它的向量数学函数库 (VML) 中有一整套向量的数学运算，包括 v?Mul，它可以满足我的需求。它适用于 c++ 数组，所以它对我来说比 GSL 更方便。

【讨论】：

【解决方案3】：

总是有 std::valarray¹，它定义了在目标支持的情况下经常编译成 SIMD 指令的元素操作（英特尔 C++ /Quse-intel-optimized-headers，G++）。

http://software.intel.com/sites/products/documentation/hpc/composerxe/en-us/cpp/mac/cref_cls/common/cppref_valarray_intro.htm

这两个编译器也会做自动向量化

在这种情况下你可以写

#define N 10000 

float a[N], b[N], c[N]; 

void f1() { 
  for (int i = 1; i < N; i++) 
  c[i] = a[i] + b[i]; 
}

并看到它编译成矢量化代码（例如使用 SSE4）

¹ 是的，它们很陈旧，通常被认为已经过时，但实际上它们都是标准的，非常适合任务。

【讨论】：

你的第一个链接好像失效了。

【解决方案4】：

在 GSL 中，gsl_vector_mul 可以解决问题。

【讨论】：