Blas产品DGEMM具有CBLASTRANS选项出乎意料的行为

BLAS product dgemm behaves unexpectedly with CblasTrans option

本文关键字：出乎意料 选项 CBLASTRANS 产品 DGEMM 具有 Blas      更新时间：2023-10-16

我想问你一个相当初学者的问题。看似简单的任务是关于矩阵a的矩阵乘法转座：c：= a'*a

我的示例是（2x3）：a：= [1 2 3;4 5 6]。因此，a'is（3x2）和c应为（3x3）。

在Row专业中，计划使用我期望的CBLASTRANS选项LDA = LDB = 3在两种情况下A和A'。

可悲的是，较低的演示程序仍然会产生完全错误的产品到目前为止，我的简单参数排列并未达到标记。事实上，由此产生的值很高，我是对结果的6元素结构感到困惑。

我在这里缺少什么？

/**
 * transposeMat.cpp, compile using: g++ -lcblas transposeMat.cpp
 */
#include <cstdlib>
#include <cblas.h>
#include <iostream>
#include <sstream>
#include <string>
using namespace std;
string matrix2string(int m, int n, double* A, CBLAS_ORDER order)
{
  ostringstream oss;
  for (int j=0;j<m;j++)
  {
    for (int k=0;k<n;k++)
    {
      switch (order)
      {
    case CblasRowMajor:
      oss << A[j*n+k];
      break;
    case CblasColMajor:
          oss << A[j+k*m];
      break;
    default:
      return "[matrix2string(..): Unknown order.]";
      }
      if (k < n-1) oss << 't';
    }
    if (j < m-1) oss << endl;
  }
  return oss.str();
}
int main(int argc, char** argv)
{
  int m=2;
  int n=3;
  // RowMajor matrix [ 1,2,3 ; 4,5,6 ]
  double A[6] = { 1,2,3,4,5,6 };
  // Using A for both xgemm-Parameters brings no luck! This is not enough though.
  double B[6] = { 1,2,3,4,5,6 }; 
  // Container for the result which will be 3x3.
  double C[9] = { 0,0,0,0,0,0,0,0,0 };
  // C:=A'A
  // Params: (Majority,TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,m,n,n,1,&A[0],n,&B[0],n,0,&C[0],n);
  //> ADDED COMMENT AFTER aka.nice ANSWERED THE QUESTION. ----------
  // 1.: "MxN" really are the dimensions of matrix C and K is the "in-between"
  //   dimension shared by the factors of the product.
  // 2.: The op(A) on the BLAS reference card actually seems to read "after
  //   the internal transpose of A".
  // 3.: Taken this into the code the above matrix B also becomes unnecessary.
  // Hence this programm runs expectedly if you
  //   replace the upper line by:
  // cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,n,n,m,1,&A[0],n,&A[0],n,0,&C[0],n);
  //< --------------------------------------------------------------
  cout << "A:" << endl << matrix2string(m,n,&A[0],CblasRowMajor).c_str() << endl <<
    "C:" << endl << matrix2string(n,n,&C[0],CblasRowMajor).c_str() << endl;
  /** Output:
  A:
  1       2       3
  4       5       6
  C:
  34      44      54
  90      117     144
  0       0       0
  */
  return EXIT_SUCCESS;
}

从Netlib：http：//www.netlib.org/blas/dgemm.m.f

看一下DGEMM

您会看到：

*  DGEMM  performs one of the matrix-matrix operations
*
*     C := alpha*op( A )*op( B ) + beta*C,

和：

*  M      - INTEGER.
*           On entry,  M  specifies  the number  of rows  of the  matrix
*           op( A )  and of the  matrix  C.  M  must  be at least  zero.
*           Unchanged on exit.

因此，如果a是（2,3），则OP（a）= a'is（3,2）。

如果您查看其他参数的定义，您将看到必须通过m = 3，n = 3，k = 2

以防有人想知道代码的外观不同：

 int m=2;
  int n=4;
  int k=3;
  // RowMajor matrix [ 1,2,3 ; 4,5,6 ]
  double A[6] = { 1,2,3,4,5,6 }; // A is 2x3
  // Using A for both xgemm-Parameters brings no luck! This is not enough though.
  double B[8] = { 1,2,3,4,5,6,7,8 };  // B 2x4
  // Container for the result which will be 3x4.
  double C[12] = { 0,0,0,0,0,0,0,0,0,0,0,0 }; // C is 3x4
cblas_dgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
    k,n,m,
    1.,&A[0],k,
    &B[0],n,0.,
    &C[0],n);
  cout << "A:" << endl << matrix2string(m,n,&A[0],CblasRowMajor).c_str() << endl <<
    "C:" << endl << matrix2string(k,n,&C[0],CblasRowMajor).c_str() << endl;
    std::exit(1);
相关文章：