[Solved] How to multiply two matrices having fractions as inputs in c [closed]

by


I wanted to use my own implementation of the Strassen optimization for matrix multiplication. It is highly optimized and hence has almost no pedagogical use but then I remembered that Wikipedia has actual C code in the Strassen-matrix-multiplication entry.

The implementation there is not the best:

  • it has no fallback to the naive algorithm if the size falls below a certain limit
  • only matrix sizes are allowed that are a power of two
  • it accumulates a lot of error if used with finite floating point data types
  • it uses a lot of memory, more than necessary
  • it needs some work to make it parallelizable

But it works, so …

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MM_OK 1
#define MM_ERR 0

void mat_print(double **M, int dim);
double **mat_init(int dim);
void mat_clear(double **M, int dim);
double **mat_rand(int dim, int seed);
int mat_mul_naive(double **A, double **B, double **C, int dim);
int mat_add(double **A, double **B, double **C, int dim);
int mat_sub(double **A, double **B, double **C, int dim);
int mat_mul_strassen(double **A, double **B, double **C, int dim);

// ALL CHECKS OMMITTED!

// Octave/Matlab format, for getting a second opinion
void mat_print(double **M, int dim)
{
  int i, j;
  printf("[");
  for (i = 0; i < dim; i++) {
    for (j = 0; j < dim; j++) {
      printf("%.10g", M[i][j]);
      if (j < dim - 1) {
        putc(',', stdout);
      }
    }
    if (i < dim - 1) {
      putc(';', stdout);
    }
  }
  printf("]\n");
}

double **mat_init(int dim)
{
  int i, j;
  double **M;
  M = malloc(dim * sizeof(double *));
  if (M == NULL) {
    return NULL;
  }
  for (i = 0; i < dim; i++) {
    M[i] = malloc(dim * sizeof(double));
    if (M[i] == NULL) {
      mat_clear(M, i);
      return NULL;
    }
    for (j = 0; j < dim; j++) {
      M[i][j] = 0;
    }
  }
  return M;
}

void mat_clear(double **M, int dim)
{
  int i;
  for (i = 0; i < dim; i++) {
    free(M[i]);
  }
  free(M);
}

double **mat_rand(int dim, int seed)
{
  int i, j;
  double **M;
  M = mat_init(dim);
  if (M == NULL) {
    return NULL;
  }
  srand(seed);
  for (i = 0; i < dim; i++) {
    for (j = 0; j < dim; j++) {
      M[i][j] = (double) rand() / (double) rand();
    }
  }
  return M;
}

int mat_mul_naive(double **A, double **B, double **C, int dim)
{
  int i, j, k;
  if (C == NULL) {
    return MM_ERR;
  }
  for (i = 0; i < dim; i++) {
    for (k = 0; k < dim; k++) {   //<- exchanged for better
      for (j = 0; j < dim; j++) { //<- memory access
        C[i][j] += (A[i][k] * B[k][j]);
      }
    }
  }
  return MM_OK;
}


// helper for mat_mul_strassen, only square matrices of same size
int mat_add(double **A, double **B, double **C, int dim)
{
  int i, j;
  if (C == NULL) {
    fprintf(stderr, "C == NULL in mat_add\n");
    return MM_ERR;
  }
  for (i = 0; i < dim; i++) {
    for (j = 0; j < dim; j++) {
      C[i][j] = A[i][j] + B[i][j];
    }
  }
  return MM_OK;
}

// helper for mat_mul_strassen, only square matrices of same size
int mat_sub(double **A, double **B, double **C, int dim)
{
  int i, j;
  if (C == NULL) {
    fprintf(stderr, "C == NULL in mat_sub\n");
    return MM_ERR;
  }
  for (i = 0; i < dim; i++) {
    for (j = 0; j < dim; j++) {
      C[i][j] = A[i][j] - B[i][j];
    }
  }
  return MM_OK;
}

// This implementation is from Wikipedia
// https://en.wikipedia.org/w/index.php?title=Strassen_algorithm&oldid=498910018#Source_code_of_the_Strassen_algorithm_in_C_language
// And it is a really slow one.
int mat_mul_strassen(double **A, double **B, double **C, int dim)
{
  if (dim == 1) {
    C[0][0] = A[0][0] * B[0][0];
    return MM_OK;
  }
#ifdef MAKE_IT_ACTUALLY_WORK_FAST
  // 256 on my machine, YMMV
  else if (dim <= 256) {
    return mat_mul_naive(A, B, C, dim);
  }
#endif
  else {
    int new_dim = dim / 2;
    double **a11, **a12, **a21, **a22;
    double **b11, **b12, **b21, **b22;
    double **c11, **c12, **c21, **c22;
    double **p1, **p2, **p3, **p4, **p5, **p6, **p7;
    double **aResult;
    double **bResult;
    int i, j;
    a11 = mat_init(new_dim);
    a12 = mat_init(new_dim);
    a21 = mat_init(new_dim);
    a22 = mat_init(new_dim);

    b11 = mat_init(new_dim);
    b12 = mat_init(new_dim);
    b21 = mat_init(new_dim);
    b22 = mat_init(new_dim);

    c11 = mat_init(new_dim);
    c12 = mat_init(new_dim);
    c21 = mat_init(new_dim);
    c22 = mat_init(new_dim);

    p1 = mat_init(new_dim);
    p2 = mat_init(new_dim);
    p3 = mat_init(new_dim);
    p4 = mat_init(new_dim);
    p5 = mat_init(new_dim);
    p6 = mat_init(new_dim);
    p7 = mat_init(new_dim);

    aResult = mat_init(new_dim);
    bResult = mat_init(new_dim);
    // the divide part of divide&conquer
    for (i = 0; i < new_dim; i++) {
      for (j = 0; j < new_dim; j++) {
        a11[i][j] = A[i][j];
        a12[i][j] = A[i][j + new_dim];
        a21[i][j] = A[i + new_dim][j];
        a22[i][j] = A[i + new_dim][j + new_dim];

        b11[i][j] = B[i][j];
        b12[i][j] = B[i][j + new_dim];
        b21[i][j] = B[i + new_dim][j];
        b22[i][j] = B[i + new_dim][j + new_dim];
      }
    }

    // Calculating p1 to p7:

    mat_add(a11, a22, aResult, new_dim);        // a11 + a22
    mat_add(b11, b22, bResult, new_dim);        // b11 + b22
    mat_mul_strassen(aResult, bResult, p1, new_dim);    // p1 = (a11+a22) * (b11+b22)

    mat_add(a21, a22, aResult, new_dim);        // a21 + a22
    mat_mul_strassen(aResult, b11, p2, new_dim);        // p2 = (a21+a22) * (b11)

    mat_sub(b12, b22, bResult, new_dim);        // b12 - b22
    mat_mul_strassen(a11, bResult, p3, new_dim);        // p3 = (a11) * (b12 - b22)

    mat_sub(b21, b11, bResult, new_dim);        // b21 - b11
    mat_mul_strassen(a22, bResult, p4, new_dim);        // p4 = (a22) * (b21 - b11)

    mat_add(a11, a12, aResult, new_dim);        // a11 + a12
    mat_mul_strassen(aResult, b22, p5, new_dim);        // p5 = (a11+a12) * (b22)  

    mat_sub(a21, a11, aResult, new_dim);        // a21 - a11
    mat_add(b11, b12, bResult, new_dim);        // b11 + b12
    mat_mul_strassen(aResult, bResult, p6, new_dim);    // p6 = (a21-a11) * (b11+b12)

    mat_sub(a12, a22, aResult, new_dim);        // a12 - a22
    mat_add(b21, b22, bResult, new_dim);        // b21 + b22
    mat_mul_strassen(aResult, bResult, p7, new_dim);    // p7 = (a12-a22) * (b21+b22)

    // calculating c21, c21, c11 e c22:

    mat_add(p3, p5, c12, new_dim);      // c12 = p3 + p5
    mat_add(p2, p4, c21, new_dim);      // c21 = p2 + p4

    mat_add(p1, p4, aResult, new_dim);  // p1 + p4
    mat_add(aResult, p7, bResult, new_dim);     // p1 + p4 + p7
    mat_sub(bResult, p5, c11, new_dim); // c11 = p1 + p4 - p5 + p7

    mat_add(p1, p3, aResult, new_dim);  // p1 + p3
    mat_add(aResult, p6, bResult, new_dim);     // p1 + p3 + p6
    mat_sub(bResult, p2, c22, new_dim); // c22 = p1 + p3 - p2 + p6

    // Grouping the results obtained in a single matrix:
    for (i = 0; i < new_dim; i++) {
      for (j = 0; j < new_dim; j++) {
        C[i][j] = c11[i][j];
        C[i][j + new_dim] = c12[i][j];
        C[i + new_dim][j] = c21[i][j];
        C[i + new_dim][j + new_dim] = c22[i][j];
      }
    }
    mat_clear(a11, new_dim);
    mat_clear(a12, new_dim);
    mat_clear(a21, new_dim);
    mat_clear(a22, new_dim);

    mat_clear(b11, new_dim);
    mat_clear(b12, new_dim);
    mat_clear(b21, new_dim);
    mat_clear(b22, new_dim);

    mat_clear(c11, new_dim);
    mat_clear(c12, new_dim);
    mat_clear(c21, new_dim);
    mat_clear(c22, new_dim);

    mat_clear(p1, new_dim);
    mat_clear(p2, new_dim);
    mat_clear(p3, new_dim);
    mat_clear(p4, new_dim);
    mat_clear(p5, new_dim);
    mat_clear(p6, new_dim);
    mat_clear(p7, new_dim);
    mat_clear(aResult, new_dim);
    mat_clear(bResult, new_dim);

  }
  return MM_OK;
}

static int ispow2(unsigned int x)
{
  while (((x & 1) == 0) && x > 1) {
    x >>= 1;
  }
  if (x == 1) {
    return 1;
  }
  return 0;
}

#include <time.h>
int main(int argc, char **argv)
{
  //int res;
  int dim1;
  double **A, **B, **C, **D;
  clock_t start, stop;

  // only square matrices of same size
  if (argc != 2) {
    fprintf(stderr, "Usage: %s dimension\n", argv[0]);
    exit(EXIT_FAILURE);
  }

  dim1 = atoi(argv[1]);
  if (!ispow2(dim1)) {
    fprintf(stderr,
            "Dimension of matrix must be a power of two and %d is not\n", dim1);
    exit(EXIT_FAILURE);
  }

  A = mat_rand(dim1, 123);
  B = mat_rand(dim1, 456);

  //mat_print(A, dim1);
  //mat_print(B, dim1);

  C = mat_init(dim1);
  start = clock();
  mat_mul_naive(A, B, C, dim1);
  stop = clock();
  printf("Naive    %.10f seconds\n", (double) (stop - start) / CLOCKS_PER_SEC);
  //mat_print(C, dim1);

  D = mat_init(dim1);
  start = clock();
  mat_mul_strassen(A, B, D, dim1);
  stop = clock();
  printf("Strassen %.10f seconds\n", (double) (stop - start) / CLOCKS_PER_SEC);
  //mat_print(D, dim1);

  mat_sub(C, D, C, dim1);
  // results will differ more and more the larger the matrix is
  //mat_print(C, dim1);

  mat_clear(A, dim1);
  mat_clear(B, dim1);
  mat_clear(C, dim1);
  mat_clear(D, dim1);
  exit(EXIT_SUCCESS);
}

solved How to multiply two matrices having fractions as inputs in c [closed]