Changeset e9124e in git


Ignore:
Timestamp:
Mar 6, 2002, 5:39:09 PM (21 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
Children:
101cd892be49d477f3f2ab71d77dbfe943baa01b
Parents:
fee17e544a54ec29e804c68a244e859ae9a1bfd3
Message:
*mschulze: check system-with for eigenval, gms


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/gaussman.lib

    rfee17e re9124e  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: gaussman.lib,v 1.74 2002-03-05 14:22:14 mschulze Exp $";
     2version="$Id: gaussman.lib,v 1.75 2002-03-06 16:39:09 mschulze Exp $";
    33category="Singularities";
    44
    … …  
    88AUTHOR:   Mathias Schulze, email: mschulze@mathematik.uni-kl.de
    99
    10 OVERVIEW: A library to compute Hodge-theoretic invariants 
     10OVERVIEW: A library to compute Hodge-theoretic invariants
    1111          of isolated hypersurface singularities
    1212
    … …  
    4242";
    4343
    44 LIB "linalg2.lib";
     44LIB "linalg.lib";
    4545
    4646///////////////////////////////////////////////////////////////////////////////
    … …  
    151151      ERROR("isolated critical point 0 expected");
    152152    }
    153   } 
     153  }
    154154  matrix B=l[2];
    155155  ideal m=kbase(g);
    … …  
    217217RETURN:
    218218@format
    219 list nf; 
     219list nf;
    220220  ideal nf[1];  projection of p to gmsbasis mod s^(K+1)
    221221  ideal nf[2];  p=nf[1]+nf[2]
    … …  
    226226"
    227227{
    228   return(system("gmsnf",p,gmsstd,gmsmatrix,(K+1)*deg(var(1))-2*gmsmaxdeg,K));
     228  if(system("with","gms"))
     229  {
     230    return(system("gmsnf",p,gmsstd,gmsmatrix,(K+1)*deg(var(1))-2*gmsmaxdeg,K));
     231  }
     232
     233  intvec v=1;
     234  v[nvars(basering)]=0;
     235
     236  int k;
     237  ideal r,q;
     238  r[ncols(p)]=0;
     239  q[ncols(p)]=0;
     240
     241  poly s;
     242  int i,j;
     243  for(k=ncols(p);k>=1;k--)
     244  {
     245    while(p[k]!=0&&deg(lead(p[k]),v)<=K)
     246    {
     247      i=1;
     248      s=lead(p[k])/lead(gmsstd[i]);
     249      while(i<ncols(gmsstd)&&s==0)
     250      {
     251        i++;
     252        s=lead(p[k])/lead(gmsstd[i]);
     253      }
     254      if(s!=0)
     255      {
     256        p[k]=p[k]-s*gmsstd[i];
     257        for(j=1;j<=nrows(gmsmatrix);j++)
     258        {
     259          p[k]=p[k]+diff(s*gmsmatrix[j,i],var(j+1))*var(1);
     260        }
     261      }
     262      else
     263      {
     264        r[k]=r[k]+lead(p[k]);
     265        p[k]=p[k]-lead(p[k]);
     266      }
     267      while(deg(lead(p[k]))>(K+1)*deg(var(1))-2*gmsmaxdeg&&
     268        deg(lead(p[k]),v)<=K)
     269      {
     270        q[k]=q[k]+lead(p[k]);
     271        p[k]=p[k]-lead(p[k]);
     272      }
     273    }
     274    q[k]=q[k]+p[k];
     275  }
     276
     277  return(list(r,q));
    229278}
    230279example
    … …  
    248297RETURN:
    249298@format
    250 list l; 
     299list l;
    251300  matrix l[1];  gmsbasis representation of p mod s^(K+1)
    252301  ideal l[2];  p=matrix(gmsbasis)*l[1]+l[2]
    … …  
    282331///////////////////////////////////////////////////////////////////////////////
    283332
    284 static proc min(ideal e)
     333static proc nmin(ideal e)
    285334{
    286335  int i;
    … …  
    297346///////////////////////////////////////////////////////////////////////////////
    298347
    299 static proc max(ideal e)
     348static proc nmax(ideal e)
    300349{
    301350  int i;
    … …  
    348397///////////////////////////////////////////////////////////////////////////////
    349398
    350 static proc basisrep(matrix A0,ideal r,module H,int k0,int K)
     399static proc tjet(matrix A0,ideal r,module H,int k0,int K)
    351400{
    352401  dbprint(printlevel-voice+2,"// compute matrix A of t");
    … …  
    364413///////////////////////////////////////////////////////////////////////////////
    365414
    366 static proc eigvals(matrix A0,ideal r,module H,int k0)
     415static proc eigenval(matrix A0,ideal r,module H,int k0)
    367416{
    368417  dbprint(printlevel-voice+2,
    369418    "// compute eigenvalues e with multiplicities m of A0");
    370419  matrix A;
    371   A,A0,r=basisrep(A0,r,H,k0,0);
     420  A,A0,r=tjet(A0,r,H,k0,0);
    372421  list l=eigenvals(A);
    373422  def e,m=l[1..2];
    … …  
    379428///////////////////////////////////////////////////////////////////////////////
    380429
    381 static proc transf(matrix A,matrix A0,ideal r,module H,module H0,ideal e,intvec m,int k0,int K,int opt)
     430static proc transform(matrix A,matrix A0,ideal r,module H,module H0,ideal e,intvec m,int k0,int K,int opt)
    382431{
    383432  int mu=ncols(gmsbasis);
    384433
    385   number e0,e1=min(e),max(e);
     434  number e0,e1=nmin(e),nmax(e);
    386435
    387436  int i,j,k;
    … …  
    431480  }
    432481
    433   A,A0,r=basisrep(A0,r,H,k0,K+k1);
     482  A,A0,r=tjet(A0,r,H,k0,K+k1);
    434483  module U0=s^k0*freemodule(mu);
    435484
    … …  
    514563ASSUME:   characteristic 0; local degree ordering;
    515564          isolated critical point 0 of t
    516 RETURN:   list l=jordan(M);  Jordan data of monodromy matrix exp(-2*pi*i*M)
    517 SEE ALSO: mondromy_lib, linalg.lib
     565RETURN:
     566@format
     567list l=jordan(M);  Jordan data of monodromy matrix exp(-2*pi*i*M)
     568  ideal l[1];
     569    number l[1][i];  eigenvalue of i-th Jordan block of M
     570  intvec l[2];
     571    int l[2][i];  size of i-th Jordan block of M
     572  intvec l[3];
     573    int l[3][i];  multiplicity of i-th Jordan block of M
     574@end format
     575SEE ALSO: mondromy_lib, linalg_lib
    518576KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; monodromy
    519577EXAMPLE:  example monodromy; shows examples
    … …  
    531589
    532590  def A0,r,H,H0,k0=saturate();
    533   e,m,A0,r=eigvals(A0,r,H,k0);
    534   A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,0,0);
     591  e,m,A0,r=eigenval(A0,r,H,k0);
     592  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,0,0);
    535593
    536594  list l=jordan(A,e,m);
    … …  
    562620@end format
    563621SEE ALSO: spectrum_lib
    564 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice; 
     622KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
    565623          mixed Hodge structure; V-filtration; spectrum
    566624EXAMPLE:  example spectrum; shows examples
    … …  
    589647  intvec spp[2];
    590648    int spp[2][i];  weight filtration index of i-th spectral pair
    591   intvec spp[3]; 
     649  intvec spp[3];
    592650    int spp[3][i];  multiplicity of i-th spectral pair
    593651@end format
    … …  
    610668///////////////////////////////////////////////////////////////////////////////
    611669
    612 proc spnf(ideal e,list #)
    613 "USAGE:   spnf(e[,m][,V]); ideal e, intvec m, list V
    614 ASSUME:  ncols(e)==size(m)==size(V); typeof(V[i])=="int"
    615 RETURN:
    616 @format
    617 list sp;  spectrum normal form of (e,m,V)
    618   ideal sp[1];  numbers in e in increasing order
    619   intvec sp[2];
    620     int sp[2][i];  multiplicity of number sp[1][i] in e
    621   list sp[3];
    622     module sp[3][i];  module associated to number sp[1][i]
    623 @end format
     670proc spnf(ideal a,list #)
     671"USAGE:   spnf(a[,m][,V]); ideal a, intvec m, list V
     672ASSUME:  ncols(a)==size(m)==size(V); typeof(V[i])=="int"
     673RETURN:  order (a[i][,V[i]]) with multiplicity m[i] lexicographically
    624674EXAMPLE: example spnf; shows examples
    625675"
    626676{
    627   int n=ncols(e);
     677  int n=ncols(a);
    628678  intvec m;
    629679  module v;
    … …  
    659709
    660710  int k;
    661   ideal e0;
     711  ideal a0;
    662712  intvec m0;
    663713  list V0;
    664   number e1;
     714  number a1;
    665715  int m1;
    666716  for(i=n;i>=1;i--)
    … …  
    672722        if(m[j]!=0)
    673723        {
    674           if(number(e[i])>number(e[j]))
     724          if(number(a[i])>number(a[j]))
    675725          {
    676             e1=number(e[i]);
    677             e[i]=e[j];
    678             e[j]=e1;
     726            a1=number(a[i]);
     727            a[i]=a[j];
     728            a[j]=a1;
    679729            m1=m[i];
    680730            m[i]=m[j];
    … …  
    687737            }
    688738          }
    689           if(number(e[i])==number(e[j]))
     739          if(number(a[i])==number(a[j]))
    690740          {
    691741            m[i]=m[i]+m[j];
    … …  
    699749      }
    700750      k++;
    701       e0[k]=e[i];
     751      a0[k]=a[i];
    702752      m0[k]=m[i];
    703753      if(size(V)>0)
    … …  
    725775  if(k>0)
    726776  {
    727     l=e0,m0;
     777    l=a0,m0;
    728778    if(size(V0)>0)
    729779    {
    … …  
    739789
    740790proc sppnf(ideal a,intvec w,list #)
    741 "USAGE:    sppnf(a,w[,m][,V]); ideal a, intvec w, intvec m, list V
    742 ASSUME:   ncols(e)=size(w)=size(m)=size(V); typeof(V[i])=="module"
    743 RETURN:
    744 @format
    745 list spp;  spectral pairs normal form of (a,w,m,V)
    746   ideal spp[1];
    747     number spp[1][i];  V-filtration index of i-th spectral pair
    748   intvec spp[2];
    749     int spp[2][i];  weight filtration index of i-th spectral pair
    750   intvec spp[3];
    751     int spp[3][i];  multiplicity of i-th spectral pair
    752   list spp[4];
    753     module spp[4][i];  vector space of i-th spectral pair
    754 @end format
    755 SEE ALSO: spectrum_lib
    756 KEYWORDS: singularities; Gauss-Manin connection; Brieskorn lattice;
    757           mixed Hodge structure; V-filtration; weight filtration;
    758           spectrum; spectral pairs
    759 EXAMPLE:  example sppnorm; shows examples
     791"USAGE:   sppnf(a,w[,m][,V]); ideal a, intvec w, intvec m, list V
     792ASSUME:  ncols(e)=size(w)=size(m)=size(V); typeof(V[i])=="module"
     793RETURN:  order (a[i][,w[i]][,V[i]]) with multiplicity m[i] lexicographically
     794EXAMPLE: example sppnorm; shows examples
    760795"
    761796{
    … …  
    806841      {
    807842        if(m[j]!=0)
    808         {
     843        {
    809844          if(number(a[i])>number(a[j])||
    810845            (number(a[i])==number(a[j])&&w[i]<w[j]))
    … …  
    887922  ideal v[1];
    888923    number v[1][i];  V-filtration index of i-th spectral pair
    889   intvec v[2]; 
     924  intvec v[2];
    890925    int v[2][i];  multiplicity of i-th spectral pair
    891   list v[3]; 
     926  list v[3];
    892927    module v[3][i];  vector space of i-th graded part in terms of v[4]
    893928  ideal v[4];  monomial vector space basis of H''/s*H''
    … …  
    922957  intvec vw[2];
    923958    int vw[2][i];  weight filtration index of i-th spectral pair
    924   intvec vw[3]; 
     959  intvec vw[3];
    925960    int vw[3][i];  multiplicity of i-th spectral pair
    926   list vw[4]; 
     961  list vw[4];
    927962    module vw[4][i];  vector space of i-th graded part in terms of vw[5]
    928963  ideal vw[5];  monomial vector space basis of H''/s*H''
    … …  
    948983
    949984  def A0,r,H,H0,k0=saturate();
    950   e,m,A0,r=eigvals(A0,r,H,k0);
    951   A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,0,1);
     985  e,m,A0,r=eigenval(A0,r,H,k0);
     986  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,0,1);
    952987
    953988  dbprint(printlevel-voice+2,"// compute weight filtration basis");
    … …  
    10521087
    10531088  def A0,r,H,H0,k0=saturate();
    1054   e,m,A0,r=eigvals(A0,r,H,k0);
    1055   int k1=int(max(e)-min(e));
    1056   A,A0,r,H0,U0,e,m=transf(A,A0,r,H,H0,e,m,k0,k0+k1,1);
     1089  e,m,A0,r=eigenval(A0,r,H,k0);
     1090  int k1=int(nmax(e)-nmin(e));
     1091  A,A0,r,H0,U0,e,m=transform(A,A0,r,H,H0,e,m,k0,k0+k1,1);
    10571092
    10581093  ring S=0,s,(ds,c);
    … …  
    12281263  ideal ev[1];
    12291264    number ev[1][i];  V-filtration index of i-th spectral pair
    1230   intvec ev[2]; 
     1265  intvec ev[2];
    12311266    int ev[2][i];  multiplicity of i-th spectral pair
    1232   list ev[3]; 
     1267  list ev[3];
    12331268    module ev[3][i];  vector space of i-th graded part in terms of ev[4]
    12341269  ideal ev[4];  monomial vector space basis of Jacobian algebra
    … …  
    12551290  for(i=ncols(m);i>=1;i--)
    12561291  {
    1257     M[i]=division(V0,coeffs(reduce(m[i]*m,U,g),m)*V0)[1];
     1292    M[i]=division(V0,coeffs(reduce(m[i]*m,g,U),m)*V0)[1];
    12581293  }
    12591294
    … …  
    16171652@format
    16181653list l;
    1619   intvec l[i];  if the spectra sp0 occur with multiplicities k 
    1620                 in a deformation of a [quasihomogeneous] singularity 
     1654  intvec l[i];  if the spectra sp0 occur with multiplicities k
     1655                in a deformation of a [quasihomogeneous] singularity
    16211656                with spectrum sp then k<=l[i]
    16221657@end format
    … …  
    16341669      {
    16351670        if(size(l)>0)
    1636         {
     1671        {
    16371672          j=1;
    16381673          while(j<size(l)&&l[j]!=l0[i])
    1639           {
     1674          {
    16401675            j++;
    16411676          }
    16421677          if(l[j]==l0[i])
    1643           {
     1678          {
    16441679            l[j][size(sp)]=k;
    16451680          }
    16461681          else
    1647           {
     1682          {
    16481683            l0[i][size(sp)]=k;
    16491684            l=l+list(l0[i]);
    16501685          }
    1651         }
     1686        }
    16521687        else
    1653         {
     1688        {
    16541689          l=l0;
    1655         }
     1690        }
    16561691      }
    16571692    }

  • Singular/LIB/linalg.lib

    rfee17e re9124e  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.25 2002-02-16 18:26:08 mschulze Exp $";
     3version="$Id: linalg.lib,v 1.26 2002-03-06 16:38:54 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 hessenberg(M);       Hessenberg form of M
     28 evnf(e[,m]);         eigenvalues normal form of (e[,m])
    2729 eigenvals(M);        eigenvalues and multiplicities of M
    2830 jordan(M);           Jordan data of M
    … …  
    12351237///////////////////////////////////////////////////////////////////////////////
    12361238
     1239static proc rowcolswap(matrix M,int i,int j)
     1240{
     1241  if(i==j)
     1242  {
     1243    return(M);
     1244  }
     1245  poly p;
     1246  for(int k=1;k<=nrows(M);k++)
     1247  {
     1248    p=M[i,k];
     1249    M[i,k]=M[j,k];
     1250    M[j,k]=p;
     1251  }
     1252  for(k=1;k<=ncols(M);k++)
     1253  {
     1254    p=M[k,i];
     1255    M[k,i]=M[k,j];
     1256    M[k,j]=p;
     1257  }
     1258  return(M);
     1259}
     1260//////////////////////////////////////////////////////////////////////////////
     1261
     1262static proc rowelim(matrix M,int i,int j,int k)
     1263{
     1264  if(jet(M[i,k],0)==0||jet(M[j,k],0)==0)
     1265  {
     1266    return(M);
     1267  }
     1268  number n=number(jet(M[i,k],0))/number(jet(M[j,k],0));
     1269  for(int l=1;l<=ncols(M);l++)
     1270  {
     1271    M[i,l]=M[i,l]-n*M[j,l];
     1272  }
     1273  for(l=1;l<=nrows(M);l++)
     1274  {
     1275    M[l,j]=M[l,j]+n*M[l,i];
     1276  }
     1277  return(M);
     1278}
     1279///////////////////////////////////////////////////////////////////////////////
     1280
     1281static proc colelim(matrix M,int i,int j,int k)
     1282{
     1283  if(jet(M[k,i],0)==0||jet(M[k,j],0)==0)
     1284  {
     1285    return(M);
     1286  }
     1287  number n=number(jet(M[k,i],0))/number(jet(M[k,j],0));
     1288  for(int l=1;l<=nrows(M);l++)
     1289  {
     1290    M[l,i]=M[l,i]-n*M[l,j];
     1291  }
     1292  for(l=1;l<=ncols(M);l++)
     1293  {
     1294    M[j,l]=M[j,l]+n*M[i,l];
     1295  }
     1296  return(M);
     1297}
     1298///////////////////////////////////////////////////////////////////////////////
     1299
     1300proc hessenberg(matrix M)
     1301"USAGE:   hessenberg(M); matrix M
     1302ASSUME:  M constant square matrix
     1303RETURN:  matrix H;  Hessenberg form of M
     1304EXAMPLE: example hessenberg; shows examples
     1305"
     1306{
     1307  if(system("with","eigenval"))
     1308  {
     1309   return(system("hessenberg",M));
     1310  }
     1311
     1312  int n=ncols(M);
     1313  int i,j;
     1314  for(int k=1;k<n-1;k++)
     1315  {
     1316    j=k+1;
     1317    while(j<n&&jet(M[j,k],0)==0)
     1318    {
     1319      j++;
     1320    }
     1321    if(jet(M[j,k],0)!=0)
     1322    {
     1323      M=rowcolswap(M,j,k+1);
     1324      for(i=j+1;i<=n;i++)
     1325      {
     1326        M=rowelim(M,i,k+1,k);
     1327      }
     1328    }
     1329  }
     1330  return(M);
     1331}
     1332example
     1333{ "EXAMPLE:"; echo=2;
     1334  ring R=0,x,dp;
     1335  matrix M[3][3]=3,2,1,0,2,1,0,0,3;
     1336  print(M);
     1337  print(hessenberg(M));
     1338}
     1339///////////////////////////////////////////////////////////////////////////////
     1340
     1341proc evnf(ideal e,list #)
     1342"USAGE:   evnf(e[,m]); ideal e, intvec m
     1343ASSUME:  ncols(e)==size(m)
     1344RETURN:  order eigenvalues e with multiplicities m
     1345EXAMPLE: example evnf; shows examples
     1346"
     1347{
     1348  int n=ncols(e);
     1349  intvec m;
     1350  int i,j;
     1351  while(i<size(#))
     1352  {
     1353    i++;
     1354    if(typeof(#[i])=="intvec")
     1355    {
     1356      m=#[i];
     1357    }
     1358  }
     1359  if(m==0)
     1360  {
     1361    for(i=n;i>=1;i--)
     1362    {
     1363      m[i]=1;
     1364    }
     1365  }
     1366
     1367  int k;
     1368  ideal e0;
     1369  intvec m0;
     1370  number e1;
     1371  int m1;
     1372  for(i=n;i>=1;i--)
     1373  {
     1374    if(m[i]!=0)
     1375    {
     1376      for(j=i-1;j>=1;j--)
     1377      {
     1378        if(m[j]!=0)
     1379        {
     1380          if(number(e[i])>number(e[j]))
     1381          {
     1382            e1=number(e[i]);
     1383            e[i]=e[j];
     1384            e[j]=e1;
     1385            m1=m[i];
     1386            m[i]=m[j];
     1387            m[j]=m1;
     1388          }
     1389          if(number(e[i])==number(e[j]))
     1390          {
     1391            m[i]=m[i]+m[j];
     1392            m[j]=0;
     1393          }
     1394        }
     1395      }
     1396      k++;
     1397      e0[k]=e[i];
     1398      m0[k]=m[i];
     1399    }
     1400  }
     1401
     1402  list l;
     1403  if(k>0)
     1404  {
     1405    l=e0,m0;
     1406  }
     1407  return(l);
     1408}
     1409example
     1410{ "EXAMPLE:"; echo=2;
     1411}
     1412///////////////////////////////////////////////////////////////////////////////
     1413
    12371414proc eigenvals(matrix M)
    12381415"USAGE:   eigenvals(M); matrix M
    … …  
    12491426"
    12501427{
    1251   return(system("eigenvals",M));
     1428  if(system("with","eigenval"))
     1429  {
     1430   return(system("eigenvals",M));
     1431  }
     1432
     1433  M=jet(hessenberg(M),0);
     1434  int n=ncols(M);
     1435  int k;
     1436  ideal e;
     1437  intvec m;
     1438  number e0;
     1439  intvec v;
     1440  list l;
     1441  int i,j;
     1442  j=1;
     1443  while(j<=n)
     1444  {
     1445    v=j;
     1446    j++;
     1447    if(j<=n)
     1448    {
     1449      while(j<n&&M[j,j-1]!=0)
     1450      {
     1451        v=v,j;
     1452        j++;
     1453      }
     1454      if(M[j,j-1]!=0)
     1455      {
     1456        v=v,j;
     1457        j++;
     1458      }
     1459    }
     1460    if(size(v)==1)
     1461    {
     1462      k++;
     1463      e[k]=M[v,v];
     1464      m[k]=1;
     1465    }
     1466    else
     1467    {
     1468      l=factorize(det(submat(M,v,v)-var(1)));
     1469      for(i=size(l[1]);i>=1;i--)
     1470      {
     1471        e0=number(jet(l[1][i]/var(1),0));
     1472        if(e0!=0)
     1473        {
     1474          k++;
     1475          e[k]=(e0*var(1)-l[1][i])/e0;
     1476          m[k]=l[2][i];
     1477        }
     1478      }
     1479    }
     1480  }
     1481  return(evnf(e,m));
    12521482}
    12531483example
    … …  
    15171747  print(jordannf(M));
    15181748}
     1749
    15191750///////////////////////////////////////////////////////////////////////////////
    15201751
Note: See TracChangeset for help on using the changeset viewer.