I have this code to multiply two matrices stored in 1D arrays.

#include <stdlib.h>
#include <assert.h>
#include <stdint.h>
#include <stdio.h>

void matrix_multiplier(int n, int m, int k,
            double* C, const double* A, const double* B) {
  
    for (int i = 0; i < n; i  ) {
        for (int j = 0; j < k; j  ) {
            int sum = 0;
            for (int p = 0; p < m; p  )
                sum = sum   A[p   i * k ] * B[j   p * n];
            C[j   i * n] = sum;
        }
    }
}



int main(){

  int n = 3;
  int m = 3;
  int k = 2;
    
  double A[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};

  double B[9] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};

  double *C = calloc(n*k, sizeof(double));


  matrix_multiplier(n, m, k, C, A, B);


  printf("\n");
  printf("input matrices\n");

   for (int j = 0; j < n;   j) {
        for (int i = 0; i < m;   i) {
             printf("%f ", A[j * m  i]);
        }
        printf("\n");
    }


  printf("\n");


  for (int j = 0; j < n;   j) {
        for (int i = 0; i < m;   i) {
             printf("%f ", B[j * m  i]);
        }
        printf("\n");
    }


   printf("output matrix\n");

   for (int j = 0; j < n;   j) {
        for (int i = 0; i < m;   i) {
             printf("%f ", C[j * m  i]);
        }
        printf("\n");
    }


     free(C);

    return 0;
}

Where A has size n x m, B has size m x k, C has size n x k.

However I'm not sure sum = sum A[p i * k ] * B[j p * n]; line is entirely correct. I tried following some pseudocode and it led me only to this. Furthermore, I'm not quite sure if this code is storing the data in row-major or column-major order, cause i do actually need it to be in row-major.

CodePudding user response：

Using integer for accumulation of doubles is a bad idea: Replace:

int sum = 0;

with

double sum = 0;

Another issue is indexing the C array. It should be:

C[j   i * k] = ...

For every increment of i by 1, the j is incremented k times.