#C5.txt:
#Feb. 6, 2014: Raleigh's power method etc.
Help:=proc(): print(` RPow(A,k,x0), Ons1(z) `): end:



#####stuff from Monday's "class"
#C4.txt: Feb. 3 2014. Distance-learning class (due to University closing due to snow storm)
Help4:=proc(): print(` MatAdd(A,B),MatMult(A,B), PadZeros(A), ,MatMultR(A,B) `) : 
print(`Tra(A), RandMat(m,n,R)`):
end:


#Tra(A): the transpose of the matrix A
Tra:=proc(A) local i,j:
[seq([seq(A[i][j],i=1..nops(A))],j=1..nops(A[1]))]:
end:


#RandMat(m,n,R): a random m by n matrix with entries between -R and R
#given as a list of lists
RandMat:=proc(m,n,R) local ra,i,j:
ra:=rand(-R..R):
[seq([seq(ra(),j=1..n)],i=1..m)]:
end:


#MatAdd(A,B): adds two matrices of the same size given as lists of lists
MatAdd:=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:


#MatMult(A,B): multiplies two matrices A and B 
MatMult:=proc(A,B) local i,j,k:


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[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:


[seq([seq( add(A[i][k]*B[k][j],k=1..nops(A[i])), 
j=1..nops(B[1]))],i=1..nops(A))]:

end:

#PadZeros(A): given  a square matrix A pads it with zero to make its
#dimension a powr of 2
PadZeros:=proc(A) local i,n,n1:

if nops(A)<>nops(A[1]) then
 print(`not a square matrix`):
 RETURN(FAIL):
fi:

n:=nops(A):

if type(log[2](n),integer) then
A:
else
n1:=2^ceil(log[2](n)):
[seq([op(A[i]),0$(n1-n)],i=1..nops(A)),[0$n1]$(n1-n)]:
fi:

end:


#MatMultR(A,B): multiplies two square matrices A and B whose dimension
#is a power of 2
MatMultR:=proc(A,B) local n,i,A11,A12,A21,A22,B11,B12,B21,B22, 
M1,M2,M3,M4,M5,M6,M7,M8,C11,C12,C21,C22:
n:=nops(A):

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

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

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

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:=MatMultR(A11,B11):
M2:=MatMultR(A12,B21):
M3:=MatMultR(A11,B12):
M4:=MatMultR(A12,B22):


M5:=MatMultR(A21,B11):
M6:=MatMultR(A22,B21):
M7:=MatMultR(A21,B12):
M8:=MatMultR(A22,B22):
 
C11:=MatAdd(M1,M2):
C12:=MatAdd(M3,M4):
C21:=MatAdd(M5,M6):
C22:=MatAdd(M7,M8):

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

end:


#####end stuff from Monday's "class"


#RPow(A,k,x0): Raleigh quotient Power method for approximating
#the dominant eigenvalue of a matrix A (given as a list of lists)
RPow:=proc(A,k,x0) local x1,n,i,j,ma:

n:=nops(A):



 x1:=x0:

for i from 1 to k do
x1:=MatMult(A,x1):
ma:=max( seq( abs(x1[j][1]), j=1..n)):
x1:=[ seq([x1[j][1]/ma], j=1..n)]:
od:

MatMult(Tra(x1), MatMult(A,x1))[1][1]/MatMult(Tra(x1),x1)[1][1]:


end:



#Ons1(): the 2 by 2 transfer matrix for the 1 by infinity
 #Ising model columns are [1,-1], rows are [1,-1]
Ons1:=proc(z) 
[[z,1/z],[1/z,z]]:
end: