Strassen's Algorithm for Matrix multiplication c# implementation

Question

I'm just doing a self-study of Algorithms & Data structures and I'd like to know if anyone has a C# (or C++) implementation of Strassen's Algorithm for Matrix Multiplication?

I'd just like to run it and see what it does and get more of an idea of how it goes to work.

Answer 1

Disclaimer: I haven't tried any of these out, but they appear to be what OP is looking for. These links were just from looking through some Google Code Search results.

I found a C# version. The project doesn't have any frills; it's just the source. However, it appears to be doing the algorithm just from my first cursory scan. In particular, you will want to look at this file.

For C++, I found some code in this google code project. The code is, of course, in English, but the wiki is all in a Cyrillic-written language (Russian?). You will want to look mostly at this file. It appears to have both a serial and and parallel version of Strassen's algorithm.

These projects may not be fully correct, but they are things at which you might want to look more closely.

Answer 2

// Recursive matrix mult by strassen's method.
// 2013-Feb-15 Fri 11:47 by moshahmed/at/gmail.

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

#define M 2
#define N (1<<M)

typedef double datatype;
#define DATATYPE_FORMAT "%4.2g"
typedef datatype mat[N][N]; // mat[2**M,2**M]  for divide and conquer mult.
typedef struct { int ra, rb, ca, cb; } corners; // for tracking rows and columns.
// A[ra..rb][ca..cb] .. the 4 corners of a matrix.

// set A[a] = I
void identity(mat A, corners a){
  int i,j;
  for(i=a.ra;i<a.rb;i++)
    for(j=a.ca;j<a.cb;j++)
      A[i][j] = (datatype) (i==j);
}

// set A[a] = k
void set(mat A, corners a, datatype k){
  int i,j;
  for(i=a.ra;i<a.rb;i++)
    for(j=a.ca;j<a.cb;j++)
      A[i][j] = k;
}

// set A[a] = [random(l..h)].
void randk(mat A, corners a, double l, double h){
  int i,j;
  for(i=a.ra;i<a.rb;i++)
    for(j=a.ca;j<a.cb;j++)
      A[i][j] = (datatype) (l + (h-l) * (rand()/(double)RAND_MAX));
}

// Print A[a]
void print(mat A, corners a, char *name) {
  int i,j;
  printf("%s = {\n",name);
  for(i=a.ra;i<a.rb;i++){
    for(j=a.ca;j<a.cb;j++)
      printf(DATATYPE_FORMAT ", ", A[i][j]);
    printf("\n");
  }
  printf("}\n");
}

// C[c] = A[a] + B[b]
void add(mat A, mat B, mat C, corners a, corners b, corners c) {
  int rd = a.rb - a.ra;
  int cd = a.cb - a.ca;
  int i,j;
  for(i = 0; i<rd; i++ ){
    for(j = 0; j<cd; j++ ){
      C[i+c.ra][j+c.ca] = A[i+a.ra][j+a.ca] + B[i+b.ra][j+b.ca];
    }
  }
}

// C[c] = A[a] - B[b]
void  sub(mat A, mat B, mat C, corners a, corners b, corners c) {
  int rd = a.rb - a.ra;
  int cd = a.cb - a.ca;
  int i,j;
  for(i = 0; i<rd; i++ ){
    for(j = 0; j<cd; j++ ){
      C[i+c.ra][j+c.ca] = A[i+a.ra][j+a.ca] - B[i+b.ra][j+b.ca];
    }
  }
}

// Return 1/4 of the matrix: top/bottom , left/right.
void find_corner(corners a, int i, int j, corners *b) {
  int rm = a.ra + (a.rb - a.ra)/2 ;
  int cm = a.ca + (a.cb - a.ca)/2 ;
  *b = a;
  if (i==0)  b->rb = rm;     // top rows
  else       b->ra = rm;     // bot rows
  if (j==0)  b->cb = cm;     // left cols
  else       b->ca = cm;     // right cols
}

// Multiply: A[a] * B[b] => C[c], recursively.
void mul(mat A, mat B, mat C, corners a, corners b, corners c) {
  corners aii[2][2], bii[2][2], cii[2][2], p;
  mat P[7], S, T;
  int i, j, m, n, k;

  // Check: A[m n] * B[n k] = C[m k]
  m = a.rb - a.ra; assert(m==(c.rb-c.ra));
  n = a.cb - a.ca; assert(n==(b.rb-b.ra));
  k = b.cb - b.ca; assert(k==(c.cb-c.ca));
  assert(m>0);

  if (n==1) {
    C[c.ra][c.ca] += A[a.ra][a.ca] * B[b.ra][b.ca];
    return;
  }

  // Create the 12 smaller matrix indexes:
  //  A00 A01   B00 B01   C00 C01
  //  A10 A11   B10 B11   C10 C11
  for(i=0;i<2;i++) {
  for(j=0;j<2;j++) {
        find_corner(a, i, j, &aii[i][j]);
        find_corner(b, i, j, &bii[i][j]);
        find_corner(c, i, j, &cii[i][j]);
      }
  }

  p.ra = p.ca = 0;
  p.rb = p.cb = m/2;

  #define LEN(A) (sizeof(A)/sizeof(A[0]))
  for(i=0; i < LEN(P); i++) set(P[i], p, 0);

  #define ST0 set(S,p,0); set(T,p,0)

  // (A00 + A11) * (B00+B11) = S * T = P0
  ST0;
  add( A, A, S, aii[0][0], aii[1][1], p);
  add( B, B, T, bii[0][0], bii[1][1], p);
  mul( S, T, P[0], p, p, p);

  // (A10 + A11) * B00 = S * B00 = P1
  ST0;
  add( A, A, S, aii[1][0], aii[1][1], p);
  mul( S, B, P[1], p, bii[0][0], p);

  // A00 * (B01 - B11) = A00 * T = P2
  ST0;
  sub( B, B, T, bii[0][1], bii[1][1], p);
  mul( A, T, P[2], aii[0][0], p, p);

  // A11 * (B10 - B00) = A11 * T = P3
  ST0;
  sub(B, B, T, bii[1][0], bii[0][0], p);
  mul(A, T, P[3], aii[1][1], p, p);

  // (A00 + A01) * B11 = S * B11 = P4
  ST0;
  add(A, A, S, aii[0][0], aii[0][1], p);
  mul(S, B, P[4], p, bii[1][1], p);

  // (A10 - A00) * (B00 + B01) = S * T = P5
  ST0;
  sub(A, A, S, aii[1][0], aii[0][0], p);
  add(B, B, T, bii[0][0], bii[0][1], p);
  mul(S, T, P[5], p, p, p);

  // (A01 - A11) * (B10 + B11) = S * T = P6
  ST0;
  sub(A, A, S, aii[0][1], aii[1][1], p);
  add(B, B, T, bii[1][0], bii[1][1], p);
  mul(S, T, P[6], p, p, p);

  // P0 + P3 - P4 + P6 = S - P4 + P6 = T + P6 = C00
  add(P[0], P[3], S, p, p, p);
  sub(S, P[4], T, p, p, p);
  add(T, P[6], C, p, p, cii[0][0]);

  // P2 + P4 = C01
  add(P[2], P[4], C, p, p, cii[0][1]);

  // P1 + P3 = C10
  add(P[1], P[3], C, p, p, cii[1][0]);

  // P0 + P2 - P1 + P5 = S - P1 + P5 = T + P5 = C11
  add(P[0], P[2], S, p, p, p);
  sub(S, P[1], T, p, p, p);
  add(T, P[5], C, p, p, cii[1][1]);

}
int main() {
  mat A, B, C;
  corners ai = {0,N,0,N};
  corners bi = {0,N,0,N};
  corners ci = {0,N,0,N};
  srand(time(0));
  // identity(A,bi); identity(B,bi);
  // set(A,ai,2); set(B,bi,2);
  randk(A,ai, 0, 2); randk(B,bi, 0, 2);
  print(A, ai, "A"); print(B, bi, "B");
  set(C,ci,0);
  // add(A,B,C, ai, bi, ci);
  mul(A,B,C, ai, bi, ci);
  print(C, ci, "C");
  return 0;
}