with(LinearAlgebra): UT:=proc(PXX,QXX) global W;local WWCI,PX,QX,N,WA,WZ,cc,counter,ne,nf,ib,i,j,ibb,vects,NS,ns,z,er; PX:=Matrix(PXX); QX:=Matrix(QXX); N:=RowDimension(PXX): counter:=0: cc:=0: P||0:=PX: Q||0:=QX: ne:=map(evalc,Eigenvectors(PX,output=list)); nf:=map(evalc,Eigenvectors(QX,output=list)); ib:={}: for i to nops(ne) do if nops(ib)>0 then break fi; for j to nops(nf) do ib:=IntersectionBasis([ ne[i,3],nf[j,3] ]); if nops(ib)>0 then break fi; od;od; if counter=0 then ibb:=ib[1] fi; while counter>; if cc>0 then WZ:=Matrix(cc,N,(i,j)-> if i=j then 1 else 0 fi); WA:=<>; else WA:=Transpose(vects); fi; WA:=map(evalc,WA); NS:=NullSpace(WA); ns:=nops(NS); try WW||counter:=map(evalc,<>); WWCI:=(WW||counter)^(-1); catch: return [rank(WW||counter)]; end try; P||counter:=map(evalc,WWCI.P||cc.WW||counter); Q||counter:=map(evalc,WWCI.Q||cc.WW||counter); PP||counter:=map(evalc,SubMatrix(P||counter,counter+1..N,counter+1..N)); QQ||counter:=map(evalc,SubMatrix(Q||counter,counter+1..N,counter+1..N)); ne:=map(evalc,Eigenvectors(PP||counter,output=list)); nf:=map(evalc,Eigenvectors(QQ||counter,output=list)); ib:={}: for i to nops(ne) do if nops(ib)>0 then break fi; for j to nops(nf) do ib:=IntersectionBasis([ ne[i,3],nf[j,3] ]); if nops(ib)>0 then break fi; od;od; z:=ZeroMatrix(counter,1); try ibb:=map(evalc,Transpose(<>)); catch: return [0]; end try; od; W:=`.`(seq(WW||k,k=1..N-1)); [map(evalf,map(evalc,map(simplify,W))), map(evalf,map(evalc,map(simplify,W^(-1).PX.W))), map(evalf,map(evalc,map(simplify,W^(-1).QX.W)))] end:"UT(R,B) will simultaneously upper-triangularize two Matrices R and B"; UTC:=proc(A,B) local T; T:=UT(Matrix(A),Matrix(B)); nops(T) end:"UTC will test for simultaneous upper-triangularizability";