c++ - BLAS 产品 dgemm 使用 CblasTrans 选项时行为异常

Question

我想问你一个相当初级的 BLAS 问题。 看似简单的任务是关于矩阵 A 与它自己的矩阵乘法 转置:C := A'*A

我的例子是 (2x3): A:=[1 2 3 ; 4 5 6]。 因此 A' 是 (3x2)，C 应该是 (3x3)。

在 Row Major 中并计划使用我期望的 CblasTrans 选项 在 A 和 A' 两种情况下，lda=ldb=3。

可悲的是，较低的演示程序仍然生成了一个完全错误的产品 到目前为止，我的简单参数排列还没有达到目标。 事实上，结果值高得离谱，我 对结果的 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;
}

Answer 1

看看来自 netlib 的 DGEMM:http://www.netlib.org/blas/dgemm.f

你会看到:

*  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' 是 (3,2)。

如果查看其他参数的定义，您会发现必须传递 M=3、N=3、K=2