source: git/Singular/LIB/jordan.lib @ 371848

///////////////////////////////////////////////////////////////////////////////
version="$Id: jordan.lib,v 1.20 2001-01-02 14:33:25 greuel Exp $";
category="Linear Algebra";
info="
LIBRARY: jordan.lib  Jordan Normal Form of a Matrix with rational Eigenvalues
AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de

PROCEDURES:
 jordan(M[,opt]);  eigenvalues, Jordan block sizes, Jordan transformation of M
 jordanmatrix(l);  Jordan matrix with eigenvalues l[1], Jordan block sizes l[2]
 jordanform(M);    Jordan normal form of constant square matrix M
 invmat(M);        inverse matrix of invertible constant matrix M
";

LIB "ring.lib";
///////////////////////////////////////////////////////////////////////////////

static proc countblocks(matrix M)
{
  int b,r,r0;

  int i=1;
  while(i<=nrows(M))
  {
    b++;
    r=nrows(M[i]);
    r0=r;

    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    while(i<r0&&i<nrows(M))
    {
      i++;
      if(i<=nrows(M))
      {
        r=nrows(M[i]);
        if(r>r0)
        {
          r0=r;
        }
      }
    }
    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");

    i++;
  }

  return(b);
}
///////////////////////////////////////////////////////////////////////////////

static proc getblock(matrix M,intvec v)
{
  matrix M0[size(v)][size(v)]=M[v,v];
  return(M0);
}
///////////////////////////////////////////////////////////////////////////////

proc jordan(matrix M,list #)
"USAGE:   jordan(M[,opt]); M constant square matrix, opt integer
ASSUME:  The eigenvalues of M are in the coefficient field.
RETURN:  The procedure returns a list jd with 3 entries:
@format       
         - jd[1] ideal, eigenvalues of M,
         - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with 
           eigenvalue jd[1][i], and
         - jd[3] a matrix, jd[3]^(-1)*M*jd[3] in Jordan normal form.
         Depending on opt, only certain entries of jd are computed.
           If opt=-1, jd[1] is computed,
           if opt= 0, jd[1] and jd[2] are computed,
           if opt= 1, jd[1], jd[2], and jd[3] are computed, and,
           if opt= 2, jd[1] and jd[3] are computed.
         By default, opt=0.
@end format
NOTE:    A non constant polynomial matrix M is replaced by its constant part.
DISPLAY: The procedure displays comments if printlevel>=1.
EXAMPLE: example jordan; shows an example.
"
{
  int n=nrows(M);
  if(n==0)
  {
    print("//empty matrix");
    return(list());
  }
  if(n!=ncols(M))
  {
    print("//no square matrix");
    return(list());
  }

  M=jet(M,0);

  dbprint(printlevel-voice+2,"//counting blocks of matrix...");
  int i=countblocks(M);
  dbprint(printlevel-voice+2,"//...blocks of matrix counted");
  if(i==1)
  {
    dbprint(printlevel-voice+2,"//matrix has 1 block");
  }
  else
  {
    dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks");
  }

  dbprint(printlevel-voice+2,"//counting blocks of transposed matrix...");
  int j=countblocks(transpose(M));
  dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted");
  if(j==1)
  {
    dbprint(printlevel-voice+2,"//transposed matrix has 1 block");
  }
  else
  {
    dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks");
  }

  if(i<j)
  {
    dbprint(printlevel-voice+2,"//transposing matrix...");
    M=transpose(M);
    dbprint(printlevel-voice+2,"//...matrix transposed");
  }

  list fd;
  matrix M0;
  poly cp;
  ideal eM,eM0;
  intvec mM,mM0;
  intvec u;
  int b,r,r0;

  i=1;
  while(i<=nrows(M))
  {
    b++;
    u=i;
    r=nrows(M[i]);
    r0=r;

    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    while(i<r0&&i<nrows(M))
    {
      i++;
      if(i<=nrows(M))
      {
        u=u,i;
        r=nrows(M[i]);
        if(r>r0)
        {
          r0=r;
        }
      }
    }
    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");

    if(size(u)==1)
    {
      dbprint(printlevel-voice+2,"//1x1-block:");
      dbprint(printlevel-voice+2,M[u[1]][u[1]]);

      if(mM[1]==0)
      {
        eM=M[u[1]][u[1]];
        mM=1;
      }
      else
      {
        eM=eM,ideal(M[u[1]][u[1]]);
        mM=mM,1;
      }
    }
    else
    {
      dbprint(printlevel-voice+2,
        "//"+string(size(u))+"x"+string(size(u))+"-block:");
      M0=getblock(M,u);
      dbprint(printlevel-voice+2,M0);

      dbprint(printlevel-voice+2,"//characteristic polynomial:");
      cp=det(module(M0-var(1)*freemodule(size(u))));
      dbprint(printlevel-voice+2,cp);

      dbprint(printlevel-voice+2,"//factorizing characteristic polynomial...");
      fd=factorize(cp,2);
      dbprint(printlevel-voice+2,"//...characteristic polynomial factorized");

      dbprint(printlevel-voice+2,"//computing eigenvalues...");
      eM0,mM0=fd[1..2];
      if(1<var(1))
      {
        for(j=ncols(eM0);j>=1;j--)
        {
          if(deg(eM0[j])>1)
          {
            print("//eigenvalues not in the coefficient field");
            return(list());
          }
          if(eM0[j][1]==0)
          {
            eM0[j]=0;
          }
          else
          {
            eM0[j]=-eM0[j][2]/(eM0[j][1]/var(1));
          }
        }
      }
      else
      {
        for(j=ncols(eM0);j>=1;j--)
        {
          if(deg(eM0[j])>1)
          {
            print("//eigenvalues not in the coefficient field");
            return(list());
          }
          if(eM0[j][2]==0)
          {
            eM0[j]=0;
          }
          else
          {
            eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1));
          }
        }
      }
      dbprint(printlevel-voice+2,"//...eigenvalues computed");

      if(mM[1]==0)
      {
        eM=eM0;
        mM=mM0;
      }
      else
      {
        eM=eM,eM0;
        mM=mM,mM0;
      }
    }

    i++;
  }

  dbprint(printlevel-voice+2,"//sorting eigenvalues...");
  poly e;
  int m;
  for(i=ncols(eM);i>=2;i--)
  {
    for(j=i-1;j>=1;j--)
    {
     if(eM[i]<eM[j])
      {
        e=eM[i];
        eM[i]=eM[j];
        eM[j]=e;
        m=mM[i];
        mM[i]=mM[j];
        mM[j]=m;
      }
    }
  }
  dbprint(printlevel-voice+2,"//...eigenvalues sorted");

  dbprint(printlevel-voice+2,"//removing multiple eigenvalues...");
  i=1;
  j=2;
  while(j<=ncols(eM))
  {
    if(eM[i]==eM[j])
    {
      mM[i]=mM[i]+mM[j];
    }
    else
    {
      i++;
      eM[i]=eM[j];
      mM[i]=mM[j];
    }
    j++;
  }
  eM=eM[1..i];
  mM=mM[1..i];
  dbprint(printlevel-voice+2,"//...multiple eigenvalues removed");

  dbprint(printlevel-voice+2,"//eigenvalues:");
  dbprint(printlevel-voice+2,eM);
  dbprint(printlevel-voice+2,"//multiplicities:");
  dbprint(printlevel-voice+2,mM);

  int opt=0;
  if(size(#)>0)
  {
    if(typeof(#[1])=="int")
    {
      opt=#[1];
    }
  }
  if(opt<0)
  {
    return(list(eM));
  }
  int k,l;
  matrix I=freemodule(n);
  matrix Mi,Ni;
  module sNi;
  list K;
  if(opt>=1)
  {
    module V,K1,K2;
    matrix v[n][1];
  }
  if(opt<=1)
  {
    list bM;
    intvec bMi;
  }

  for(i=ncols(eM);i>=1;i--)
  {
    Mi=M-eM[i]*I;

    dbprint(printlevel-voice+2,
      "//computing kernels of powers of matrix minus eigenvalue "
      +string(eM[i]));
    K=list(module());
    for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
    {
      sNi=syz(Ni);
      K=K+list(sNi);
    }
    dbprint(printlevel-voice+2,"//...kernels computed");

    if(opt<=1)
    {
      dbprint(printlevel-voice+2,
        "//computing Jordan block sizes for eigenvalue "
        +string(eM[i])+"...");
      bMi=0;
      bMi[size(K[2])]=0;
      for(j=size(K);j>=2;j--)
      {
        for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--)
        {
          bMi[k]=bMi[k]+1;
        }
      }
      bM=list(bMi)+bM;
      dbprint(printlevel-voice+2,"//...Jordan block sizes computed");
    }

    if(opt>=1)
    {
      dbprint(printlevel-voice+2,
        "//computing Jordan basis vectors for eigenvalue "
        +string(eM[i])+"...");
      if(size(K)>1)
      {
        for(j,K1=2,0;j<=size(K)-1;j++)
        {
          K2=K[j];
          K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1]))));
          K1=K2;
        }
        K[j]=interred(reduce(K[j],std(K1)));
      }
      for(j=size(K);j>=2;j--)
      {
        for(k=size(K[j]);k>=1;k--)
        {
          v=K[j][k];
          for(l=j;l>=1;l--)
          {
            V=module(v)+V;
            v=Mi*v;
          }
        }
      }
      dbprint(printlevel-voice+2,"//...Jordan basis vectors computed");
    }
  }

  list jd=eM;
  if(opt<=1)
  {
    jd[2]=bM;
  }
  if(opt>=1)
  {
    jd[3]=V;
  }
  return(jd);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  jordan(M);
}
///////////////////////////////////////////////////////////////////////////////

proc jordanmatrix(list jd)
"USAGE:   jordanmatrix(jd); jd list of ideal and list of intvecs
RETURN:  The procedure returns the Jordan matrix J, satisfying:
@*       - eigenvalues of J are given by jd[1]
@*       - j-th Jordan block of J with eigenvalue jd[1][i] has size jd[2][i][j]
DISPLAY: The procedure displays comments if printlevel>=1.
EXAMPLE: example jordanmatrix; shows an example.
"
{
  if(size(jd)<2)
  {
    print("//not enough entries in argument list");
    matrix J[1][0];
    return(J);
  }
  def eJ,bJ=jd[1..2];
  if(typeof(eJ)!="ideal")
  {
    print("//first entry in argument list not an ideal");
    matrix J[1][0];
    return(J);
  }
  if(typeof(bJ)!="list")
  {
    print("//second entry in argument list not a list");
    matrix J[1][0];
    return(J);
  }
  if(size(eJ)<size(bJ))
  {
    int s=size(eJ);
  }
  else
  {
    int s=size(bJ);
  }

  int i,j,k,n;
  for(i=s;i>=1;i--)
  {
    if(typeof(bJ[i])!="intvec")
    {
      print("//second entry in argument list not a list of intvecs");
      matrix J[1][0];
      return(J);
    }
    else
    {
      for(j=size(bJ[i]);j>=1;j--)
      {
        k=bJ[i][j];
        if(k>0)
        {
          n=n+k;
        }
      }
    }
  }

  int l;
  matrix J[n][n];
  for(i,l=1,1;i<=s;i++)
  {
    for(j=1;j<=size(bJ[i]);j++)
    {
      k=bJ[i][j];
      if(k>0)
      {
        while(k>=2)
        {
          J[l,l]=eJ[i];
          J[l,l+1]=1;
          k,l=k-1,l+1;
        }
        J[l,l]=eJ[i];
        l++;
      }
    }
  }

  return(J);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  list l;
  l[1]=ideal(2,3);
  l[2]=list(intvec(1),intvec(2));
  print(jordanmatrix(l));
}
///////////////////////////////////////////////////////////////////////////////

proc jordanform(matrix M)
"USAGE:   jordanform(M); M constant square matrix
ASSUME:  The eigenvalues of M are in the coefficient field.
RETURN:  matrix, the Jordan normal form of M.
NOTE:    A non constant polynomial matrix M is replaced by its constant part.
DISPLAY: The procedure displays more comments for higher printlevel.
EXAMPLE: example jordanform; shows an example.
"
{
  return(jordanmatrix(jordan(M)));
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  print(jordanform(M));
}
///////////////////////////////////////////////////////////////////////////////

proc invmat(matrix M)
"USAGE:   invmat(M); M constant square matrix
ASSUME:  M is invertible.
RETURN:  The procedure returns the inverse matrix of M.
NOTE:    A non constant polynomial matrix M is replaced by its constant part.
EXAMPLE: example invmat; shows an example.
"
{
  if(nrows(M)==ncols(M))
  {
    matrix invM=lift(jet(M,0),freemodule(nrows(M)));
  }
  else
  {
    ERROR("// ** no square matrix");
    //matrix invM[1][0];
  }
  return(invM);
}
example
{ "EXAMPLE:"; echo=2;
  ring R=0,x,dp;
  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
  print(M);
  print(invmat(M));
}
///////////////////////////////////////////////////////////////////////////////