BLAS 库如何直接与 boost 多阵列一起使用？答案

【问题标题】：How can BLAS library be directly used with boost multiarrays?BLAS 库如何直接与 boost 多阵列一起使用？
【发布时间】：2020-12-11 14:41:57
【问题描述】：

给定两个矩阵A 和B，C 是它们相乘的结果。

#include "boost/multi_array.hpp"
typedef boost::multi_array<double, 2> matrix;
int m=5;
int n=6;
int k=7;

matrix A(boost::extents[m][k]);
matrix B(boost::extents[k][n]);
matrix C(boost::extents[m][n]);

如何调用blas 库中的dgemm 函数来计算A 和B 的矩阵乘积？我知道uBLAS 部分boost 库、armadillo、MTL 4、eigen 和其他一些为blas 函数提供方便包装的库。这里的问题是如何直接在多数组上调用dgemm。

【问题讨论】：

标签： c++11 boost matrix-multiplication

【解决方案1】：

您可以访问连续的元素存储。

原型是

void cblas_dgemm( CBLAS_LAYOUT layout, CBLAS_TRANSPOSE TransA, CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc )

那么，让我们填写吧：

cblas_dgemm(
    CBLAS_LAYOUT::CblasRowMajor,
    CBLAS_TRANSPOSE::CblasNoTrans,
    CBLAS_TRANSPOSE::CblasNoTrans,
    m, n, k,
    1.0, // alpha
    A.data(), A.shape()[1],
    B.data(), B.shape()[1],
    0.0, // beta
    C.data(), C.shape()[1]);

那是使用docs here 获得一些关于 alpha/beta 的指导：

演示

#include <boost/multi_array.hpp>
typedef boost::multi_array<double, 2> matrix;

#include <iostream>
namespace io { // for debug output
    auto& dump(std::ostream& os, double v) { return os << v; }

    template <typename R> auto& dump(std::ostream& os, R const& r) {
        std::string_view sep = "";
        os << "{";
        for (auto const& el : r) { dump(os << sep, el); sep = ","; }
        return os << "}";
    }
}

std::ostream& operator<<(std::ostream& os, matrix const& m) { return io::dump(os, m); }

// demo
#include <cblas-netlib.h>
#include <numeric> // iota

int main() {
    constexpr auto m=5, n=6, k=7;

    matrix A(boost::extents[m][k]);
    matrix B(boost::extents[k][n]);
    matrix C(boost::extents[m][n]);

    std::iota(A.data(), A.data() + A.num_elements(), 0);
    std::iota(B.data(), B.data() + B.num_elements(), 50);
    std::iota(C.data(), C.data() + C.num_elements(), 100);

    std::cout << "A: " << A << "\nB: " << B << "\n";
    assert(A.storage_order().all_dims_ascending());

    /*
     * void cblas_dgemm(
     *   CBLAS_LAYOUT layout,
     *   CBLAS_TRANSPOSE TransA,
     *   CBLAS_TRANSPOSE TransB,
     *   const int M, const int N, const int K,
     *   const double alpha,
     *   const double *A, const int lda,
     *   const double *B, const int ldb,
     *   const double beta,
     *   double *C, const int ldc )
     */

    cblas_dgemm(
        CBLAS_LAYOUT::CblasRowMajor,
        CBLAS_TRANSPOSE::CblasNoTrans,
        CBLAS_TRANSPOSE::CblasNoTrans,
        m, n, k,
        1.0, // alpha
        A.data(), A.shape()[1],
        B.data(), B.shape()[1],
        0.0, // beta
        C.data(), C.shape()[1]);

    std::cout << "C:\n" << C << "\n";
}

哪些打印：

A: {{0,1,2,3,4,5,6},{7,8,9,10,11,12,13},{14,15,16,17,18,19,20},{21,22,23,24,25,26,27},{28,29,30,31,32,33,34}}
B: {{50,51,52,53,54,55},{56,57,58,59,60,61},{62,63,64,65,66,67},{68,69,70,71,72,73},{74,75,76,77,78,79},{80,81,82,83,84,85},{86,87,88,89,90,91}}
C:
{{1596,1617,1638,1659,1680,1701},{4928,4998,5068,5138,5208,5278},{8260,8379,8498,8617,8736,8855},{11592,11760,11928,12096,12264,12432},{14924,15141,15358,15575,15792,16009}}

这与Wolfram Alpha 结帐：

【讨论】：

太棒了，这完全解决了我的问题，希望对其他人也有用！