Changeset 65b27c in git


Ignore:
Timestamp:
Mar 5, 2001, 7:28:49 PM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
Children:
457d8d6b8765ca2fc446fb19f1402006dffcebaa
Parents:
41f495bc8a4fba73437044f667b69ab85490137f
Message:
*mschulze: jordan and gaussman use eigenval now


git-svn-id: file:///usr/local/Singular/svn/trunk@5311 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
Singular
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gaussman.lib

    r41f495 r65b27c  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: gaussman.lib,v 1.33 2001-02-08 14:20:29 mschulze Exp $";
     2version="$Id: gaussman.lib,v 1.34 2001-03-05 18:28:46 mschulze Exp $";
    33category="Singularities";
    44
    … …  
    145145
    146146    dbprint(printlevel-voice+2,"//gaussman::monodromy: compute C");
    147     C=coeffs(reduce(w,sJ,U),m);
     147    C=coeffs(reduce(w,U,sJ),m);
    148148    A0=A0+C*var(1)^k;
    149149
    … …  
    177177        "//gaussman::monodromy: compute A on saturation");
    178178      l=division(H*var(1),A0*H+var(1)^2*diff(matrix(H),var(1)));
    179       A=jet(l[1],N-1,l[2]);
     179      A=jet(l[1],l[2],N-1);
    180180      if(mide==0)
    181181      {
    … …  
    183183          "//gaussman::monodromy: compute eigenvalues e and"+
    184184          "multiplicities b of A");
    185         l=jordan(A);
     185        l=system("eigenval",jet(A,0));
    186186        ideal e=l[1];
    187         intvec b;
    188         for(i=ncols(e);i>=1;i--)
    189         {
    190           b[i]=sum(l[2][i]);
    191         }
     187        intvec b=l[2];
    192188        dbprint(printlevel-voice+2,"//gaussman::monodromy: e="+string(e));
    193189        dbprint(printlevel-voice+2,"//gaussman::monodromy: b="+string(b));
    … …  
    397393
    398394    dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute C");
    399     C=coeffs(reduce(w,sJ,U),m);
     395    C=coeffs(reduce(w,U,sJ),m);
    400396    A=A+C*var(1)^k;
    401397
    … …  
    455451    "//gaussman::vfiltration: transform H0 to saturation");
    456452  l=division(H,H0);
    457   H0=jet(l[1],N-1,l[2]);
     453  H0=jet(l[1],l[2],N-1);
    458454
    459455  dbprint(printlevel-voice+2,
    … …  
    490486    "//gaussman::vfiltration: compute A on saturation");
    491487  l=division(H*var(1),A*H+var(1)^2*diff(matrix(H),var(1)));
    492   A=jet(l[1],N-1,l[2]);
     488  A=jet(l[1],l[2],N-1);
    493489
    494490  dbprint(printlevel-voice+2,"//gaussman::vfiltration: compute matrix M of A");
    … …  
    520516  dbprint(printlevel-voice+2,
    521517    "//gaussman::vfiltration: compute eigenvalues eA of A");
    522   ideal eA=jordan(A,-1)[1];
     518  ideal eA=system("eigenval",jet(A,0))[1];
    523519  dbprint(printlevel-voice+2,"//gaussman::vfiltration: eA="+string(eA));
    524520
    … …  
    731727  for(i=ncols(m);i>=1;i--)
    732728  {
    733     M[i]=lift(V,coeffs(reduce(m[i]*m,sJ,U),m)*V);
     729    M[i]=lift(V,coeffs(reduce(m[i]*m,U,sJ),m)*V);
    734730  }
    735731

  • Singular/LIB/linalg.lib

    r41f495 r65b27c  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.10 2001-01-11 17:12:59 Singular Exp $";
     3version="$Id: linalg.lib,v 1.11 2001-03-05 18:28:45 mschulze Exp $";
    44category="Linear Algebra";
    55info="
    … …  
    2525 U_D_O(A);          P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
    2626 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
     27 jordan(M);         Jordan data of constant square matrix M
     28 jordanbasis(M);    Jordan basis of constant square matrix M
     29 jordanmatrix(e,b); Jordan matrix with eigenvalues e and Jordan block sizes b
     30 jordanform(M);     Jordan matrix of constant square matrix M
    3031";
    3132
    … …  
    12331234///////////////////////////////////////////////////////////////////////////////
    12341235
    1235 static proc countblocks(matrix M)
    1236 {
    1237   int b,r,r0;
    1238 
    1239   int i=1;
    1240   while(i<=nrows(M))
     1236proc jordan(matrix M,list #)
     1237"USAGE:   jordan(M); matrix M
     1238ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
     1239RETURN:
     1240@format
     1241<list> l:
     1242  <ideal> l[1]: eigenvalues of M
     1243  <list> l[2]:
     1244    <intvec> l[2][i]:
     1245       <int> l[2][i][j]: size of Jordan block j of M with eigenvalue l[1][i]
     1246@end format
     1247EXAMPLE: example jordan; shows an example
     1248"
     1249{
     1250  if(nrows(M)==0)
    12411251  {
    1242     b++;
    1243     r=nrows(M[i]);
    1244     r0=r;
    1245 
    1246     dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    1247     while(i<r0&&i<nrows(M))
     1252    ERROR("non empty expected");
     1253  }
     1254  if(ncols(M)!=nrows(M))
     1255  {
     1256    ERROR("square matrix expected");
     1257  }
     1258
     1259  M=jet(M,0);
     1260
     1261  list l=system("eigenval",M);
     1262  def e,m=l[1..2];
     1263
     1264  int i;
     1265  for(i=1;i<=ncols(e);i++)
     1266  {
     1267    if(deg(e[i]>0))
    12481268    {
    1249       i++;
    1250       if(i<=nrows(M))
     1269      ERROR("eigenvalues in coefficient field expected");
     1270      return(list());
     1271    }
     1272  }
     1273
     1274  int j,k;
     1275  matrix N,NN;
     1276  module K0;
     1277  list K;
     1278  intvec b0;
     1279  list b;
     1280
     1281  for(i=ncols(e);i>0;i--)
     1282  {
     1283    N=M-e[i]*freemodule(ncols(M));
     1284
     1285    K0=0;
     1286    NN=N;
     1287    K=module();
     1288    while(size(K0)<m[i])
     1289    {
     1290      K0=syz(NN);
     1291      K=K+list(K0);
     1292      NN=NN*N;
     1293    }
     1294
     1295    b0=0;
     1296    b0[size(K[2])]=0;
     1297    for(j=size(K);j>1;j--)
     1298    {
     1299      for(k=size(b0);k>size(b0)+size(K[j-1])-size(K[j]);k--)
    12511300      {
    1252         r=nrows(M[i]);
    1253         if(r>r0)
    1254         {
    1255           r0=r;
    1256         }
    1257       }
    1258     }
    1259     dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
    1260 
    1261     i++;
    1262   }
    1263 
    1264   return(b);
    1265 }
    1266 ///////////////////////////////////////////////////////////////////////////////
    1267 
    1268 static proc getblock(matrix M,intvec v)
    1269 {
    1270   matrix M0[size(v)][size(v)]=M[v,v];
    1271   return(M0);
    1272 }
    1273 ///////////////////////////////////////////////////////////////////////////////
    1274 
    1275 proc jordan(matrix M,list #)
    1276 "USAGE:   jordan(M[,opt]); M constant square matrix, opt integer
    1277 ASSUME:  The eigenvalues of M are in the coefficient field.
    1278 RETURN:  The procedure returns a list jd with 3 entries:
    1279 @format
    1280          - jd[1] ideal, eigenvalues of M,
    1281          - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with
    1282            eigenvalue jd[1][i], and
    1283          - jd[3] a matrix, jd[3]^(-1)*M*jd[3] in Jordan normal form.
    1284          Depending on opt, only certain entries of jd are computed.
    1285            If opt=-1, jd[1] is computed,
    1286            if opt= 0, jd[1] and jd[2] are computed,
    1287            if opt= 1, jd[1], jd[2], and jd[3] are computed, and,
    1288            if opt= 2, jd[1] and jd[3] are computed.
    1289          By default, opt=0.
    1290 @end format
    1291 NOTE:    A non constant polynomial matrix M is replaced by its constant part.
    1292 DISPLAY: The procedure displays comments if printlevel>=1.
    1293 EXAMPLE: example jordan; shows an example.
    1294 "
    1295 {
    1296   int n=nrows(M);
    1297   if(n==0)
    1298   {
    1299     print("//empty matrix");
    1300     return(list());
    1301   }
    1302   if(n!=ncols(M))
    1303   {
    1304     print("//no square matrix");
    1305     return(list());
    1306   }
    1307 
    1308   M=jet(M,0);
    1309 
    1310   dbprint(printlevel-voice+2,"//counting blocks of matrix...");
    1311   int i=countblocks(M);
    1312   dbprint(printlevel-voice+2,"//...blocks of matrix counted");
    1313   if(i==1)
    1314   {
    1315     dbprint(printlevel-voice+2,"//matrix has 1 block");
    1316   }
    1317   else
    1318   {
    1319     dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks");
    1320   }
    1321 
    1322   dbprint(printlevel-voice+2,"//counting blocks of transposed matrix...");
    1323   int j=countblocks(transpose(M));
    1324   dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted");
    1325   if(j==1)
    1326   {
    1327     dbprint(printlevel-voice+2,"//transposed matrix has 1 block");
    1328   }
    1329   else
    1330   {
    1331     dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks");
    1332   }
    1333 
    1334   if(i<j)
    1335   {
    1336     dbprint(printlevel-voice+2,"//transposing matrix...");
    1337     M=transpose(M);
    1338     dbprint(printlevel-voice+2,"//...matrix transposed");
    1339   }
    1340 
    1341   list fd;
    1342   matrix M0;
    1343   poly cp;
    1344   ideal eM,eM0;
    1345   intvec mM,mM0;
    1346   intvec u;
    1347   int b,r,r0;
    1348 
    1349   i=1;
    1350   while(i<=nrows(M))
    1351   {
    1352     b++;
    1353     u=i;
    1354     r=nrows(M[i]);
    1355     r0=r;
    1356 
    1357     dbprint(printlevel-voice+2,"//searching for block "+string(b)+"...");
    1358     while(i<r0&&i<nrows(M))
    1359     {
    1360       i++;
    1361       if(i<=nrows(M))
    1362       {
    1363         u=u,i;
    1364         r=nrows(M[i]);
    1365         if(r>r0)
    1366         {
    1367           r0=r;
    1368         }
    1369       }
    1370     }
    1371     dbprint(printlevel-voice+2,"//...block "+string(b)+" found");
    1372 
    1373     if(size(u)==1)
    1374     {
    1375       dbprint(printlevel-voice+2,"//1x1-block:");
    1376       dbprint(printlevel-voice+2,M[u[1]][u[1]]);
    1377 
    1378       if(mM[1]==0)
    1379       {
    1380         eM=M[u[1]][u[1]];
    1381         mM=1;
    1382       }
    1383       else
    1384       {
    1385         eM=eM,ideal(M[u[1]][u[1]]);
    1386         mM=mM,1;
    1387       }
    1388     }
    1389     else
    1390     {
    1391       dbprint(printlevel-voice+2,
    1392         "//"+string(size(u))+"x"+string(size(u))+"-block:");
    1393       M0=getblock(M,u);
    1394       dbprint(printlevel-voice+2,M0);
    1395 
    1396       dbprint(printlevel-voice+2,"//characteristic polynomial:");
    1397       cp=det(module(M0-var(1)*freemodule(size(u))));
    1398       dbprint(printlevel-voice+2,cp);
    1399 
    1400       dbprint(printlevel-voice+2,"//factorizing characteristic polynomial...");
    1401       fd=factorize(cp,2);
    1402       dbprint(printlevel-voice+2,"//...characteristic polynomial factorized");
    1403 
    1404       dbprint(printlevel-voice+2,"//computing eigenvalues...");
    1405       eM0,mM0=fd[1..2];
    1406       if(1<var(1))
    1407       {
    1408         for(j=ncols(eM0);j>=1;j--)
    1409         {
    1410           if(deg(eM0[j])>1)
    1411           {
    1412             print("//eigenvalues not in the coefficient field");
    1413             return(list());
    1414           }
    1415           if(eM0[j][1]==0)
    1416           {
    1417             eM0[j]=0;
    1418           }
    1419           else
    1420           {
    1421             eM0[j]=-eM0[j][2]/(eM0[j][1]/var(1));
    1422           }
    1423         }
    1424       }
    1425       else
    1426       {
    1427         for(j=ncols(eM0);j>=1;j--)
    1428         {
    1429           if(deg(eM0[j])>1)
    1430           {
    1431             print("//eigenvalues not in the coefficient field");
    1432             return(list());
    1433           }
    1434           if(eM0[j][2]==0)
    1435           {
    1436             eM0[j]=0;
    1437           }
    1438           else
    1439           {
    1440             eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1));
    1441           }
    1442         }
    1443       }
    1444       dbprint(printlevel-voice+2,"//...eigenvalues computed");
    1445 
    1446       if(mM[1]==0)
    1447       {
    1448         eM=eM0;
    1449         mM=mM0;
    1450       }
    1451       else
    1452       {
    1453         eM=eM,eM0;
    1454         mM=mM,mM0;
    1455       }
    1456     }
    1457 
    1458     i++;
    1459   }
    1460 
    1461   dbprint(printlevel-voice+2,"//sorting eigenvalues...");
    1462   poly e;
    1463   int m;
    1464   for(i=ncols(eM);i>=2;i--)
    1465   {
    1466     for(j=i-1;j>=1;j--)
    1467     {
    1468      if(eM[i]<eM[j])
    1469       {
    1470         e=eM[i];
    1471         eM[i]=eM[j];
    1472         eM[j]=e;
    1473         m=mM[i];
    1474         mM[i]=mM[j];
    1475         mM[j]=m;
    1476       }
    1477     }
    1478   }
    1479   dbprint(printlevel-voice+2,"//...eigenvalues sorted");
    1480 
    1481   dbprint(printlevel-voice+2,"//removing multiple eigenvalues...");
    1482   i=1;
    1483   j=2;
    1484   while(j<=ncols(eM))
    1485   {
    1486     if(eM[i]==eM[j])
    1487     {
    1488       mM[i]=mM[i]+mM[j];
    1489     }
    1490     else
    1491     {
    1492       i++;
    1493       eM[i]=eM[j];
    1494       mM[i]=mM[j];
    1495     }
    1496     j++;
    1497   }
    1498   eM=eM[1..i];
    1499   mM=mM[1..i];
    1500   dbprint(printlevel-voice+2,"//...multiple eigenvalues removed");
    1501 
    1502   dbprint(printlevel-voice+2,"//eigenvalues:");
    1503   dbprint(printlevel-voice+2,eM);
    1504   dbprint(printlevel-voice+2,"//multiplicities:");
    1505   dbprint(printlevel-voice+2,mM);
    1506 
    1507   int opt=0;
    1508   if(size(#)>0)
    1509   {
    1510     if(typeof(#[1])=="int")
    1511     {
    1512       opt=#[1];
    1513     }
    1514   }
    1515   if(opt<0)
    1516   {
    1517     return(list(eM));
    1518   }
    1519   int k,l;
    1520   matrix I=freemodule(n);
    1521   matrix Mi,Ni;
    1522   module sNi;
    1523   list K;
    1524   if(opt>=1)
    1525   {
    1526     module V,K1,K2;
    1527     matrix v[n][1];
    1528   }
    1529   if(opt<=1)
    1530   {
    1531     list bM;
    1532     intvec bMi;
    1533   }
    1534 
    1535   for(i=ncols(eM);i>=1;i--)
    1536   {
    1537     Mi=M-eM[i]*I;
    1538 
    1539     dbprint(printlevel-voice+2,
    1540       "//computing kernels of powers of matrix minus eigenvalue "
    1541       +string(eM[i]));
    1542     K=list(module());
    1543     for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi)
    1544     {
    1545       sNi=syz(Ni);
    1546       K=K+list(sNi);
    1547     }
    1548     dbprint(printlevel-voice+2,"//...kernels computed");
    1549 
    1550     if(opt<=1)
    1551     {
    1552       dbprint(printlevel-voice+2,
    1553         "//computing Jordan block sizes for eigenvalue "
    1554         +string(eM[i])+"...");
    1555       bMi=0;
    1556       bMi[size(K[2])]=0;
    1557       for(j=size(K);j>=2;j--)
    1558       {
    1559         for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--)
    1560         {
    1561           bMi[k]=bMi[k]+1;
    1562         }
    1563       }
    1564       bM=list(bMi)+bM;
    1565       dbprint(printlevel-voice+2,"//...Jordan block sizes computed");
    1566     }
    1567 
    1568     if(opt>=1)
    1569     {
    1570       dbprint(printlevel-voice+2,
    1571         "//computing Jordan basis vectors for eigenvalue "
    1572         +string(eM[i])+"...");
    1573       if(size(K)>1)
    1574       {
    1575         for(j,K1=2,0;j<=size(K)-1;j++)
    1576         {
    1577           K2=K[j];
    1578           K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1]))));
    1579           K1=K2;
    1580         }
    1581         K[j]=interred(reduce(K[j],std(K1)));
    1582       }
    1583       for(j=size(K);j>=2;j--)
    1584       {
    1585         for(k=size(K[j]);k>=1;k--)
    1586         {
    1587           v=K[j][k];
    1588           for(l=j;l>=1;l--)
    1589           {
    1590             V=module(v)+V;
    1591             v=Mi*v;
    1592           }
    1593         }
    1594       }
    1595       dbprint(printlevel-voice+2,"//...Jordan basis vectors computed");
    1596     }
    1597   }
    1598 
    1599   list jd=eM;
    1600   if(opt<=1)
    1601   {
    1602     jd[2]=bM;
    1603   }
    1604   if(opt>=1)
    1605   {
    1606     jd[3]=V;
    1607   }
    1608   return(jd);
     1301        b0[k]=b0[k]+1;
     1302      }
     1303    }
     1304    b=list(b0)+b;
     1305  }
     1306
     1307  l[2]=b;
     1308  return(l);
    16091309}
    16101310example
    … …  
    16171317///////////////////////////////////////////////////////////////////////////////
    16181318
    1619 proc jordanmatrix(list jd)
    1620 "USAGE:   jordanmatrix(jd); jd list of ideal and list of intvecs
    1621 RETURN:  The procedure returns the Jordan matrix J, satisfying:
    1622 @*       - eigenvalues of J are given by jd[1]
    1623 @*       - j-th Jordan block of J with eigenvalue jd[1][i] has size jd[2][i][j]
    1624 DISPLAY: The procedure displays comments if printlevel>=1.
    1625 EXAMPLE: example jordanmatrix; shows an example.
     1319proc jordanbasis(matrix M)
     1320"USAGE:   jordan(M); matrix M
     1321ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
     1322RETURN:  <matrix> U: inverse(U)*M*U Jordan normal form of M
     1323EXAMPLE: example jordanbasis; shows an example
    16261324"
    16271325{
    1628   if(size(jd)<2)
     1326  if(nrows(M)==0)
    16291327  {
    1630     print("//not enough entries in argument list");
    1631     matrix J[1][0];
    1632     return(J);
    1633   }
    1634   def eJ,bJ=jd[1..2];
    1635   if(typeof(eJ)!="ideal")
     1328    ERROR("non empty matrix expected");
     1329  }
     1330  if(ncols(M)!=nrows(M))
    16361331  {
    1637     print("//first entry in argument list not an ideal");
    1638     matrix J[1][0];
    1639     return(J);
    1640   }
    1641   if(typeof(bJ)!="list")
     1332    ERROR("square matrix expected");
     1333  }
     1334
     1335  M=jet(M,0);
     1336
     1337  list l=system("eigenval",M);
     1338  def e,m=l[1..2];
     1339  kill l;
     1340
     1341  int i;
     1342  for(i=1;i<=ncols(e);i++)
    16421343  {
    1643     print("//second entry in argument list not a list");
    1644     matrix J[1][0];
    1645     return(J);
    1646   }
    1647   if(size(eJ)<size(bJ))
     1344    if(deg(e[i]>0))
     1345    {
     1346      ERROR("eigenvalues in coefficient field expected");
     1347      return(freemodule(ncols(M)));
     1348    }
     1349  }
     1350
     1351  int j,k,l;
     1352  matrix N,NN;
     1353  module K0,K1;
     1354  list K;
     1355  matrix u[ncols(M)][1];
     1356  module U;
     1357
     1358  for(i=ncols(e);i>0;i--)
    16481359  {
    1649     int s=size(eJ);
    1650   }
    1651   else
     1360    N=M-e[i]*freemodule(ncols(M));
     1361
     1362    K0=0;
     1363    NN=N;
     1364    K=module();
     1365    while(size(K0)<m[i])
     1366    {
     1367      K0=syz(NN);
     1368      K=K+list(K0);
     1369      NN=NN*N;
     1370    }
     1371
     1372    K1=0;
     1373    for(j=2;j<size(K);j++)
     1374    {
     1375      K0=K[j];
     1376      K[j]=interred(reduce(K[j],std(K1+module(N*K[j+1]))));
     1377      K1=K0;
     1378    }
     1379    K[j]=interred(reduce(K[j],std(K1)));
     1380
     1381    for(j=size(K);j>1;j--)
     1382    {
     1383      for(k=size(K[j]);k>0;k--)
     1384      {
     1385        u=K[j][k];
     1386        for(l=j;l>0;l--)
     1387        {
     1388          U=module(u)+U;
     1389          u=N*u;
     1390        }
     1391      }
     1392    }
     1393  }
     1394
     1395  return(U);
     1396}
     1397example
     1398{ "EXAMPLE:"; echo=2;
     1399  ring R=0,x,dp;
     1400  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
     1401  print(M);
     1402  matrix U=jordanbasis(M);
     1403  print(U);
     1404  print(inverse(U)*M*U);
     1405}
     1406///////////////////////////////////////////////////////////////////////////////
     1407
     1408proc jordanmatrix(ideal e,list b)
     1409"USAGE:   jordanmatrix(e,b); ideal e, list b
     1410RETURN:  <matrix> J: Jordan matrix with eigenvalues e and block sizes b
     1411EXAMPLE: example jordanmatrix; shows an example
     1412"
     1413{
     1414  if(size(e)!=size(b))
    16521415  {
    1653     int s=size(bJ);
     1416    ERROR("arguments of equal size expected");
    16541417  }
    16551418
    16561419  int i,j,k,n;
    1657   for(i=s;i>=1;i--)
     1420  for(i=size(e);i>=1;i--)
    16581421  {
    1659     if(typeof(bJ[i])!="intvec")
     1422    if(typeof(b[i])!="intvec")
    16601423    {
    1661       print("//second entry in argument list not a list of intvecs");
    1662       matrix J[1][0];
    1663       return(J);
     1424      ERROR("second argument of type list of intvec expected");
    16641425    }
    16651426    else
    16661427    {
    1667       for(j=size(bJ[i]);j>=1;j--)
     1428      for(j=size(b[i]);j>=1;j--)
    16681429      {
    1669         k=bJ[i][j];
     1430        k=b[i][j];
    16701431        if(k>0)
    16711432        {
    … …  
    16751436    }
    16761437  }
    1677 
    1678   int l;
    16791438  matrix J[n][n];
    1680   for(i,l=1,1;i<=s;i++)
     1439
     1440  int l=1;
     1441  for(i=1;i<=size(e);i++)
    16811442  {
    1682     for(j=1;j<=size(bJ[i]);j++)
     1443    for(j=1;j<=size(b[i]);j++)
    16831444    {
    1684       k=bJ[i][j];
     1445      k=b[i][j];
    16851446      if(k>0)
    16861447      {
    16871448        while(k>=2)
    16881449        {
    1689           J[l,l]=eJ[i];
     1450          J[l,l]=e[i];
    16901451          J[l,l+1]=1;
    1691           k,l=k-1,l+1;
     1452          k--;
     1453          l++;
    16921454        }
    1693         J[l,l]=eJ[i];
     1455        J[l,l]=e[i];
    16941456        l++;
    16951457      }
    … …  
    17021464{ "EXAMPLE:"; echo=2;
    17031465  ring R=0,x,dp;
    1704   list l;
    1705   l[1]=ideal(2,3);
    1706   l[2]=list(intvec(1),intvec(2));
    1707   print(jordanmatrix(l));
     1466  ideal e=ideal(2,3);
     1467  list b=list(intvec(1),intvec(2));
     1468  print(jordanmatrix(e,b));
    17081469}
    17091470///////////////////////////////////////////////////////////////////////////////
    17101471
    17111472proc jordanform(matrix M)
    1712 "USAGE:   jordanform(M); M constant square matrix
    1713 ASSUME:  The eigenvalues of M are in the coefficient field.
    1714 RETURN:  matrix, the Jordan normal form of M.
    1715 NOTE:    A non constant polynomial matrix M is replaced by its constant part.
    1716 DISPLAY: The procedure displays more comments for higher printlevel.
    1717 EXAMPLE: example jordanform; shows an example.
     1473"USAGE:   jordanform(M); matrix M
     1474ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
     1475RETURN:  <matrix> J: Jordan normal form of M
     1476EXAMPLE: example jordanform; shows an example
    17181477"
    17191478{
    1720   return(jordanmatrix(jordan(M)));
     1479  list l=jordan(M);
     1480  return(jordanmatrix(l[1],l[2]));
    17211481}
    17221482example
    … …  
    17271487  print(jordanform(M));
    17281488}
    1729 
    17301489///////////////////////////////////////////////////////////////////////////////
    17311490

  • Singular/Makefile.in

    r41f495 r65b27c  
    112112    GMPrat.cc multicnt.cc npolygon.cc semic.cc spectrum.cc splist.cc \
    113113    libparse.cc mod_raw.cc sing_win.cc\
    114     pcv.cc units.cc kbuckets.cc sbuckets.cc\
     114    pcv.cc units.cc eigenval.cc kbuckets.cc sbuckets.cc\
    115115    mpr_inout.cc mpr_base.cc mpr_numeric.cc \
    116116    prCopy.cc p_Mult_q.cc \
    … …  
    174174        ndbm.h dbm_sl.h polys-impl.h libparse.h \
    175175        GMPrat.h multicnt.h npolygon.h semic.h spectrum.h splist.h multicnt.h \
    176         pcv.h units.h mod_raw.h kbuckets.h sbuckets.h\
     176        pcv.h units.h eigenval.h mod_raw.h kbuckets.h sbuckets.h\
    177177        mpr_global.h mpr_inout.h mpr_base.h mpr_numeric.h \
    178178        feOpt.h fegetopt.h distrib.h walk.h \

  • Singular/extra.cc

    r41f495 r65b27c  
    22*  Computer Algebra System SINGULAR      *
    33*****************************************/
    4 /* $Id: extra.cc,v 1.163 2001-02-26 15:08:42 levandov Exp $ */
     4/* $Id: extra.cc,v 1.164 2001-03-05 18:28:49 mschulze Exp $ */
    55/*
    66* ABSTRACT: general interface to internals of Singular ("system" command)
    … …  
    106106#endif
    107107#endif /* not HAVE_DYNAMIC_LOADING */
     108
     109// eigenvalues of constant square matrices
     110#ifdef HAVE_EIGENVAL
     111#include "eigenval.h"
     112#endif
    108113
    109114// see clapsing.cc for a description of the `FACTORY_*' options
    … …  
    543548#endif
    544549#endif /* HAVE_DYNAMIC_LOADING */
     550/*==================== eigenval =============================*/
     551    if(strcmp(sys_cmd,"tridiag")==0)
     552    {
     553      return tridiag(res,h);
     554    }
     555    else
     556    if(strcmp(sys_cmd,"eigenval")==0)
     557    {
     558      return eigenval(res,h);
     559    }
     560    else
    545561/*==================== contributors =============================*/
    546562   if(strcmp(sys_cmd,"contributors") == 0)
    … …  
    548564     res->rtyp=STRING_CMD;
    549565     res->data=(void *)omStrDup(
    550        "Olaf Bachmann, Hubert Grassmann, Kai Krueger, Wolfgang Neumann, Thomas Nuessler, Wilfred Pohl, Jens Schmidt, Thomas Siebert, Ruediger Stobbe, Moritz Wenk, Tim Wichmann");
     566       "Olaf Bachmann, Hubert Grassmann, Kai Krueger, Wolfgang Neumann, Thomas Nuessler, Wilfred Pohl, Jens Schmidt, Mathias Schulze, Thomas Siebert, Ruediger Stobbe, Moritz Wenk, Tim Wichmann");
    551567     return FALSE;
    552568   }

  • Singular/mod2.h.in

    r41f495 r65b27c  
    66 *          DO NOT EDIT!
    77 *
    8  *  Version: $Id: mod2.h.in,v 1.99 2001-02-02 11:32:27 mschulze Exp $
     8 *  Version: $Id: mod2.h.in,v 1.100 2001-03-05 18:28:49 mschulze Exp $
    99 *******************************************************************/
    1010#ifndef MOD2_H
    … …  
    171171#define HAVE_UNITS
    172172
     173// eigenvalues of constant square matrices
     174#define HAVE_EIGENVAL
     175
    173176/* Define to use old mechanismen for saving currRing with procedures
    174177 * Does not work with HAVE_NAMESPACES enabled
Note: See TracChangeset for help on using the changeset viewer.