function test_strass()
clear all
close all
cases = [100, 128, 500, 512, 800, 1024, 1453];
count = zeros(7,1)
for ii = 1:7
    n = cases(ii)
    cost1(ii) = (2*n-1)*n^2
    cost2(ii) = 2*n^2.8;
%   n = 2^ii;
%   n = 16;
    A = zeros(n,n); B = A;
    A = magic(n);
    for i = 1:n
        for j = 1:n
            B(i,j) = cos((i-1)*(j-1)*pi/(n^2));
        end
    end
    countin = 0;
    C1 = A*B;
    [C2,count(ii,1)] = strass(A,B,100,countin);
    count(ii)
    ndiff = norm(C2-C1)/norm(C1)
end
xx = (1:ii)';
whos
semilogy(xx,count,'r-',xx,cost1',xx,cost2','g-')
figure
plot(xx, ndiff)

function [C,countout] = strass(A,B,nmin,countin)
% A "divide-and-conquer" recursive algorithm for matrix multiplication
% as n grows, scales like n^2.807 compared to n^3 for direct multiplication
[n,k] = size(A);
if n <= nmin
    C = A*B;
    countout = countin + (2*n-1)*n^2;
else
    m = ceil(n/2);
    if 2*m ~= n
        A = [A, zeros(n,   1)];
        A = [A; zeros(1, n+1)];
        B = [B, zeros(n,   1)];
        B = [B; zeros(1, n+1)];
    end
    u = 1:m; v = m+1:2*m;
    [P1,count1] = strass(A(u,u)+A(v,v), B(u,u)+B(v,v), nmin, countin);
    [P2,count2] = strass(A(v,u)+A(v,v), B(u,u)       , nmin, count1);
    [P3,count3] = strass(A(u,u)       , B(u,v)-B(v,v), nmin, count2);
    [P4,count4] = strass(A(v,v)       , B(v,u)-B(u,u), nmin, count3);
    [P5,count5] = strass(A(u,u)+A(u,v),        B(v,v), nmin, count4);
    [P6,count6] = strass(A(v,u)-A(u,u), B(u,u)+B(u,v), nmin, count5);
    [P7,count7] = strass(A(u,v)-A(v,v), B(v,u)+B(v,v), nmin, count6);
    C0(u,u) = P1 + P4 - P5 + P7;
    C0(u,v) = P3 + P5;
    C0(v,u) = P2 + P4;
    C0(v,v) = P1 + P3 - P2 + P6;
    countout = count7 + 18*m^2;
    C = C0(1:n,1:n);
end
