load("eigen");

rowdim(A):=length(A)$
coldim(A):=length(A[1])$
size(A):=matrix([length(A),length(A[1])])$

cfloat(x):=float(realpart(x))+float(imagpart(x))*%i$
cbfloat(x):=bfloat(realpart(x))+bfloat(imagpart(x))*%i$

mfloat(x):=block([rx,cx,xx,j1,j2],
 rx:rowdim(x),
 cx:coldim(x),
 xx:copymatrix(x),
for j1 : 1 thru rx do
for j2 : 1 thru cx do 
xx[j1,j2]:cfloat(xx[j1,j2]),
return(xx)
)$


mbfloat(x):=block([rx,cx,xx,j1,j2],
 rx:rowdim(x),
 cx:coldim(x),
 xx:copymatrix(x),
for j1 : 1 thru rx do
for j2 : 1 thru cx do 
xx[j1,j2]:cbfloat(xx[j1,j2]),
return(xx)
)$

gyouop(x,i,j,c):=block([xx],xx:copymatrix(x),
                       xx[i]:x[i]+c*x[j],return(xx))$
retuop(x,i,j,c):=block([xx],xx:transpose
                      (gyouop(transpose(x),i,j,c)),return(xx))$

rowmul(x,i,a):=block([xx],
xx:copymatrix(x),
xx[i]:a*x[i],
return(xx))$

colmul(x,i,a):=transpose(rowmul(transpose(x),i,a))$

tr(x):=block([rx,cx,j1,xx],
 rx:rowdim(x),
 cx:coldim(x),
xx:0,
if rx=cx then (for j1 : 1 thru rx do 
xx:xx+x[j1,j1],
return(xx))
else "not square matrix"
)$


kanyaku(x):=block([cx,rk,xx,j1,j2,j3],
 cx:coldim(x),
 rk:rank(x),
 xx:echelon(x),
if rk<2 then return(xx) else 
for j1 : 2 thru rk do 
block(
for j2 : j1 thru cx do 
if xx[j1,j2]=1 then 
(for j3 : 1 thru j1-1 do 
if xx[j3,j2]#0 then 
 xx[j3]:xx[j3]-xx[j3,j2]*xx[j1],
return(xx))),
return(xx))$

eigenspace(x):=block([eigv,mul,n,r,j1,j2,j3,j4,i1,i2,k0,k1,k2,g,y,z,w,ww,v,vv],
n:length(x),eigv:eigenvalues(x),k0:length(eigv[1]),
for j1 from 1 thru k0 do 
block(kill(g),kill(w),kill(ww),
y:kanyaku(x-eigv[1][j1]*ident(n)),r:rank(y),z:-copymatrix(y),
(for j2 from r+1 thru n do z:submatrix(r+1,z)),
 for j2 from r step -1 thru 1 do
block(for j3 from 1 thru n  
do if y[j2,j3]=1 then (z:submatrix(z,j3),g[j2]:j3,return(z))),
i1:1, i2:1, for j4  from 1 thru n do 
(if j4=g[i1] then (w[j4]:z[i1],i1:i1+1) else (w[j4]:ident(n-r)[i2],i2:i2+1)),
for k1 from 1 thru n do for k2 from 1 thru n-r do ww[k1,k2]:w[k1][k2],
vv:genmatrix(ww,n,n-r,1,1),
print(eigv[1][j1],",",eigv[2][j1],"-ple,",n-r,"-dim,",vv),
if j1=1 then v:vv else v:addcol(v,vv)),return(v))$

Gschmidt(x):=block([y,j],y:gramschmidt(x),
    for j from 1 thru length(y) do
    (y[j]:y[j]/sqrt(conjugate(y[j]).transpose(y[j]))),return(y))$