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.
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.
// 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;
}
