# hw4.txt, 2014-02-09, Jeff Ames

with(LinearAlgebra):

MatMultRG := proc(A, B)
    local i, m, n, M:
    m := nops(A):
    n := nops(B[1]):
    M := MatMult(PadZeros(A), PadZeros(B)):
    for i from 1 to m do
        M[i] := M[i][..n]:
    od:
    M[..m]:
end:

# Test MatMultRG versus Maple's built-in matrix multiplication n times
# on matrices of size sz x sz, with values in [-range, range].
TestMatMultRG := proc(sz, range, n)
    local i, A, B, ours, maples, result:
    result := true:
    for i from 1 to n do
        A := RandMat(sz, sz, range):
        B := RandMat(sz, sz, range):
        ours := Matrix(MatMultRG(A, B)):
        maples := Matrix(A) . Matrix(B):
        if not Equal(ours, maples) then
            result := false:
        fi:
    od:
    result:
end:

MatSubt := proc(A, B)
    MatAdd(A, MinusMat(B)):
end:

MinusMat := proc(A)
    local n:
    n := nops(A):
    [seq([op(1..n, -A[i])], i = 1..n)]:
end:

Strassen := proc(A, B)
    local i, j, k, n,
          A11, A12, A21, A22,
          B11, B12, B21, B22,
          M1, M2, M3, M4, M5, M6, M7,
          C11, C12, C21, C22:

    if nops({seq(nops(A[i]), i=1..nops(A))}) <> 1 then
        print(`Rows do not all have the same lengths, hence`, A, `is not a matrix `):
        return FAIL:
    fi:

    if nops({seq(nops(B[i]), i=1..nops(B))}) <> 1 then
        print(`Rows do not all have the same lengths, hence`, B, `is not a matrix `):
        return FAIL:
    fi:

    if nops(A[1]) <> nops(B) then
        print(`Number of columns of`, A, `is not the same as the number of rows of`, B):
        print(`Hence can't multiply them.`):
        return FAIL:
    fi:

    n := nops(A):

    if not type(log[2](n), integer) then
        return FAIL:
    fi:

    if n = 1 then
        # Base case: 1x1 matrices, so just use normal multiplication
        [[A[1][1] * B[1][1]]]:
    else
        # Recursive case: break matrices into 4 submatrices and
        # multiply using Strassen's algorithm

        A11 := [seq([op(1..n/2, A[i])], i = 1..n/2)]:
        A12 := [seq([op(n/2+1..n, A[i])], i = 1..n/2)]:
        A21 := [seq([op(1..n/2, A[i])], i = n/2+1..n)]:
        A22 := [seq([op(n/2+1..n, A[i])], i = n/2+1..n)]:

        B11 := [seq([op(1..n/2, B[i])], i = 1..n/2)]:
        B12 := [seq([op(n/2+1..n, B[i])], i = 1..n/2)]:
        B21 := [seq([op(1..n/2, B[i])], i = n/2+1..n)]:
        B22 := [seq([op(n/2+1..n, B[i])], i = n/2+1..n)]:

        M1 := Strassen(MatAdd(A11, A22), MatAdd(B11, B22)):
        M2 := Strassen(MatAdd(A21, A22), B11):
        M3 := Strassen(A11, MatSubt(B12, B22)):
        M4 := Strassen(A22, MatSubt(B21, B11)):
        M5 := Strassen(MatAdd(A11, A12), B22):
        M6 := Strassen(MatSubt(A21, A11), MatAdd(B11, B12)):
        M7 := Strassen(MatSubt(A12, A22), MatAdd(B21, B22)):

        C11 := MatAdd(MatSubt(MatAdd(M1, M4), M5), M7):
        C12 := MatAdd(M3, M5):
        C21 := MatAdd(M2, M4):
        C22 := MatAdd(MatAdd(MatSubt(M1, M2), M3), M6):

        [seq([op(C11[i]), op(C12[i])], i = 1..nops(C11)),
         seq([op(C21[i]), op(C22[i])], i = 1..nops(C21))]:
    fi:
end:

# Test Strassen versus Maple's built-in matrix multiplication n times
# on matrices of size sz x sz, with values in [-range, range].
TestStrassen := proc(sz, range, n)
    local i, A, B, ours, maples, result:

    result := true:

    A := [[1]]:
    B := [[1]]:
    if Strassen(A, B) <> [[1]] then result := false fi:

    A := [[1,0],[0,1]]:
    B := [[1,0],[0,1]]:
    if Strassen(A, B) <> [[1,0],[0,1]] then result := false fi:

    for i from 1 to n do
        A := RandMat(sz, sz, range):
        B := RandMat(sz, sz, range):
        ours := Matrix(Strassen(A, B)):
        maples := Matrix(A) . Matrix(B):
        if not Equal(ours, maples) then
            result := false:
        fi:
    od:
    result:
end:

RandMat := proc(n, R)
    local ra,i,j:
    ra := rand(0..R):
    [seq([seq(ra(), j = 1..n)], i = 1..n)]:
end:

# Test Strassen versus MatMultR n times on matrices of size sz x sz,
# with values in [0, range]. Return the total time for Strassen and
# MatMultR as a 2-element array.
BenchmarkStrassen := proc(sz, range, n)
    local i, stime, mtime, A, B:
    stime := 0:
    mtime := 0:
    A := RandMat(sz, range):
    B := RandMat(sz, range):
    for i from 1 to n do
        stime := stime + time(Strassen(A, B)):
        mtime := mtime + time(MatMultR(A, B)):
    od:
    [stime, mtime]:
end:

# Running "for i from 1 to 12 do BenchmarkStrassen(2^i, 50, 1) od;" gives these results:
#
# 2x2 => [0., 0.]
# 4x4 => [0.001, 0.001]
# 8x8 => [0.012, 0.006]
# 16x16 => [0.101, 0.052]
# 32x32 => [0.636, 0.504]
# 64x64 => [4.530, 3.819]
# 128x128 => [32.014, 30.845]
# 256x256 => [230.258, 256.891]
#
# So Strassen is slower for all of these, but it is possibly faster as
# n increases beyond some threshold.