Changeset 092430 in git


Ignore:
Timestamp:
Jun 21, 1999, 10:45:41 AM (24 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
0f0974e0182296258c671cafb367b12244db65c1
Parents:
d694de170d576d5f00d6f1047b53246e20433998
Message:
*** empty log message ***


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/jordan.lib

    rd694de r092430  
    11///////////////////////////////////////////////////////////////////////////////
    22
    3 version="$Id: jordan.lib,v 1.11 1999-06-16 17:24:35 mschulze Exp $";
     3version="$Id: jordan.lib,v 1.12 1999-06-21 08:45:41 mschulze Exp $";
    44info="
    55LIBRARY: jordan.lib  PROCEDURES TO COMPUTE THE JORDAN NORMAL FORM
    … …  
    1717static proc countblocks(matrix M)
    1818{
    19   int n;
     19  int b,r,r0;
    2020
    2121  int i=1;
    22   int r,r0;
    2322  while(i<=nrows(M))
    2423  {
    25     n++;
     24    b++;
    2625    r=nrows(M[i]);
    2726    r0=r;
    2827
     28    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    2929    while(i<r0&&i<nrows(M))
    3030    {
    … …  
    3939      }
    4040    }
     41    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
    4142
    4243    i++;
    4344  }
    4445
    45   return(n);
     46  return(b);
    4647}
    4748///////////////////////////////////////////////////////////////////////////////
    … …  
    7273"
    7374{
    74   // check parameter M
    7575  int n=nrows(M);
    7676  if(n==0)
    … …  
    8585  }
    8686
    87   // replace M by its constant part
    8887  M=jet(M,0);
    8988
    90   // try to increase number of blocks by transposing M
    91   if(countblocks(M)<countblocks(transpose(M)))
    92   {
     89  dbprint(printlevel-voice+2,"//counting blocks of matrix...");
     90  int i=countblocks(M);
     91  dbprint(printlevel-voice+2,"//...blocks of matrix counted");
     92  if(i==1)
     93  {
     94    dbprint(printlevel-voice+2,"//matrix has 1 block");
     95  }
     96  else
     97  {
     98    dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks");
     99  }
     100
     101  dbprint(printlevel-voice+2,"//counting blocks of transposed matrix...");
     102  int j=countblocks(transpose(M));
     103  dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted");
     104  if(j==1)
     105  {
     106    dbprint(printlevel-voice+2,"//transposed matrix has 1 block");
     107  }
     108  else
     109  {
     110    dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks");
     111  }
     112
     113  if(i<j)
     114  {
     115    dbprint(printlevel-voice+2,"//transposing matrix...");
    93116    M=transpose(M);
    94   }
    95 
    96   // compute eigenvalues eM and multiplicities mM of M
    97   list fac;
     117    dbprint(printlevel-voice+2,"//...matrix transposed");
     118  }
     119
     120  list fd;
    98121  matrix M0;
     122  poly cp;
    99123  ideal eM,eM0;
    100124  intvec mM,mM0;
    101125  intvec u;
    102   int i=1;
    103   int j,r,r0;
     126  int b,r,r0;
     127
     128  i=1;
    104129  while(i<=nrows(M))
    105130  {
     131    b++;
    106132    u=i;
    107133    r=nrows(M[i]);
    108134    r0=r;
    109135
    110     // find next block
     136    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    111137    while(i<r0&&i<nrows(M))
    112138    {
    … …  
    122148      }
    123149    }
    124 
    125     // if 1x1-block
     150    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
     151
    126152    if(size(u)==1)
    127153    {
     154      dbprint(printlevel-voice+2,"//1x1-block:");
     155      dbprint(printlevel-voice+2,M[u[1]][u[1]]);
     156
    128157      if(mM[1]==0)
    129158      {
    … …  
    139168    else
    140169    {
    141       // factorize characteristic polynomial of block M0
     170      dbprint(printlevel-voice+2,
     171        "//"+string(size(u))+"x"+string(size(u))+"-block:");
    142172      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
     173      dbprint(printlevel-voice+2,M0);
     174
     175      dbprint(printlevel-voice+2,"//characteristic polynomial:");
     176      cp=det(module(M0-var(1)*freemodule(size(u))));
     177      dbprint(printlevel-voice+2,cp);
     178
     179      dbprint(printlevel-voice+2,"//factorizing characteristic polynomial...");
     180      fd=factorize(cp,2);
     181      dbprint(printlevel-voice+2,"//...characteristic polynomial factorized");
     182
     183      dbprint(printlevel-voice+2,"//computing eigenvalues...");
     184      eM0,mM0=fd[1..2];
    147185      if(1<var(1))
    148186      {
    … …  
    169207        }
    170208      }
     209      dbprint(printlevel-voice+2,"//...eigenvalues computed");
    171210
    172211      if(mM[1]==0)
    … …  
    185224  }
    186225
    187   // sort eigenvalues
     226  dbprint(printlevel-voice+2,"//sorting eigenvalues...");
    188227  poly e;
    189228  int m;
    … …  
    203242    }
    204243  }
    205 
    206   // remove multiple eigenvalues
     244  dbprint(printlevel-voice+2,"//...eigenvalues sorted");
     245
     246  dbprint(printlevel-voice+2,"//removing multiple eigenvalues...");
    207247  i=1;
    208248  j=2;
    … …  
    217257      i++;
    218258      eM[i]=eM[j];
     259      mM[i]=mM[j];
    219260    }
    220261    j++;
    … …  
    222263  eM=eM[1..i];
    223264  mM=mM[1..i];
    224 
    225   // get optional parameter opt
     265  dbprint(printlevel-voice+2,"//...multiple eigenvalues removed");
     266
     267  dbprint(printlevel-voice+2,"//eigenvalues:");
     268  dbprint(printlevel-voice+2,eM);
     269  dbprint(printlevel-voice+2,"//multiplicities:");
     270  dbprint(printlevel-voice+2,mM);
     271
    226272  int opt=0;
    227273  if(size(#)>0)
    … …  
    236282    return(list(eM));
    237283  }
    238 
    239   // define variables
    240284  int k,l;
    241   matrix E=freemodule(n);
     285  matrix I=freemodule(n);
    242286  matrix Mi,Ni;
    243287  module sNi;
    … …  
    256300  for(i=ncols(eM);i>=1;i--)
    257301  {
    258     Mi=M-eM[i]*E;
    259 
    260     // compute kernels K of powers of Mi
     302    Mi=M-eM[i]*I;
     303
     304    dbprint(printlevel-voice+2,
     305      "//computing kernels of powers of matrix minus eigenvalue "
     306      +string(eM[i]));
    261307    K=list(module());
    262308    for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
    … …  
    265311      K=K+list(sNi);
    266312    }
     313    dbprint(printlevel-voice+2,"//...kernels computed");
    267314
    268315    if(opt<=1)
    269316    {
    270       // compute Jordan block size vector bMi corresponding to eigenvalue eM[i]
     317      dbprint(printlevel-voice+2,
     318        "//computing Jordan block sizes for eigenvalue "
     319        +string(eM[i])+"...");
    271320      bMi=0;
    272321      bMi[size(K[2])]=0;
    … …  
    279328      }
    280329      bM=list(bMi)+bM;
     330      dbprint(printlevel-voice+2,"//...Jordan block sizes computed");
    281331    }
    282332
    283333    if(opt>=1)
    284334    {
    285       // compute Jordan basis vectors K corresponding to eigenvalue eM[i]
     335      dbprint(printlevel-voice+2,
     336        "//computing Jordan basis vectors for eigenvalue "
     337        +string(eM[i])+"...");
    286338      if(size(K)>1)
    287339      {
    … …  
    306358        }
    307359      }
     360      dbprint(printlevel-voice+2,"//...Jordan basis vectors computed");
    308361    }
    309362  }
    … …  
    336389"
    337390{
    338   // check parameters
    339391  if(size(jd)<2)
    340392  {
    … …  
    365417  }
    366418
    367   // get size of Jordan matrix
    368419  int i,j,k,n;
    369420  for(i=s;i>=1;i--)
    … …  
    388439  }
    389440
    390   // create Jordan matrix
    391441  int l;
    392442  matrix J[n][n];
Note: See TracChangeset for help on using the changeset viewer.