#########################################################################################################################################################
#            																		#
#     The following algorithm computes a canonical set of Jordan chains of a regular matrix polynomial 							#
#			 L(mu) = A_l mu^l + ... + A_1 mu + A_0												#
#     associated with a given eigenvalue lambda.													#
#																			#
#     This is the implementation in Maple of the algorithm presented in the thesis of J. C. Zuniga Anaya (2005)						#
#		"Numerical algorithms for polynomial matrices with applications in control"								#
#########################################################################################################################################################     


#CanonicalSetJC(L,mu,lambda)
#INPUT: 	L- matrix polynomial
#       	mu- variable 
#       	lambda- an eigenvalue of L
#OUTPUT:	A list of elements of the form [k,r,[V_1,...,V_r]] where k represents a partial multiplicity associated with lambda, r is the number of partial multiplicities equal to k, and V_i is a vector of size k*n (n being the size of L) containing a Jordan chain associated with lambda. All Jordan chains returned in the output form a canonical set of Jordan chains at lambda.

CanonicalSetJC:=proc(L,mu,lambda)
	local PM, n , p, succdiff, A, A0, Q, r0, rho0, V, T0, i, j, k, VV, CS, Mul, T, Q1, c, r1, W, rho1, nb, B, JC, M, X, JC1, alpha, t0;
	
	t0:=time();

	PM := [];

	n := LinearAlgebra:-RowDimension(L);
	p := max(map(degree,L,mu));

	succdiff := [L,seq(simplify( 1/k!*map(diff,L,mu$k) ), k=1..p)];

	A0 := map(simplify,convert(evalm(subs(mu=lambda,evalm(L))),Matrix));  
	Q := LinearAlgebra:-NullSpace(A0);                                      
	r0 := n-nops(Q);                 #r0:=rank(A0);	
	rho0 := r0; 
	V := convert([op(Q)],Matrix);    
	
	A := succdiff[2];
	T0 := convert(evalm(subs(mu=lambda,evalm(A))),Matrix);

	i := 2;

	while r0<>n do
		Mul := map(simplify,T0.V);	            
		T := <Mul | A0>; 
		Q1 := LinearAlgebra:-NullSpace(T);  
		c := LinearAlgebra:-ColumnDimension(T);
		r1 := c-nops(Q1);     #r1:=rank(T);
		W := convert([op(Q1)],Matrix);
		rho1 := rho0+r1;   
		nb := r1-r0;
		if nb<>0 then
			JC := [];
			for j from 1 to (i-1)*n-rho0 do
				if LinearAlgebra:-Equal(LinearAlgebra:-SubMatrix(V,1..n,j),Matrix(n,1,shape=zero))=false then 
					JC := [op(JC),map(simplify,convert(evalm(LinearAlgebra:-SubMatrix(V,1..-1,j)),Matrix))];
				end if;
			end do;
			PM := [op(PM),[i-1,nb,JC]];
		end if;
		r0 := r1;
		B := LinearAlgebra:-DiagonalMatrix([V,LinearAlgebra:-IdentityMatrix(n)]);
		V := B.W;
		rho0 := rho1;
		if i<=p then	
			A := succdiff[i+1];		
		else 	 
			A := LinearAlgebra:-ZeroMatrix(n);          
		end if;
		T0 := convert(<convert(evalm(subs(mu=lambda,eval(A))),Matrix)|T0>,Matrix);
		i := i+1;
		Q := Q1;
	end do;

	alpha := [seq([seq([seq(RandomTools:-Generate(integer(range = 1 .. 10)),k=1.. nops(PM[i][3]))],j=1..PM[i][2])],i=1..nops(PM))];
	JC1 := [seq([seq(map(simplify,add(alpha[i][j][k]*LinearAlgebra:-SubMatrix(PM[i][3][k],[1..n],[1]), k=1.. nops(PM[i][3]))),j=1..PM[i][2])],i=1..nops(PM))];
	X := Matrix([seq(seq(JC1[i][j],j=1..PM[i][2]),i=1..nops(PM))]);
	
	while LinearAlgebra:-Rank(X) <> LinearAlgebra:-ColumnDimension(X) do
		alpha := [seq([seq([seq(RandomTools:-Generate(integer(range = 1 .. 10)),k=1.. nops(PM[i][3]))],j=1..PM[i][2])],i=1..nops(PM))];
		JC1 := [seq([seq(add(alpha[i][j][k]*LinearAlgebra:-SubMatrix(PM[i][3][k],[1..n],[1]), k=1.. nops(PM[i][3])),j=1..PM[i][2])],i=1..nops(PM))];
		X := Matrix([seq(seq(JC1[i][j],j=1..PM[i][2]),i=1..nops(PM))]);
	end do;
	
	CS := [seq([PM[i][1],PM[i][2], [seq(map(simplify,add( alpha[i][j][k]*convert(PM[i][3][k],Vector),k=1..nops(PM[i][3]))),j=1..PM[i][2])]], i=1..nops(PM))];
	return (CS, time()-t0);
end proc: