Fortran gemm function on C-contiguous matrices

Question

I'm trying to do matrix multiplication using the fortran BLAS gemm function, see here.

The signature of this function is, all the parameters' meaning could be found in the above link.

call sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)

My problem is that, I want to use C-contiguous arrays instead of Fortran-contiguous ones, and I've been playing with the above sgemm for quite a while, still badly confused.

Please help me to walk through some concrete examples.

All my input arrays are C-contigous.

a = [[0,1],
     [2,3]]
b = [[0,1,2],
     [3,4,5]]
# pre-alloc memory for c
c = [[0,0,0],
     [0,0,0]]

# compute c = a * b, which should be as follows
# c = [[3,4,5],
#      [9,14,19]]

# since sgemm assumes Fortran-contiguous, so I thought it would be
sgemm('T', 'T', 2, 3, 2, 1.0, a, 2, b, 3, 0, c, 2)
      ~~~~~~~   ~~~~~~~         ~~~   ~~~      ~~~
    trans both   m,n,k          lda   ldb      ldc

# HOWEVER, c is not what I expected, 
c = [[3,9,4],
     [14,5,19]] 

Apparently sgemm stores the elements in Fortran-contiguous order, how to solve this problem? Also I don't quite understand how those m,n,k,lda,ldb are determined if the transa/transb='T' or 'N', hope you could give me a detailed explanation.

NOTE

I'm using this gemm function exported from scipy.linalg.cython_blas, which means, I have no other choice instead of playing this Fortran ordering stuff.

Answer 1

If you want to use row-major matrix instead of Fortran style col-major, you could use the CBLAS API gemm. You can choose matrix storage layout with the first parameter.

https://software.intel.com/en-us/node/520775

---

Or you could still use Fortran API. As changing matrix layout is equivalent to matrix transpose. However you are calculating the transposed C in a wrong way.

Your code calculate C in col-major, but you need a C in row-major. So you need to calculate C^T in col-major by Fortran API, which is equivalent to C in row-major.

It should be

C^T = B^T * A^T

Basically you need to exchange A and B, and corresponding parameters. For more details about those parameters, you could see this answer.

Transpose matrix multiplication in cuBLAS howto