BLAS product dgemm behaves unexpectedly with CblasTrans option

Question

I would like to ask you a rather beginner BLAS question. The seemingly simple task is about a matrix multiplication of matrix A with its own transpose: C := A'*A

My example is (2x3): A:=[1 2 3 ; 4 5 6]. Hence A' is (3x2) and C should be (3x3).

In Row Major and planning to use the CblasTrans option I'd expect lda=ldb=3 in both cases A and A'.

Sadly the lower demo program still generates a totally wrong product and simple parameter permutation on my part did not hit the mark so far. As a matter of fact the resulting values are ridiculously high and I am baffled by the 6-element structure of the result.

What am I missing here?

/**
 * 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

Take a look at DGEMM from netlib: http://www.netlib.org/blas/dgemm.f

You will see that:

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

and that:

*  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.

Thus, if A is (2,3), then op(A)=A' is (3,2).

If you look at definition for other arguments, you will see that you must pass M=3, N=3, K=2

Answer 2

In case someone was wondering how the code would look like with different dimensions:

 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);