##----------------------------------------------------------------------
##    Title: 	PencilReduction
##    Created:	March 2010
##    Authors: 	Luu Ba Thang <lbathang@sophia.inria.fr>
##    Advisor:  Laurent Buse <Laurent.Buse@inria.fr>  
## ----------------------------------------------------------------------
##  Description: Some functions dealing with linearization of a polynomial matrix, 
##               extracting the regular part of a non square pencil of matrices and computing eigenvalues.

##-------------------------------------------------------------------------
##   Version last modified 04/2010
##-------------------------------------------------------------------------




__ := NULL: 
addtohelp := proc() 
  global __; 
  __ := __, args: 
end:


with(LinearAlgebra):
with(MatrixPolynomialAlgebra): 
with(PolynomialTools):
with(Groebner):

######## Extracting a regular part of a pencil matrix A-tB ##########

#----------------------------------------------------------------------
addtohelp(FirstCompanionMatrix):
addtohelp(SecondCompanionMatrix):
#----------------------------------------------------------------------
LinearizationMatrix := proc (M, var) 
   local r, s, A, B, U, V, i; 
    s := Degree(M, var); 
    r := RowDimension(M); 
    A := Transpose(Coeff(M, var, 0)); 
    i := 1; 
    while i < s do 
     B := Transpose(Coeff(M, var, i)); 
     A := `<|>`(A, B); i := i+1 
    end do; 
  RETURN(A);
end proc:

## HELP FirstCompanionMatrix; SecondCompanionMatrix
## CALLING SEQUENCE:
##     FirstCompanionMatrix(M,var); SecondCompanionMatrix(M,var)
## DESCRIPTION:
## - Given a univariate polynomial matrix M.
## - var - : variable
## RETURN: A, B: Companion matrices of M. 
## EXAMPLES: M := Matrix(2, 2, {(1, 1) = t^2, (1, 2) = t+1, (2, 1) = 1, (2, 2) = t})
## FirstCompanionMatrix(M, var)
##     Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0,
##                   (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1,
##                   (3, 1) = 0, (3, 2) = 1, (3, 3) = 0, (3, 4) = 0,
##                   (4, 1) = 1, (4, 2) = 0, (4, 3) = 1, (4, 4) = 1})
##  SecondCompanionMatrix(M, var)
##      Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0,
##                    (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, 
##                    (3, 1) = 0, (3, 2) = 0, (3, 3) = -1, (3, 4) = 0, 
##                    (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0})
FirstCompanionMatrix := proc (M, var) 
   local m, d1, d, U, V, W, A; 
       m := RowDimension(M); 
       d := Degree(M, var); 
      d1 := (d-1)*m; 
       U := Matrix(d1, m); 
       V := IdentityMatrix(d1); 
       W := `<|>`(U, V); 
       A := `<,>`(W, LinearizationMatrix(M, var)); 
RETURN(A);
end proc:

SecondCompanionMatrix := proc (M, var) 
    local m, n, d, d1, U, V, W, T, B1, B2, B; 
      m := RowDimension(M); 
      n := ColumnDimension(M); 
      d := Degree(M, var); 
     d1 := (d-1)*m; 
      U := Matrix(n, d1); 
      V := IdentityMatrix(d1); 
      W := Matrix(d1, m); 
      T := -Transpose(Coeff(M, var, d));
     B1 := `<|>`(V, W);
     B2 := `<|>`(U, T); 
      B := `<,>`(B1, B2); 
    RETURN(B); 
end proc:

ReductionStep := proc (A, B) 
   local m, n, p, r, s, B1, B2, B3, B4, B5, A1, A2, A3, A4, A5, C, C1, C2, P, L, U, P1, L1, U1, E1, E2, R, R1, E, F, M, N; 
      m := RowDimension(B); 
      n := ColumnDimension(B); 
      r := Rank(B); 
      M := A; N := B; 
    while r < n do 
       B1 := Transpose(N); 
       P, L, U, R := LUDecomposition(B1, method = 'RREF'); 
       C := Multiply(Multiply(P, L), U); 
       C1 := Transpose(C); 
       C2 := MatrixInverse(C1,method = LU); 
       B2 := Multiply(N, C2); 
       A1 := Multiply(M, C2); 
       A2 := DeleteColumn(A1, 1 .. r); 
        s := Rank(A2); 
       P1, L1, U1, R1 := LUDecomposition(A2, method = 'RREF'); 
       E1 := Multiply(Multiply(P1, L1), U1); 
       E2 := MatrixInverse(E1,method = LU); 
       A3 := Multiply(E2, A1); 
       B3 := Multiply(Multiply(E2, N), C2); 
      if 0 < s then A4 := DeleteRow(A3, 1 .. s); 
                    B4 := DeleteRow(B3, 1 .. s) 
      else A4 := A3; 
           B4 := B3 
        end if; 
       A5 := DeleteColumn(A4, r+1 .. n); 
       B5 := DeleteColumn(B4, r+1 .. n); 
        M := A5; N := B5; r := Rank(B5); 
        n := ColumnDimension(B5); 
        m := RowDimension(B5) 
        end do; 
 RETURN(A5, B5); 
end proc:

#----------------------------------------------------------------------
addtohelp(RegularPencil):
#----------------------------------------------------------------------

## HELP RegularPencil
## CALLING SEQUENCE:
##     RegularPencil(A,B)
## DESCRIPTION:
## - Given a pencil matrix A-t*B
## - A, B: Matrices of size (m,n) with m>=n. 
## RETURN: A pair matrix (A',B') such that A'-t*B' is regular pencil matrix. 
## EXAMPLES: M := Matrix(2, 2, {(1, 1) = t^2, (1, 2) = t+1, (2, 1) = 1, (2, 2) = t})
##  A:= FirstCompanionMatrix(M, var)
##  B:= SecondCompanionMatrix(M, var)
## RegularPencil(A,B)
## Matrix(3, 3, {(1, 1) = 0, (1, 2) = 0, (1, 3) = -1, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (3, 1) = 1, (3, 2) = 0, (3, 3) = -1}), 
## Matrix(3, 3, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 1, (3, 1) = 0, (3, 2) = -1, (3, 3) = 0})

RegularPencil := proc (A, B) 
  local m, n, p, r, s, i, j, B1, B2, B3, B4, B5, B6, A1, A2, A3, A4, A5, A6, C, C1, C2, P, L, U, P1, L1, U1, E1, E2, R, R1, E, F, M, N; 
    m := RowDimension(B); 
    n := ColumnDimension(B); 
    r := Rank(B); M := A; N := B; 
   while r < n do 
        B1 := Transpose(N); 
        P, L, U, R := LUDecomposition(B1, method = 'RREF'); 
        C := Multiply(Multiply(P, L), U); 
       C1 := Transpose(C); 
       C2 := MatrixInverse(C1,method = LU); 
       B2 := Multiply(N, C2); 
       A1 := Multiply(M, C2); 
       A2 := DeleteColumn(A1, 1 .. r); 
        s := Rank(A2); 
       P1, L1, U1, R1 := LUDecomposition(A2, method = 'RREF'); 
      E1 := Multiply(Multiply(P1, L1), U1); 
      E2 := MatrixInverse(E1, method = LU); 
      A3 := Multiply(E2, A1); 
      B3 := Multiply(Multiply(E2, N), C2); 
      if 0 < s then 
          A4 := DeleteRow(A3, 1 .. s); 
          B4 := DeleteRow(B3, 1 .. s) 
      else A4 := A3; B4 := B3 
      end if; 
      A5 := DeleteColumn(A4, r+1 .. n); 
      B5 := DeleteColumn(B4, r+1 .. n); 
       M := A5; 
       N := B5; 
       r := Rank(B5); 
       n := ColumnDimension(B5); 
       m := RowDimension(B5) end do; 
      if r = m then RETURN(M, N) 
      else A6 := Transpose(M); 
           B6 := Transpose(N); 
         M, N := ReductionStep(A6, B6) 
       end if; 
  RETURN(M, N);
end proc:
#----------------------------------------------------------------------
addtohelp(ComputationEigenvalues):
#----------------------------------------------------------------------

## HELP ComputationEigenvalues
## CALLING SEQUENCE:
##     ComputationEigenvalues(M,var)
## DESCRIPTION:
## - Given a univariate polynomial matrix M(t) size (m,n) with m<=n.
## - var -: Variable t
## RETURN: List of value t such that rank M(t) drop.
## EXAMPLES:M := Matrix(4, 5, {(1, 1) = t^3, (1, 2) = 0, (1, 3) = t^2-1, (1, 4) = 2*t, (1, 5) = 0, 
##                              (2, 1) = 0, (2, 2) = t^2+1, (2, 3) = 2*t^2, (2, 4) = 0, (2, 5) = 0, 
##                              (3, 1) = 0, (3, 2) = t+1, (3, 3) = 0, (3, 4) = 1, (3, 5) = 1, 
##                              (4, 1) = 0, (4, 2) = 0, (4, 3) = t^2-2, (4, 4) = 0, (4, 5) = 0})
##  ComputationEigenvalues(M,t) 
##  [0, sqrt(2), -sqrt(2), I, -I]
ComputationEigenvalues := proc (M, var) 
    local A, B, C, C1, C2; 
      A := FirstCompanionMatrix(M, var); 
      B := SecondCompanionMatrix(M, var); 
 C1, C2 := RegularPencil(A, B); 
      C := Eigenvalues(C1, C2, output = 'list'); 
 RETURN(C);
end proc:

[__];