To execute an example as given below , one has to define a proper Maple environment:restart;with(LinearAlgebra): with(MatrixPolynomialAlgebra):Next, one has to activate all Maple procedures, defined in the following section.

LCMDenomMatrixPolynom:=proc(W,x) local k,m,nc,D0,D1,cofmax,mat; mat := convert(W(x),Matrix): nc := ColumnDimension(mat): D0:=[]:for k from 1 to nc do for m from 1 to nc do D0:=[op(D0),denom(mat[k,m])]; end do end do; D1:=lcm(op(D0)): cofmax:=coeff(D1,x,degree(D1,x)): return(expand(D1/cofmax));end;GetPolesandZeros:=proc(pf,zf) local poles, zeros, mu; poles:=[solve(pf(x),x)]: print('poles'=poles); zeros:=[solve(zf(x),x)]: print('zeros'=zeros); mu:=convert([op({op(convert(poles,list)), op(convert(zeros,list))})],Vector[row]): print(`Set of different poles and zeros`=mu); return([convert(poles,Vector[row]), convert(zeros,Vector[row]),mu]); end proc;GetMultiplicity:=proc(p::Vector,mu::Vector) local dimp,dimgv,mP,k,m; dimp:=Dimension(p): dimgv:=Dimension(mu): mP:=Vector[row](dimgv,0): for k from 1 to dimgv do for m from 1 to dimp do if evalb(p[m]=mu[k]) then mP[k]:=mP[k]+1; end if: end do:end do: return(mP);end;GetAllMOrderings:=proc(mA::Vector,mZ::Vector,mu::Vector) local h,kordering,orderings, newperm; newperm:=true:h:=0:kordering:=0:while (kordering=0) do h:=h+1:orderings:=GetMOrderingsH(mA,mZ,mu,h,newperm): newperm:=false: kordering:=orderings[1]: end do: return(h,orderings);end;GetMOrderingsH:=proc(muA,muZ,mu,h,newperm) local nmu,nperm,kk,k,perm,orderingP,orderingZ; global allperm; nmu:= Dimension(mu): if newperm=true then allperm:= combinat[permute](nmu): end if: nperm:= combinat[numbperm](nmu): kk:=0: for k from 1 to nperm do perm:=allperm[k]: if (TestOrderingMAZ(muA,muZ,perm,h)) then kk:=kk+1: orderingP[kk]:=GetOrderedVector(muA,mu,convert(perm,list)): orderingZ[kk]:=GetOrderedVector(muZ,mu,convert(perm,list)): end if: end do: if kk>0 then return(kk,orderingP,orderingZ) else return(kk); end if: end proc;TestOrderingMAZ:= proc(ma,mz,perms,h) local s, bol, t; s:=Dimension(mz): bol:=true: t:=0: while (t<s) and bol do t:=t+1: bol:=not((add(mz[perms[m1]],m1=1..t)- add(ma[perms[m2]],m2=1..(t-1)))>(h)): end do;return(bol);end;GetOrderedVector:=proc(ma,mu,ordering) local no,n1,tk,nc, k, m, mup, mpp, orderedV; mup:=mu[ordering]: mpp:=ma[ordering]: no:=Dimension(mu): n1:=add(ma[k],k=1..no): nc:=0: orderedV:=Vector[row](n1,0): for k from 1 to no do tk:=mpp[k]: for m from 1 to tk do nc:=nc+1: orderedV[nc]:=mup[k]: end do: end do: return(orderedV);end proc;Tcomplete:=proc(oZ::Vector,oA::Vector) local k,nc,S; nc:=Dimension(oA): if not (Dimension(oZ)=nc) then error "Input vectors should have equal length"; end if: S:=Scol(oZ,oA,0): for k from 2 to nc do S:=<S|Scol(oZ,oA,k-1)>:end do: return(Transpose(S));end;Scol:=proc(oZ,oA,j) local pol,k,nc,vj,x; nc:=Dimension(oA);vj:=Vector(nc):pol:=1: if (j>0) then for k from 1 to j do pol:=pol*(x-oA[k]) end do: end if: if ((j+1)<nc) then for k from (j+1) to (nc-1) do pol:=pol*(x-oZ[k+1]): end do end if: for k from 1 to nc do vj[k]:=coeff(pol,x,k-1);end do: return(vj);end;MakeFactorization:=proc(A,B,C,x) local WW, nv, m, Im, k, R, oa; m:=RowDimension(C);nv:=Dimension(A):Im:=IdentityMatrix(m): WW:=Vector(nv): for k from 1 to nv[1] do R:=MakeRmatrix(B,C,k): WW[k]:=unapply(map(factor,Im+ScalarMultiply(R,1/(x-A[k,k]))),x): end do: return(WW);end;MakeRmatrix:=proc(B,C,k) local nc, mat, i, j; nc:=RowDimension(C); mat:=Matrix(nc,nc,0); for i from 1 to nc do for j from 1 to nc do mat[i,j]:=C[i,k]*B[k,j]; end do end do: return(mat);end;Factors2Transfer := proc(AllFactors,x) local Wtest,DimS,k,n,Wdum,ResultW; DimS:=ColumnDimension(AllFactors[1](x)): n:=Dimension(AllFactors):Wtest:=IdentityMatrix(DimS): for k from 1 to n do Wtest:=map(simplify,Wtest.AllFactors[k](x)): end do: Wdum:=convert(map(factor,map(simplify,evalm(Wtest))),Matrix): ResultW:=unapply(Wdum,x): return(ResultW);end proc;#End of procedures

Define rational matrix function W:W:=x-><<1,0>|<1/((x-1)*(x-3)^4*(x-4)*(x-5)^3*(x-6)), (x-2)^2*(x-6)/((x-1)*(x-3)*(x-5))>>: 'W(lambda)'=W(lambda);q:=unapply(LCMDenomMatrixPolynom(W,x),x):ppoles:=q:pzeros:=unapply(simplify(W(x)[2,2]*q(x)),x):'p(lambda)'=sort(collect(ppoles(lambda),lambda),lambda);'(p^(x))(lambda)'=sort(collect(pzeros(lambda),lambda),lambda);res1:=GetPolesandZeros(ppoles,pzeros):poles:=res1[1]:zeros:=res1[2]:mu:=res1[3]:npoles:=Dimension(poles);nzeros:=Dimension(zeros); nmu:=Dimension(mu);muA:=GetMultiplicity(poles,mu);muZ:=GetMultiplicity(zeros,mu);Get orderings:AllMOrderings:=GetAllMOrderings(muA,muZ,mu):AllMOrderings;orderedA:=AllMOrderings[3][6]:`alpha`=orderedA;orderedZ:=AllMOrderings[4][6]:`zeta`=orderedZ;Construct companion based realizationr:=unapply(simplify(W(lambda)[1,2]*ppoles(lambda)),lambda): Amin:=Transpose(CompanionMatrix(ppoles(lambda),lambda)); Bmin:=Matrix(npoles,2):Bmin[npoles,2]:=1: Cmin:=Matrix(2,npoles,0): for k from 1 to npoles do Cmin[1,k]:=coeff(r(lambda),lambda,k-1): Cmin[2,k]:=coeff(pzeros(lambda)-ppoles(lambda),lambda,k-1): end do: Dmin:=IdentityMatrix(2,2): Amincross:=Transpose(CompanionMatrix(pzeros(lambda),lambda));Triangularization:Tr:=Tcomplete(orderedZ,orderedA):`T`=Tr;Trinv:=MatrixInverse(Tr): Atr:=Tr.Amin.Trinv:`A`=Atr; Btr:=Tr.Bmin:Ctr:=Cmin.Trinv: Atrcross:=Tr.(Amincross).Trinv:'A^(x)'=Atrcross;Elementary factors:Allfactors:=map(simplify,MakeFactorization(Atr,Btr,Ctr,lambda)):afactors:=Vector[row](npoles,0): for k from 1 to npoles do afactors[k]:=Allfactors[k](lambda): end do: print(`Elementary factors`=afactors);Wtest:=Factors2Transfer(Allfactors,lambda): print(`Product elementary factors`=Wtest(lambda));

MakeCompleteFactors:=proc(Amin,Bmin,Cmin,orderP,orderZ,x) local Tr,Trinv,Atr,Btr,Ctr,Allfactors,afactors,k,np; np:=Dimension(orderP): Tr:=Tcomplete(orderZ,orderP): Trinv:=MatrixInverse(Tr): Atr:=Tr.Amin.Trinv: Btr:=Tr.Bmin:Ctr:=Cmin.Trinv: Allfactors:=map(simplify,MakeFactorization(Atr,Btr,Ctr,x)): afactors:=Vector[row](np,0): for k from 1 to np do afactors[k]:=Allfactors[k](x): end do: print(`Elementary factors`=afactors);end proc:Example how to use MakeCompleteFactorsorderP:=AllMOrderings[3][5]:orderZ:=AllMOrderings[4][5]: print(`ordering poles`=orderP);print(`ordering zeros`=orderZ); MakeCompleteFactors(Amin,Bmin,Cmin,orderP,orderZ,lambda);