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<N-1 do 

   cc:=counter;
   counter:=counter+1; 
   vects:=<<ibb>>;

 if cc>0 then 
		WZ:=Matrix(cc,N,(i,j)-> if i=j then 1 else 0 fi);
		WA:=<<WZ,vects>>;  
	else 
		WA:=Transpose(vects); 
  fi;

	WA:=map(evalc,WA);

            NS:=NullSpace(WA);
	ns:=nops(NS);


	try

            	WW||counter:=map(evalc,<<Transpose(WA) | seq(NS[a],a=1..ns)>>);
		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(<<z,Matrix(ib[1])>>));
	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";