Changeset 6188357 in git


Ignore:
Timestamp:
Jan 3, 2001, 2:24:49 AM (23 years ago)
Author:
Gert-Martin Greuel <greuel@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
9cc5d37df563ce9febe8a449f39cf9986d705e1e
Parents:
cda1ad905ad654f1520fe50c63f3628c521558b0
Message:
* GMG: charpoly umbenannt in charpoly_B, neue proc charpoly zugefuegt
*     mit jordan.lib vereinigt. Proc invmat entfernt.
*     Kosmetik,


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/linalg.lib

    rcda1ad r6188357  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.7 2000-12-22 14:12:32 greuel Exp $";
     3version="$Id: linalg.lib,v 1.8 2001-01-03 01:24:49 greuel Exp $";
    44category="Linear Algebra";
    55info="
    66LIBRARY:  linalg.lib    Algorithmic Linear Algebra
    7 AUTHOR:   Ivor Saynisch (ivs@math.tu-cottbus.de)
     7AUTHORS:  Ivor Saynisch (ivs@math.tu-cottbus.de)
     8@*        Mathias Schulze (mschulze@mathematik.uni-kl.de)
    89
    910PROCEDURES:
     
    1718 charpoly(A,v);     characteristic polynomial of A ( using busadj(A) )
    1819 adjoint(A);        adjoint of A ( using busadj(A) )
    19  det_B(A);          determinant of A. ( uses busadj(A) )
     20 det_B(A);          determinant of A ( using busadj(A) )
    2021 gaussred(A);       gaussian reduction: P*A=U*S, S a row reduced form of A
    2122 gaussred_pivot(A); gaussian reduction: P*A=U*S, uses row pivoting
    2223 gauss_nf(A);       gaussian normal form of A
    2324 mat_rk(A);         rank of constant matrix A
    24  U_D_O(A);          P*A=U*D*O, P,D,U,O = permutaion,diag,lower-,upper-triang
     25 U_D_O(A);          P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
    2526 pos_def(A,i);      test symmetric matrix for positive definiteness
     27 jordan(M[,opt]);   eigenvalues, Jordan block sizes, transformation matrix
     28 jordanmatrix(l);   Jordan matrix with eigenvalues, Jordan block sizes
     29 jordanform(M);     Jordan normal form of constant square matrix M
    2630";
    2731
    2832LIB "matrix.lib";
    2933LIB "ring.lib";
    30 
     34LIB "elim.lib";
    3135//////////////////////////////////////////////////////////////////////////////
    3236// help functions
     
    3943}
    4044
    41 static
    42 proc const_mat(matrix A)
     45static proc const_mat(matrix A)
    4346"RETURN:   1 (0) if A is (is not) a constant matrix"
    4447{
     
    99102/////////////////////////////////////////////////////////////////////////////
    100103
    101 proc inverse(matrix A)
     104proc inverse(matrix A, list #)
    102105"USAGE:    inverse(A);  A = constant, square matrix
    103106RETURN:   matrix, the inverse matrix of A
    104 NOTE:     parameters and minpoly are allowed, uses std applied to A enlarged
    105           bei unit matrix
     107NOTE:     parameters and minpoly are allowed. The procedure applies std with
     108          ordering (c,dp) to A-enlarged-by-unit-matrix
    106109EXAMPLE:  example inverse; shows an example"
    107110{
     
    113116
    114117  if (ncols(A)!=n){
    115     "input is not a square matrix";
     118    "// ** input is not a square matrix";
    116119    matrix invmat;
    117120    return(invmat);
    118121  }
    119122  if(!const_mat(A)){
    120     "input is not a constant matrix";
     123    "// ** input is not a constant matrix";
    121124    matrix invmat;
    122125    return(invmat);
     
    130133  matrix B=transpose(concat(A,unitmat(n)));
    131134  matrix D=transpose(std(B));
     135
     136//option(intStrategy);
     137//option(nointStrategy);
     138//if (#[1]==1) {matrix D=transpose(interred(B));}
     139
    132140  matrix D1=submat(D,1..n,1..n);
    133141  matrix D2=submat(D,1..n,(n+1)..2*n);
     
    146154  for (i=1;i<=n;i=i+1){
    147155    if (deg(B[n+i,i])<0){
    148       "matrix is not invertible";
     156      "// ** matrix is not invertible";
    149157      setring BR;
    150158      kill @R;kill @@R;
     
    181189
    182190  if (ncols(A)!=n){
    183     "input is not a square matrix";;
     191    "// ** input is not a square matrix";;
    184192    return(A);
    185193  }
    186194 
    187195  if(!const_mat(A)){
    188     "input is not a constant matrix";
     196    "// ** input is not a constant matrix";
    189197    return(A);
    190198  }
    191199 
    192200  if(deg(std(A-transpose(A))[1])!=-1){
    193     "input is not a symmetric matrix";
     201    "// ** input is not a symmetric matrix";
    194202    return(A);
    195203  }
     
    222230
    223231  if(!const_mat(A)){
    224     "input is not a constant matrix";
     232    "// ** input is not a constant matrix";
    225233    matrix B;
    226234    return(B);
     
    232240    for(j=1;j<i;j=j+1) {
    233241      k=inner_product(B[j],B[j]);
    234       if (k==0) { "Error: vector with zero length"; return(matrix(B)); }
     242      if (k==0) { "Error: vector of length zero"; return(matrix(B)); }
    235243      B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j];
    236244    }
     
    250258proc diag_test(matrix A)
    251259"USAGE:   diag_test(A); A = const square matrix
    252 NOTE:     Der Test ist nur fuer zerfallende Matrizen aussagefaehig.
    253           Parameterringe werden nicht unterstuetzt (benutzt factorize,gcd).
    254 RETURN:   int,  1 falls A diagonalisierbar, 0 nicht diagonalisierbar
    255                -1 keine Aussage moeglich, da A nicht zerfallend.
     260RETURN:   int,  1 if A is diagonalisable, 0 if not
     261               -1 no statement possible, since A does not split.
     262NOTE:     The test works only for split matrices, i.e if eigenvalues of A
     263          are in the ground field.
     264          Does not work with parameters (uses factorize,gcd).
    256265EXAMPLE:  example diag_test; shows an example"
    257266{
     
    305314  for(i=1;i<=nf;i++){
    306315    if(lead(Xfactors[1,i])>=@t^2){
    307       //"Die Matrix ist nicht zerfallend";
     316      //" matrix does not split";
    308317      setring BR;
    309318      return(-1);
     
    322331   
    323332      if(0!=reduce(b,g)){
    324         //"Die Matrix ist nicht diagonalisierbar!";
     333        //" matrix not diagonalizable";
    325334        setring BR;
    326335        return(0);
     
    346355//////////////////////////////////////////////////////////////////////////////
    347356proc busadj(matrix A)
    348 "USAGE:   busadj(A);  A = square matrix
    349 RETURN:  list L;
    350          L[1] contains the (n+1) coefficients of the characteristic polynomial
    351          X of A, i.e.
     357"USAGE:   busadj(A);  A = square matrix (of size nxn)
     358RETURN:  list L:
     359@format
     360         L[1] contains the (n+1) coefficients of the characteristic
     361              polynomial X of A, i.e.
    352362              X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n
    353363         L[2] contains the n (nxn)-matrices Hk which are the coefficients of
    354               the busadjoint bA=adjoint(E*t-A) of A, i.e.
     364              the busadjoint bA = adjoint(E*t-A) of A, i.e.
    355365              bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0,  ( Hk=L[2][k+1] )
     366@end format
    356367EXAMPLE: example busadj; shows an example"
    357368{
     
    400411//////////////////////////////////////////////////////////////////////////////
    401412proc charpoly(matrix A, list #)
    402 "USAGE:   charpoly(A,[v]); A = square matrix, v = name of a ringvariable
     413"USAGE:   charpoly(A[,v]); A square matrix, v string, name of a variable
     414RETURN:  poly, the characteristic polynomial det(E*v-A)
     415         (default: v=name of last variable)
     416NOTE:    A must be independent of the variable v. The computation uses det.
     417         If printlevel>0, det(E*v-A) is diplayed recursively.
     418EXAMPLE: example charpoly; shows an example"
     419{
     420  int n = nrows(A);
     421  int z = nvars(basering);
     422  int i,j;
     423  string v;
     424  poly X;
     425  if(ncols(A) != n)
     426  {     
     427    "// input is not a square matrix";
     428    return(X);
     429  }     
     430  //---------------------- test for correct variable -------------------------
     431  if( size(#)==0 ){
     432    #[1] = varstr(z);
     433  }
     434  if( typeof(#[1]) == "string") { v = #[1]; }
     435  else
     436  {
     437    "// 2nd argument must be a name of a variable not contained in the matrix";
     438    return(X);
     439  } 
     440  j=-1;
     441  for(i=1; i<=z; i++)
     442  {
     443    if(varstr(i)==v){j=i;}
     444  }
     445  if(j==-1)
     446  {
     447    "// "+v+" is not a variable in the basering";
     448    return(X);
     449  }
     450  if ( size(select1(module(A),j)) != 0 )
     451  {
     452    "// matrix must not contain the variable "+v;
     453    "// change to a ring with an extra variable, if necessary."
     454    return(X);
     455  }
     456  //-------------------------- compute charpoly ------------------------------
     457  X = det(var(j)*unitmat(n)-A);
     458  if( printlevel-voice+2 >0) { showrecursive(X,var(j));}
     459  return(X);
     460}
     461example
     462{ "EXAMPLE"; echo=2;
     463  ring r=0,(x,t),dp;
     464  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
     465  print(A);
     466  charpoly(A,"t");
     467}
     468
     469//////////////////////////////////////////////////////////////////////////////
     470proc charpoly_B(matrix A, list #)
     471"USAGE:   charpoly_B(A[,v]); A square matrix, v string, name of a variable
     472RETURN:  poly, the characteristic polynomial det(E*v-A)
     473         (default: v=name of last variable)
    403474NOTE:    A must be constant in the variable v. The computation uses busadj(A).
    404 RETURN:  poly, the characteristic polynomial det(E*v-A)
    405 EXAMPLE: example charpoly; shows an example"
     475EXAMPLE: example charpoly_B; shows an example"
    406476{
    407477  int i,j;
     
    420490  //test for correct variable
    421491  if( size(#)==0 ){
    422     #[1] = varstr(1);
     492    #[1] = varstr(nvars(BR));
    423493  }
    424494  if( typeof(#[1]) == "string"){
     
    453523  for(i=1; i<=n; i++){
    454524    if(deg(lead(A)[i])>=1){
    455       "matrix must not contain the variable "+v;
     525      "// matrix must not contain the variable "+v;
    456526      kill @R;
    457527      setring BR;
     
    475545  matrix A[3][3]=1,x2,x,x2,6,4,x,4,1;
    476546  print(A);
    477   charpoly(A,"t");
     547  charpoly_B(A,"t");
    478548}
    479549
     
    481551proc adjoint(matrix A)
    482552"USAGE:    adjoint(A);  A = square matrix
     553RETURN:   adjoint matrix of A, i.e. Adj*A=det(A)*E
    483554NOTE:     computation uses busadj(A)
    484 RETURN:   adjoint
    485555EXAMPLE:  example adjoint; shows an example"
    486556{
     
    511581//////////////////////////////////////////////////////////////////////////////
    512582proc inverse_B(matrix A)
    513 "USAGE:     inverse_B(A);  A = square matrix
    514 RETURN:    list Inv with Inv[1]=matrix I and Inv[2]=poly p
    515            and I*A = unitmat(n)*p;
    516 NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
    517 EXAMPLE:   example inverse_B; shows an example"
     583"USAGE:    inverse_B(A);  A = square matrix
     584RETURN:   list Inv with
     585          - Inv[1] = matrix I and
     586          - Inv[2] = poly p
     587          such that I*A = unitmat(n)*p;
     588NOTE:     p=1 if 1/det(A) is computable and  p=det(A) if not;
     589          the computation uses busadj.
     590EXAMPLE:  example inverse_B; shows an example"
    518591{
    519592  int i;
     
    545618{ "EXAMPLE"; echo=2;
    546619  ring r=0,(x,y),lp;
    547   matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
     620  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
    548621  print(A);
    549622  list Inv=inverse_B(A);
     
    580653"USAGE:     inverse_L(A);  A = square matrix
    581654RETURN:    list Inv representing a left inverse of A, i.e
    582            Inv[1] = matrix I and Inv[2]=poly p such that
    583            I*A = unitmat(n)*p;
    584 NOTE:      p=1 if 1/det(A) is computable and  p=det(A) if not
     655           - Inv[1] = matrix I and
     656           - Inv[2] = poly p
     657           such tha I*A = unitmat(n)*p;
     658NOTE:      p=1 if 1/det(A) is computable and p=det(A) if not;
     659           the computation uses lift
    585660EXAMPLE:   example inverse_L; shows an example"
    586661{
     
    621696{ "EXAMPLE"; echo=2;
    622697  ring r=0,(x,y),lp;
    623   matrix A[3][3]=x,y,1,1,x2,y,x+y,6,7;
     698  matrix A[3][3]=x,y,1,1,x2,y,x,6,0;
    624699  print(A);
    625700  list Inv=inverse_L(A);
     
    631706//////////////////////////////////////////////////////////////////////////////
    632707proc gaussred(matrix A)
    633 "USAGE:     gaussred(A);   A any constant matrix
    634 RETURN:    list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
    635            gives a row reduced matrix S, a permutation matrix P and a
    636            normalized lower triangular matrix U, with P*A=U*S           
    637 EXAMPLE:   example gaussred; shows an example"
     708"USAGE:   gaussred(A);   A any constant matrix
     709RETURN:  list Z:  Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A)
     710         gives a row reduced matrix S, a permutation matrix P and a
     711         normalized lower triangular matrix U, with P*A=U*S   
     712NOTE:    This procedure is designed for teaching purposes mainly.
     713         The straight forward implementation in the interpreted library
     714         is not very efficient (no standard basis computation).         
     715EXAMPLE: example gaussred; shows an example"
    638716{
    639717  int i,j,l,k,jp,rang;
     
    812890proc gauss_nf(matrix A)
    813891"USAGE:    gauss_nf(A); A any constant matrix
    814 RETURN:   matrix; gauss normal form of A
     892RETURN:   matrix; gauss normal form of A (uses gaussred)
    815893EXAMPLE:  example gauss_nf; shows an example"
    816894{
     
    838916  list Z;
    839917  if(!const_mat(A)){
    840     "//input is not a constant matrix";
     918    "// input is not a constant matrix";
    841919    return(-1);
    842920  }
     
    857935           gives a permutation matrix P,
    858936           a normalized lower triangular matrix U ,
    859            a diagonal matrix D, and a normalized upper triangular matrix O
     937           a diagonal matrix D, and
     938           a normalized upper triangular matrix O
    860939           with P*A=U*D*O
    861 NOTE:      Z[1]=-1 means that A is not regular
     940NOTE:      Z[1]=-1 means that A is not regular (proc uses gaussred)
    862941EXAMPLE:   example U_D_O; shows an example"
    863942{
     
    920999proc pos_def(matrix A)
    9211000"USAGE:     pos_def(A); A = constant, symmetric square matrix       
    922 RETURN:    int; 1 if A is positive definit , 0 if not, -1 if unknown
     1001RETURN:    int:
     1002           1  if A is positive definit ,
     1003           0  if not,
     1004           -1 if unknown
    9231005EXAMPLE:   example pos_def; shows an example"
    9241006{
     
    9811063proc linsolve(matrix A, matrix b)
    9821064"USAGE:     linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix
    983 RETURN:    a solution of inhomogeneous linear system A*X = b
     1065RETURN:    a 1xm matrix X, solution of inhomogeneous linear system A*X = b
     1066           return the 0-matrix if system is not solvable
     1067NOTE:      uses gaussred
    9841068EXAMPLE:   example linsolve; shows an example"
    9851069{
     
    10891173//////////////////////////////////////////////////////////////////////////////
    10901174
     1175///////////////////////////////////////////////////////////////////////////////
     1176//    PROCEDURES for Jordan Normal form
     1177//    AUTHOR:  Mathias Schulze, email: mschulze@mathematik.uni-kl.de
     1178///////////////////////////////////////////////////////////////////////////////
     1179
     1180static proc countblocks(matrix M)
     1181{
     1182  int b,r,r0;
     1183
     1184  int i=1;
     1185  while(i<=nrows(M))
     1186  {
     1187    b++;
     1188    r=nrows(M[i]);
     1189    r0=r;
     1190
     1191    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
     1192    while(i<r0&&i<nrows(M))
     1193    {
     1194      i++;
     1195      if(i<=nrows(M))
     1196      {
     1197        r=nrows(M[i]);
     1198        if(r>r0)
     1199        {
     1200          r0=r;
     1201        }
     1202      }
     1203    }
     1204    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
     1205
     1206    i++;
     1207  }
     1208
     1209  return(b);
     1210}
     1211///////////////////////////////////////////////////////////////////////////////
     1212
     1213static proc getblock(matrix M,intvec v)
     1214{
     1215  matrix M0[size(v)][size(v)]=M[v,v];
     1216  return(M0);
     1217}
     1218///////////////////////////////////////////////////////////////////////////////
     1219
     1220proc jordan(matrix M,list #)
     1221"USAGE:   jordan(M[,opt]); M constant square matrix, opt integer
     1222ASSUME:  The eigenvalues of M are in the coefficient field.
     1223RETURN:  The procedure returns a list jd with 3 entries:
     1224@format       
     1225         - jd[1] ideal, eigenvalues of M,
     1226         - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with
     1227           eigenvalue jd[1][i], and
     1228         - jd[3] a matrix, jd[3]^(-1)*M*jd[3] in Jordan normal form.
     1229         Depending on opt, only certain entries of jd are computed.
     1230           If opt=-1, jd[1] is computed,
     1231           if opt= 0, jd[1] and jd[2] are computed,
     1232           if opt= 1, jd[1], jd[2], and jd[3] are computed, and,
     1233           if opt= 2, jd[1] and jd[3] are computed.
     1234         By default, opt=0.
     1235@end format
     1236NOTE:    A non constant polynomial matrix M is replaced by its constant part.
     1237DISPLAY: The procedure displays comments if printlevel>=1.
     1238EXAMPLE: example jordan; shows an example.
     1239"
     1240{
     1241  int n=nrows(M);
     1242  if(n==0)
     1243  {
     1244    print("//empty matrix");
     1245    return(list());
     1246  }
     1247  if(n!=ncols(M))
     1248  {
     1249    print("//no square matrix");
     1250    return(list());
     1251  }
     1252
     1253  M=jet(M,0);
     1254
     1255  dbprint(printlevel-voice+2,"//counting blocks of matrix...");
     1256  int i=countblocks(M);
     1257  dbprint(printlevel-voice+2,"//...blocks of matrix counted");
     1258  if(i==1)
     1259  {
     1260    dbprint(printlevel-voice+2,"//matrix has 1 block");
     1261  }
     1262  else
     1263  {
     1264    dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks");
     1265  }
     1266
     1267  dbprint(printlevel-voice+2,"//counting blocks of transposed matrix...");
     1268  int j=countblocks(transpose(M));
     1269  dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted");
     1270  if(j==1)
     1271  {
     1272    dbprint(printlevel-voice+2,"//transposed matrix has 1 block");
     1273  }
     1274  else
     1275  {
     1276    dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks");
     1277  }
     1278
     1279  if(i<j)
     1280  {
     1281    dbprint(printlevel-voice+2,"//transposing matrix...");
     1282    M=transpose(M);
     1283    dbprint(printlevel-voice+2,"//...matrix transposed");
     1284  }
     1285
     1286  list fd;
     1287  matrix M0;
     1288  poly cp;
     1289  ideal eM,eM0;
     1290  intvec mM,mM0;
     1291  intvec u;
     1292  int b,r,r0;
     1293
     1294  i=1;
     1295  while(i<=nrows(M))
     1296  {
     1297    b++;
     1298    u=i;
     1299    r=nrows(M[i]);
     1300    r0=r;
     1301
     1302    dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
     1303    while(i<r0&&i<nrows(M))
     1304    {
     1305      i++;
     1306      if(i<=nrows(M))
     1307      {
     1308        u=u,i;
     1309        r=nrows(M[i]);
     1310        if(r>r0)
     1311        {
     1312          r0=r;
     1313        }
     1314      }
     1315    }
     1316    dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
     1317
     1318    if(size(u)==1)
     1319    {
     1320      dbprint(printlevel-voice+2,"//1x1-block:");
     1321      dbprint(printlevel-voice+2,M[u[1]][u[1]]);
     1322
     1323      if(mM[1]==0)
     1324      {
     1325        eM=M[u[1]][u[1]];
     1326        mM=1;
     1327      }
     1328      else
     1329      {
     1330        eM=eM,ideal(M[u[1]][u[1]]);
     1331        mM=mM,1;
     1332      }
     1333    }
     1334    else
     1335    {
     1336      dbprint(printlevel-voice+2,
     1337        "//"+string(size(u))+"x"+string(size(u))+"-block:");
     1338      M0=getblock(M,u);
     1339      dbprint(printlevel-voice+2,M0);
     1340
     1341      dbprint(printlevel-voice+2,"//characteristic polynomial:");
     1342      cp=det(module(M0-var(1)*freemodule(size(u))));
     1343      dbprint(printlevel-voice+2,cp);
     1344
     1345      dbprint(printlevel-voice+2,"//factorizing characteristic polynomial...");
     1346      fd=factorize(cp,2);
     1347      dbprint(printlevel-voice+2,"//...characteristic polynomial factorized");
     1348
     1349      dbprint(printlevel-voice+2,"//computing eigenvalues...");
     1350      eM0,mM0=fd[1..2];
     1351      if(1<var(1))
     1352      {
     1353        for(j=ncols(eM0);j>=1;j--)
     1354        {
     1355          if(deg(eM0[j])>1)
     1356          {
     1357            print("//eigenvalues not in the coefficient field");
     1358            return(list());
     1359          }
     1360          if(eM0[j][1]==0)
     1361          {
     1362            eM0[j]=0;
     1363          }
     1364          else
     1365          {
     1366            eM0[j]=-eM0[j][2]/(eM0[j][1]/var(1));
     1367          }
     1368        }
     1369      }
     1370      else
     1371      {
     1372        for(j=ncols(eM0);j>=1;j--)
     1373        {
     1374          if(deg(eM0[j])>1)
     1375          {
     1376            print("//eigenvalues not in the coefficient field");
     1377            return(list());
     1378          }
     1379          if(eM0[j][2]==0)
     1380          {
     1381            eM0[j]=0;
     1382          }
     1383          else
     1384          {
     1385            eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1));
     1386          }
     1387        }
     1388      }
     1389      dbprint(printlevel-voice+2,"//...eigenvalues computed");
     1390
     1391      if(mM[1]==0)
     1392      {
     1393        eM=eM0;
     1394        mM=mM0;
     1395      }
     1396      else
     1397      {
     1398        eM=eM,eM0;
     1399        mM=mM,mM0;
     1400      }
     1401    }
     1402
     1403    i++;
     1404  }
     1405
     1406  dbprint(printlevel-voice+2,"//sorting eigenvalues...");
     1407  poly e;
     1408  int m;
     1409  for(i=ncols(eM);i>=2;i--)
     1410  {
     1411    for(j=i-1;j>=1;j--)
     1412    {
     1413     if(eM[i]<eM[j])
     1414      {
     1415        e=eM[i];
     1416        eM[i]=eM[j];
     1417        eM[j]=e;
     1418        m=mM[i];
     1419        mM[i]=mM[j];
     1420        mM[j]=m;
     1421      }
     1422    }
     1423  }
     1424  dbprint(printlevel-voice+2,"//...eigenvalues sorted");
     1425
     1426  dbprint(printlevel-voice+2,"//removing multiple eigenvalues...");
     1427  i=1;
     1428  j=2;
     1429  while(j<=ncols(eM))
     1430  {
     1431    if(eM[i]==eM[j])
     1432    {
     1433      mM[i]=mM[i]+mM[j];
     1434    }
     1435    else
     1436    {
     1437      i++;
     1438      eM[i]=eM[j];
     1439      mM[i]=mM[j];
     1440    }
     1441    j++;
     1442  }
     1443  eM=eM[1..i];
     1444  mM=mM[1..i];
     1445  dbprint(printlevel-voice+2,"//...multiple eigenvalues removed");
     1446
     1447  dbprint(printlevel-voice+2,"//eigenvalues:");
     1448  dbprint(printlevel-voice+2,eM);
     1449  dbprint(printlevel-voice+2,"//multiplicities:");
     1450  dbprint(printlevel-voice+2,mM);
     1451
     1452  int opt=0;
     1453  if(size(#)>0)
     1454  {
     1455    if(typeof(#[1])=="int")
     1456    {
     1457      opt=#[1];
     1458    }
     1459  }
     1460  if(opt<0)
     1461  {
     1462    return(list(eM));
     1463  }
     1464  int k,l;
     1465  matrix I=freemodule(n);
     1466  matrix Mi,Ni;
     1467  module sNi;
     1468  list K;
     1469  if(opt>=1)
     1470  {
     1471    module V,K1,K2;
     1472    matrix v[n][1];
     1473  }
     1474  if(opt<=1)
     1475  {
     1476    list bM;
     1477    intvec bMi;
     1478  }
     1479
     1480  for(i=ncols(eM);i>=1;i--)
     1481  {
     1482    Mi=M-eM[i]*I;
     1483
     1484    dbprint(printlevel-voice+2,
     1485      "//computing kernels of powers of matrix minus eigenvalue "
     1486      +string(eM[i]));
     1487    K=list(module());
     1488    for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
     1489    {
     1490      sNi=syz(Ni);
     1491      K=K+list(sNi);
     1492    }
     1493    dbprint(printlevel-voice+2,"//...kernels computed");
     1494
     1495    if(opt<=1)
     1496    {
     1497      dbprint(printlevel-voice+2,
     1498        "//computing Jordan block sizes for eigenvalue "
     1499        +string(eM[i])+"...");
     1500      bMi=0;
     1501      bMi[size(K[2])]=0;
     1502      for(j=size(K);j>=2;j--)
     1503      {
     1504        for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--)
     1505        {
     1506          bMi[k]=bMi[k]+1;
     1507        }
     1508      }
     1509      bM=list(bMi)+bM;
     1510      dbprint(printlevel-voice+2,"//...Jordan block sizes computed");
     1511    }
     1512
     1513    if(opt>=1)
     1514    {
     1515      dbprint(printlevel-voice+2,
     1516        "//computing Jordan basis vectors for eigenvalue "
     1517        +string(eM[i])+"...");
     1518      if(size(K)>1)
     1519      {
     1520        for(j,K1=2,0;j<=size(K)-1;j++)
     1521        {
     1522          K2=K[j];
     1523          K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1]))));
     1524          K1=K2;
     1525        }
     1526        K[j]=interred(reduce(K[j],std(K1)));
     1527      }
     1528      for(j=size(K);j>=2;j--)
     1529      {
     1530        for(k=size(K[j]);k>=1;k--)
     1531        {
     1532          v=K[j][k];
     1533          for(l=j;l>=1;l--)
     1534          {
     1535            V=module(v)+V;
     1536            v=Mi*v;
     1537          }
     1538        }
     1539      }
     1540      dbprint(printlevel-voice+2,"//...Jordan basis vectors computed");
     1541    }
     1542  }
     1543
     1544  list jd=eM;
     1545  if(opt<=1)
     1546  {
     1547    jd[2]=bM;
     1548  }
     1549  if(opt>=1)
     1550  {
     1551    jd[3]=V;
     1552  }
     1553  return(jd);
     1554}
     1555example
     1556{ "EXAMPLE:"; echo=2;
     1557  ring R=0,x,dp;
     1558  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
     1559  print(M);
     1560  jordan(M);
     1561}
     1562///////////////////////////////////////////////////////////////////////////////
     1563
     1564proc jordanmatrix(list jd)
     1565"USAGE:   jordanmatrix(jd); jd list of ideal and list of intvecs
     1566RETURN:  The procedure returns the Jordan matrix J, satisfying:
     1567@*       - eigenvalues of J are given by jd[1]
     1568@*       - j-th Jordan block of J with eigenvalue jd[1][i] has size jd[2][i][j]
     1569DISPLAY: The procedure displays comments if printlevel>=1.
     1570EXAMPLE: example jordanmatrix; shows an example.
     1571"
     1572{
     1573  if(size(jd)<2)
     1574  {
     1575    print("//not enough entries in argument list");
     1576    matrix J[1][0];
     1577    return(J);
     1578  }
     1579  def eJ,bJ=jd[1..2];
     1580  if(typeof(eJ)!="ideal")
     1581  {
     1582    print("//first entry in argument list not an ideal");
     1583    matrix J[1][0];
     1584    return(J);
     1585  }
     1586  if(typeof(bJ)!="list")
     1587  {
     1588    print("//second entry in argument list not a list");
     1589    matrix J[1][0];
     1590    return(J);
     1591  }
     1592  if(size(eJ)<size(bJ))
     1593  {
     1594    int s=size(eJ);
     1595  }
     1596  else
     1597  {
     1598    int s=size(bJ);
     1599  }
     1600
     1601  int i,j,k,n;
     1602  for(i=s;i>=1;i--)
     1603  {
     1604    if(typeof(bJ[i])!="intvec")
     1605    {
     1606      print("//second entry in argument list not a list of intvecs");
     1607      matrix J[1][0];
     1608      return(J);
     1609    }
     1610    else
     1611    {
     1612      for(j=size(bJ[i]);j>=1;j--)
     1613      {
     1614        k=bJ[i][j];
     1615        if(k>0)
     1616        {
     1617          n=n+k;
     1618        }
     1619      }
     1620    }
     1621  }
     1622
     1623  int l;
     1624  matrix J[n][n];
     1625  for(i,l=1,1;i<=s;i++)
     1626  {
     1627    for(j=1;j<=size(bJ[i]);j++)
     1628    {
     1629      k=bJ[i][j];
     1630      if(k>0)
     1631      {
     1632        while(k>=2)
     1633        {
     1634          J[l,l]=eJ[i];
     1635          J[l,l+1]=1;
     1636          k,l=k-1,l+1;
     1637        }
     1638        J[l,l]=eJ[i];
     1639        l++;
     1640      }
     1641    }
     1642  }
     1643
     1644  return(J);
     1645}
     1646example
     1647{ "EXAMPLE:"; echo=2;
     1648  ring R=0,x,dp;
     1649  list l;
     1650  l[1]=ideal(2,3);
     1651  l[2]=list(intvec(1),intvec(2));
     1652  print(jordanmatrix(l));
     1653}
     1654///////////////////////////////////////////////////////////////////////////////
     1655
     1656proc jordanform(matrix M)
     1657"USAGE:   jordanform(M); M constant square matrix
     1658ASSUME:  The eigenvalues of M are in the coefficient field.
     1659RETURN:  matrix, the Jordan normal form of M.
     1660NOTE:    A non constant polynomial matrix M is replaced by its constant part.
     1661DISPLAY: The procedure displays more comments for higher printlevel.
     1662EXAMPLE: example jordanform; shows an example.
     1663"
     1664{
     1665  return(jordanmatrix(jordan(M)));
     1666}
     1667example
     1668{ "EXAMPLE:"; echo=2;
     1669  ring R=0,x,dp;
     1670  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
     1671  print(M);
     1672  print(jordanform(M));
     1673}
     1674///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.