#hw4.txt, February 3, 2014, Frank Wagner
Help:=proc(): print(` MatMultRG(A,B) `): end:


###########Problem 1###########
 
#Small sample of the final part of Garvan's Maple booklet.
#with(linalg):
#A:=matrix(3,4,[1,1,3,-3,5,5,13,-7,3,1,7,-11]);
#A1:=addrow(A,1,2,-5);
#A2:=addrow(A1,1,3,-3);
#A3:=mulrow(A2,3,-1/2);
#A4:=swaprow(A3,2,3);
#A5:=mulrow(A4,3,-1/2);
#A:=matrix(3,3,[177,77,-28,-546,-236,84,-364,-154,51]):
#egenvals(A);
#P:=matrix(3,3,[1,2,-,=3,0,4,-2,13,11]);
#evalm(inverse(P)&*A&*P);
#jordan(A, 'P');
#print(P);



###########Problem 3###########

MatMultRG:=proc(A,B) local n,m1,m,p,A1,B1,M1,i:
n:=nops(A):

if nops(B)<>n then
   RETURN(FAIL):
fi:

if type(log[2](n),integer) then
   MatMultR(A,B):
fi:

m1:=ceil(log[2](n)):
m:=2^(m1):
p:=m-n:

A1:=[seq([op(1..n,A[i]),seq(0,i=1..p)],i=1..n),seq([seq(0,i=1..m)],i=1..p)]:
B1:=[seq([op(1..n,B[i]),seq(0,i=1..p)],i=1..n),seq([seq(0,i=1..m)],i=1..p)]:

M1:=MatMultR(A1,B1):

[seq([op(1..n,M1[i])],i=1..n)]:

end:



###########Problem 4###########

MatSubt:=proc(A,B) local i,j:

if nops(A)<>nops(B) then
 print(`Not the same number of rows`):
 RETURN(FAIL):
fi:

if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all 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 the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A)>0 and nops(A[1])<>nops(B[1]) then
 print(`Not the same number of columns`):
 RETURN(FAIL):
fi:


[seq([seq(A[i][j]-B[i][j],j=1..nops(A[i]))],i=1..nops(A))]:

end:



MatMinus:=proc(A):

[seq([seq(-1*A[i][j],j=1..nops(A[i]))],i=1..nops(A))]:

end:



Strassen:=proc(A,B) local nar,nac,nbr,nbc,n,m1,m,A1,B1,i,A11,A12,A21,A22,B11,B12,B21,B22,
M1,M2,M3,M4,M5,M6,M7,C11,C12,C21,C22,C,A111,A112,A113,A114,A115,B111,B112,B113,B114,B115:
nar:=nops(A):
nac:=nops(A[1]):
nbr:=nops(B):
nbc:=nops(B[1]):
n:=max(nar,nac,nbr,nbc):

if n=1 then
   A[1]*B[1]:
fi:

if nbr<>nac then
   RETURN(FAIL):
fi:

m1:=ceil(log[2](n)):
m:=2^(m1):

A1:=[seq([op(1..nac,A[i]),seq(0,i=1..m-nac)],i=1..nar),seq([seq(0,i=1..m)],i=1..m-nar)]:
B1:=[seq([op(1..nbc,B[i]),seq(0,i=1..m-nbc)],i=1..nbr),seq([seq(0,i=1..m)],i=1..m-nbr)]:

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

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

#Here is where I get stuck; I tried a recursion on this process with a "while" statement, but kept
#crashing Maple on every attempt. So, I just did one process of Strassen (instead of breaking the
#matrix down to numbers) and didn't save any time on the process.

A111:=MatAdd(A11,A22):
B111:=MatAdd(B11,B22):
A112:=MatAdd(A21,A22):
B112:=MatSubt(B12,B22):
B113:=MatSubt(B21,B11):
A113:=MatAdd(A11,A12):
A114:=MatSubt(A21,A11):
B114:=MatAdd(B11,B12):
A115:=MatSubt(A12,A22):
B115:=MatAdd(B21,B22):

M1:=MatMultR(A111,B111):
M2:=MatMultR(A112,B11):
M3:=MatMultR(A11,B112):
M4:=MatMultR(A22,B113):
M5:=MatMultR(A113,B22):
M6:=MatMultR(A114,B114):
M7:=MatMultR(A115,B115):


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


C:=[seq([op(1..m/2,C11[i]),op(1..m/2,C12[i])],i=1..m/2),seq([op(1..m/2,C21[i]),op(1..m/2,C22[i])],i=1..m/2)]:


C:=[seq([op(1..nbc,C[i])],i=1..nar)]:

end: