Changeset b24184 in git for Singular/LIB


Ignore:
Timestamp:
Jun 16, 1999, 7:24:35 PM (25 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
e38016b4a56783e188de1f4d208aab897b2f6f25
Parents:
64b0be57e2ca9b4b5ca36cb9e9d11a9e6e2c0eb8
Message:
*** empty log message ***


git-svn-id: file:///usr/local/Singular/svn/trunk@3140 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/jordan.lib

    r64b0be5 rb24184  
    11///////////////////////////////////////////////////////////////////////////////
    22
    3 version="$Id: jordan.lib,v 1.10 1999-06-08 09:15:20 mschulze Exp $";
     3version="$Id: jordan.lib,v 1.11 1999-06-16 17:24:35 mschulze Exp $";
    44info="
    55LIBRARY: jordan.lib  PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM
    … …  
    1515///////////////////////////////////////////////////////////////////////////////
    1616
     17static proc countblocks(matrix M)
     18{
     19  int n;
     20
     21  int i=1;
     22  int r,r0;
     23  while(i<=nrows(M))
     24  {
     25    n++;
     26    r=nrows(M[i]);
     27    r0=r;
     28
     29    while(i<r0&&i<nrows(M))
     30    {
     31      i++;
     32      if(i<=nrows(M))
     33      {
     34        r=nrows(M[i]);
     35        if(r>r0)
     36        {
     37          r0=r;
     38        }
     39      }
     40    }
     41
     42    i++;
     43  }
     44
     45  return(n);
     46}
     47///////////////////////////////////////////////////////////////////////////////
     48
     49static proc getblock(matrix M,intvec v)
     50{
     51  matrix M0[size(v)][size(v)]=M[v,v];
     52  return(M0);
     53}
     54///////////////////////////////////////////////////////////////////////////////
     55
    1756proc jordan(matrix M,list #)
    18 "USAGE:   jordan(M[,opt]); with M constant square matrix and opt integer.
     57"USAGE:   jordan(M[,opt]); M constant square matrix, opt integer
    1958ASSUME:  The eigenvalues of M are in the coefficient field.
    20 RETURN:  The procedure returns a list l with 3 entries of type
     59RETURN:  The procedure returns a list jd with 3 entries of type
    2160         ideal, list of intvecs, matrix with
    22          l[1] eigenvalues of M,
    23          l[2][i][j] size of j-th Jordan block with eigenvalue l[1][i], and
    24          l[3]^(-1)*M*l[3] in Jordan normal form.
    25          Depending on opt, only certain entries of l are computed.
    26          If opt=-1, l[1] is computed,
    27          if opt= 0, l[1] and l[2] are computed,
    28          if opt= 1, l[1],l[2], and l[3] are computed, and,
    29          if opt= 2, l[1] and l[3] are computed.
     61         jd[1] eigenvalues of M,
     62         jd[2][i][j] size of j-th Jordan block with eigenvalue jd[1][i], and
     63         jd[3]^(-1)*M*jd[3] in Jordan normal form.
     64         Depending on opt, only certain entries of jd are computed.
     65         If opt=-1, jd[1] is computed,
     66         if opt= 0, jd[1] and jd[2] are computed,
     67         if opt= 1, jd[1],jd[2], and jd[3] are computed, and,
     68         if opt= 2, jd[1] and jd[3] are computed.
    3069         By default, opt=0.
    3170NOTE:    A non constant polynomial matrix M is replaced by its constant part.
    … …  
    3372"
    3473{
    35   // test if M is a square matrix
     74  // check parameter M
    3675  int n=nrows(M);
     76  if(n==0)
     77  {
     78    print("//empty matrix");
     79    return(list());
     80  }
    3781  if(n!=ncols(M))
    3882  {
    … …  
    4488  M=jet(M,0);
    4589
    46   // change to polynomial ring for factorization
    47   def br=basering;
    48   changeord("pr","dp");
    49   matrix M=imap(br,M);
    50 
    51   // factorize characteristic polynomial of M
    52   list l=factorize(det(M-var(1)*freemodule(n)),2);
    53 
    54   // get multiplicities mM of M
    55   def eM,mM=l[1..2];
    56   kill l;
    57 
    58   // test if eigenvalues of M are in the coefficient field
    59   int i;
    60   for(i=size(eM);i>=1;i--)
    61   {
    62     if(deg(eM[i])>1)
    63     {
    64       print("//eigenvalues not in the coefficient field");
    65       setring br;
    66       kill pr;
    67       return(list());
    68     }
    69   }
    70 
    71   // get eigenvalues eM of M
    72   map inv=pr,-var(1);
    73   eM=jet(simplify(inv(eM),1),0);
    74   setring br;
    75   ideal eM=fetch(pr,eM);
    76   kill pr;
    77 
    78   // sort eigenvalues eM
    79   int j;
     90  // try to increase number of blocks by transposing M
     91  if(countblocks(M)<countblocks(transpose(M)))
     92  {
     93    M=transpose(M);
     94  }
     95
     96  // compute eigenvalues eM and multiplicities mM of M
     97  list fac;
     98  matrix M0;
     99  ideal eM,eM0;
     100  intvec mM,mM0;
     101  intvec u;
     102  int i=1;
     103  int j,r,r0;
     104  while(i<=nrows(M))
     105  {
     106    u=i;
     107    r=nrows(M[i]);
     108    r0=r;
     109
     110    // find next block
     111    while(i<r0&&i<nrows(M))
     112    {
     113      i++;
     114      if(i<=nrows(M))
     115      {
     116        u=u,i;
     117        r=nrows(M[i]);
     118        if(r>r0)
     119        {
     120          r0=r;
     121        }
     122      }
     123    }
     124
     125    // if 1x1-block
     126    if(size(u)==1)
     127    {
     128      if(mM[1]==0)
     129      {
     130        eM=M[u[1]][u[1]];
     131        mM=1;
     132      }
     133      else
     134      {
     135        eM=eM,ideal(M[u[1]][u[1]]);
     136        mM=mM,1;
     137      }
     138    }
     139    else
     140    {
     141      // factorize characteristic polynomial of block M0
     142      M0=getblock(M,u);
     143      fac=factorize(det(M0-var(1)*freemodule(size(u))),2);
     144      eM0,mM0=fac[1..2];
     145
     146      // compute eigenvalues eM0 of M0
     147      if(1<var(1))
     148      {
     149        for(j=ncols(eM0);j>=1;j--)
     150        {
     151          if(deg(eM0[j])>1)
     152          {
     153            print("//eigenvalues not in the coefficient field");
     154            return(list());
     155          }
     156          eM0[j]=-(eM0[j][2]/var(1))/eM0[j][1];
     157        }
     158      }
     159      else
     160      {
     161        for(j=ncols(eM0);j>=1;j--)
     162        {
     163          if(deg(eM0[j])>1)
     164          {
     165            print("//eigenvalues not in the coefficient field");
     166            return(list());
     167          }
     168          eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1));
     169        }
     170      }
     171
     172      if(mM[1]==0)
     173      {
     174        eM=eM0;
     175        mM=mM0;
     176      }
     177      else
     178      {
     179        eM=eM,eM0;
     180        mM=mM,mM0;
     181      }
     182    }
     183
     184    i++;
     185  }
     186
     187  // sort eigenvalues
    80188  poly e;
    81189  int m;
    82   for(i=size(eM);i>=2;i--)
     190  for(i=ncols(eM);i>=2;i--)
    83191  {
    84192    for(j=i-1;j>=1;j--)
    85193    {
    86       if(eM[i]<eM[j])
     194     if(eM[i]<eM[j])
    87195      {
    88196        e=eM[i];
    … …  
    95203    }
    96204  }
    97   kill e,m;
     205
     206  // remove multiple eigenvalues
     207  i=1;
     208  j=2;
     209  while(j<=ncols(eM))
     210  {
     211    if(eM[i]==eM[j])
     212    {
     213      mM[i]=mM[i]+mM[j];
     214    }
     215    else
     216    {
     217      i++;
     218      eM[i]=eM[j];
     219    }
     220    j++;
     221  }
     222  eM=eM[1..i];
     223  mM=mM[1..i];
    98224
    99225  // get optional parameter opt
    … …  
    128254  }
    129255
    130   // do for all eigenvalues eM[i]
    131256  for(i=ncols(eM);i>=1;i--)
    132257  {
    … …  
    183308    }
    184309  }
    185   kill l;
    186 
    187   // create return list l
    188   list l=eM;
     310
     311  list jd=eM;
    189312  if(opt<=1)
    190313  {
    191     l[2]=bM;
     314    jd[2]=bM;
    192315  }
    193316  if(opt>=1)
    194317  {
    195     l[3]=V;
    196   }
    197 
    198   return(l);
     318    jd[3]=V;
     319  }
     320  return(jd);
    199321}
    200322example
    … …  
    207329///////////////////////////////////////////////////////////////////////////////
    208330
    209 proc jordanmatrix(list l)
    210 "USAGE:   jordanmatrix(l); with l list of ideal and list of intvecs.
    211 RETURN:  The procedure returns the Jordan matrix J with eigenvalues l[1] and
    212          size l[2][i][j] of j-th Jordan block with eigenvalue l[1][i].
     331proc jordanmatrix(list jd)
     332"USAGE:   jordanmatrix(jd); jd list of ideal and list of intvecs
     333RETURN:  The procedure returns the Jordan matrix J with eigenvalues jd[1] and
     334         size jd[2][i][j] of j-th Jordan block with eigenvalue jd[1][i].
    213335EXAMPLE: example jordanmatrix; shows an example.
    214336"
    215337{
    216338  // check parameters
    217   if(size(l)<2)
     339  if(size(jd)<2)
    218340  {
    219341    print("//not enough entries in argument list");
    … …  
    221343    return(J);
    222344  }
    223   def eJ,bJ=l[1..2];
    224   kill l;
     345  def eJ,bJ=jd[1..2];
    225346  if(typeof(eJ)!="ideal")
    226347  {
    … …  
    302423
    303424proc jordanform(matrix M)
    304 "USAGE:   jordanform(M); with M constant square matrix.
     425"USAGE:   jordanform(M); M constant square matrix
    305426ASSUME:  The eigenvalues of M are in the coefficient field.
    306427RETURN:  The procedure returns the Jordan normal form of M.
    … …  
    321442
    322443proc invmat(matrix M)
    323 "USAGE:   invmat(M); with M constant square matrix.
     444"USAGE:   invmat(M); M constant square matrix
    324445ASSUME:  M is invertible.
    325446RETURN:  The procedure returns the inverse matrix of M.
Note: See TracChangeset for help on using the changeset viewer.