matrix multiplication as column major

Question

I am trying to multiply as column major and I can't seem to find the right formula! I want to have the matrices as 1D.

Let's say I have these matrices:

A=

1 3

2 4

and B=

5 2 1

6 3 7

The above matrices are assumed that are stored already in column major order.

I am trying:

int main(int argc, const char* argv[]) {
  int rows=2;
  int cols=3;


  int A[rows*rows];
  int B[rows*cols];
  int res[rows*cols];

  A[0]=1;
  A[1]=3;
  A[2]=2;
  A[3]=4;


  B[0]=5;
  B[1]=2;
  B[2]=1;
  B[3]=6;
  B[4]=3;
  B[5]=7;

  /*A[0]=1;
  A[1]=2;
  A[2]=3;
  A[3]=4;


  B[0]=5;
  B[1]=6;
  B[2]=2;
  B[3]=3;
  B[4]=1;
  B[5]=7;
  */


  //multiplication as column major
 for (int i=0;i<rows;i++){
   for (int j=0;j<cols;j++){
     res[i+j*rows]=0;
     for (int k=0;k<rows;k++){
    res[i+j*rows]+=A[i+k*rows]*B[k+j*cols];


   }

   }
}



  for (int i=0;i<rows*cols;i++){
     printf("\n\nB[%d]=%d\t",i,res[i]);

  }

  return 0; 

 }

I am not getting the correct results.

Also,I can't understand (in the case where the matrices are stored in column major already) ,how to index the matrices A and B.

 A[0]=1;

 A[1]=3;

...

or

  A[0]=1;

  A[1]=2;
...

I don't want to transpose the matrices and then use row major.

I want to handle the data as column major.

Because the indices ,if stored as column major,will be different (hence,will matter in order to do the multiplication).

Answer 1

There are two things that lead to your confusion here.

First, the data in your contiguous one-dimensional vector is not in column-major order as you say, but in row-major order, as is the usual layout of two-dimensional contiguous arrays in C. The linear one-dimensional indices of row i and column j in a matrix with M rows and N columns (MxN) are:

A[i*N + j]      // row major
A[i + M*j]      // column major

The "major" refers to the dimension of the outer loop when traversing the array sequentially with two nested loops:

n = 0;
for (i = 0; i < M; i++) {
    for (j = 0; j < N; j++) {
        printf("%8d", A[n++]);
    }
    printf("\n");
}

Second, you use the two dimensions rows and columns which are the dimensions of the resulting matrix, which is confusing, because the number of columns in A is rows.

In fact, there are three different dimensions involved in matrix multiplication when you multiply an MxL matrix A with an LxN matrix B to get an MxN matrix C. In your case, M and L happen to be both 2:

           L (k)   |     N (j)
                   |
                   |   5   2   1
   L (k)           |
                   |   6   3   7
                   |
  -----------------+-------------
                   |
            1   3  |  23  11  22
   M (i)           |
            2   4  |  34  16  30
                   |

The letters in parentheses are the variables the code below uses to iterate over the respective dimension.

Now you can multiply your matrices in row-major format:

    #define M 2
    #define N 3
    #define L 2

    int A[M * L] = {1, 3, 2, 4};
    int B[L * N] = {5, 2, 1, 6, 3, 7};
    int res[M * N];

    int i, j, k;

    for (i = 0; i < M; i++) {
        for (j = 0; j < N; j++) {
            res[j + i * N] = 0;

            for (k = 0; k < L; k++) {
                res[j + i * N] += A[k + i * L] * B[j + k * N];
            }
        }
    }

    for (i = 0; i < M * N; i++) printf("[%d] = %d\n", i, res[i]);

or in column-major format:

    #define M 2
    #define N 3
    #define L 2

    int A[M * L] = {1, 2, 3, 4};
    int B[L * N] = {5, 6, 2, 3, 1, 7};
    int res[M * N];

    int i, j, k;

    for (i = 0; i < M; i++) {
        for (j = 0; j < N; j++) {
            res[j * M + i] = 0;

            for (k = 0; k < L; k++) {
                res[j * M + i] += A[k * M + i] * B[j * L + k];
            }
        }
    }

    for (i = 0; i < M * N; i++) printf("[%d] = %d\n", i, res[i]);

Both input and output are in the respective matrix representation and differ in the two cases, of course.

Answer 2

What do you think about

res[i+j*rows]+=A[i+k*rows]*B[k+j*cols];  

what it will do?

It will access array res, A and B out of bound when i and k becomes 1 and j becomes 2.

res[1+2*2]+=A[1+1*2]*B[1+2*3] = res[5]+=A[4]*B[7];  

This will invoke undefined behavior and you may get either expected or unexpected result.

I think you need this:

res[i*rows+j] += A[i*rows + k] * B[j + k*cols];