Changeset a0c62d in git for Singular/LIB/linalg.lib


Ignore:
Timestamp:
Aug 25, 2001, 4:56:41 PM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
Children:
e48054420f3b08092bdc82ec27f35dadd698c988
Parents:
0ff6b53f29ede7e671edc683ddee8ffff811f123
Message:
*mschulze: new proc spnormalize


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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/linalg.lib

    r0ff6b5 ra0c62d  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.18 2001-08-25 10:51:07 mschulze Exp $";
     3version="$Id: linalg.lib,v 1.19 2001-08-25 14:56:41 mschulze Exp $";
    44category="Linear Algebra";
    55info="
    … …  
    99
    1010PROCEDURES:
    11  inverse(A);        matrix, the inverse of A
    12  inverse_B(A);      list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
    13  inverse_L(A);      list(matrix Inv,poly p),Inv*A=p*En ( using lift )
    14  sym_gauss(A);      symmetric gaussian algorithm
    15  orthogonalize(A);  Gram-Schmidt orthogonalization
    16  diag_test(A);      test whether A can be diagnolized
    17  busadj(A);         coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
    18  charpoly(A,v);     characteristic polynomial of A ( using busadj(A) )
    19  adjoint(A);        adjoint of A ( using busadj(A) )
    20  det_B(A);          determinant of A ( using busadj(A) )
    21  gaussred(A);       gaussian reduction: P*A=U*S, S a row reduced form of A
    22  gaussred_pivot(A); gaussian reduction: P*A=U*S, uses row pivoting
    23  gauss_nf(A);       gaussian normal form of A
    24  mat_rk(A);         rank of constant matrix A
    25  U_D_O(A);          P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
    26  pos_def(A,i);      test symmetric matrix for positive definiteness
    27  hessenberg(M);     transforms M to Hessenberg form
    28  eigenval(M);       eigenvalues of M with multiplicities
    29  jordan(M);         Jordan data of constant square matrix M
    30  jordanbasis(M);    Jordan basis of constant square matrix M
    31  jordanmatrix(e,b); Jordan matrix with eigenvalues e and Jordan block sizes b
    32  jordanform(M);     Jordan matrix of constant square matrix M
    33  commutator(A);     matrix of commutator B->[A,B]=AB-BA
     11 inverse(A);          matrix, the inverse of A
     12 inverse_B(A);        list(matrix Inv,poly p),Inv*A=p*En ( using busadj(A) )
     13 inverse_L(A);        list(matrix Inv,poly p),Inv*A=p*En ( using lift )
     14 sym_gauss(A);        symmetric gaussian algorithm
     15 orthogonalize(A);    Gram-Schmidt orthogonalization
     16 diag_test(A);        test whether A can be diagnolized
     17 busadj(A);           coefficients of Adj(E*t-A) and coefficients of det(E*t-A)
     18 charpoly(A,v);       characteristic polynomial of A ( using busadj(A) )
     19 adjoint(A);          adjoint of A ( using busadj(A) )
     20 det_B(A);            determinant of A ( using busadj(A) )
     21 gaussred(A);         gaussian reduction: P*A=U*S, S a row reduced form of A
     22 gaussred_pivot(A);   gaussian reduction: P*A=U*S, uses row pivoting
     23 gauss_nf(A);         gaussian normal form of A
     24 mat_rk(A);           rank of constant matrix A
     25 U_D_O(A);            P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang
     26 pos_def(A,i);        test symmetric matrix for positive definiteness
     27 hessenberg(M);       transforms M to Hessenberg form
     28 spnormalize(a[,m]);  normalize spectrum
     29 eigenval(M);         eigenvalues of M with multiplicities
     30 jordan(M);           Jordan data of constant square matrix M
     31 jordanbasis(M);      Jordan basis of constant square matrix M
     32 jordanmatrix(e,b);   Jordan matrix with eigenvalues e and Jordan block sizes b
     33 jordanform(M);       Jordan matrix of constant square matrix M
     34 commutator(A);       matrix of commutator B->[A,B]=AB-BA
    3435";
    3536
    … …  
    13341335///////////////////////////////////////////////////////////////////////////////
    13351336
    1336 static proc addval(list l,poly e,int m)
     1337static proc spappend(list l,number e,int m)
    13371338{
    13381339  if(size(l)==0)
    13391340  {
    1340     return(list(ideal(e),intvec(m)));
    1341   }
    1342   int n=ncols(l[1]);
    1343   for(int i=n;i>=1;i--)
    1344   {
    1345     if(l[1][i]==e)
     1341    l=list(ideal(e),intvec(m));
     1342  }
     1343  else
     1344  {
     1345    int n=ncols(l[1]);
     1346    l[1][n+1]=e;
     1347    l[2][n+1]=m;
     1348  }
     1349  return(l);
     1350}
     1351///////////////////////////////////////////////////////////////////////////////
     1352
     1353proc spnormalize(ideal e,list #)
     1354"USAGE:    spnormalize(a,w[,m]);
     1355RETURN:
     1356@format
     1357list Sp: normalized spectrum (a,m)
     1358  ideal Sp[1]: numbers in a in increasing order
     1359  intvec Sp[2]:
     1360    int Sp[2][i]: multiplicity of number Sp[1][i] in a
     1361@end format
     1362EXAMPLE:  example spnormalize; shows examples
     1363"
     1364{
     1365  intvec m;
     1366  int i,j;
     1367  if(size(#)==0)
     1368  {
     1369    for(i=ncols(e);i>=1;i--)
    13461370    {
    1347       l[2][i]=l[2][i]+m;
    1348       return(l);
    1349     }
    1350   }
    1351   l[1][n+1]=e;
    1352   l[2][n+1]=m;
     1371      m[i]=1;
     1372    }
     1373  }
     1374  else
     1375  {
     1376    m=#[1];
     1377  }
     1378
     1379  list l;
     1380  number e0;
     1381  int m0;
     1382  for(i=ncols(e);i>=1;i--)
     1383  {
     1384    if(m[i]!=0)
     1385    {
     1386      for(j=i-1;j>=1;j--)
     1387      {
     1388        if(m[j]!=0)
     1389        {
     1390          if(number(e[i])>number(e[j]))
     1391          {
     1392            e0=number(e[i]);
     1393            e[i]=e[j];
     1394            e[j]=e0;
     1395            m0=m[i];
     1396            m[i]=m[j];
     1397            m[j]=m0;
     1398          }
     1399          if(number(e[i])==number(e[j]))
     1400          {
     1401            m[i]=m[i]+m[j];
     1402            m[j]=0;
     1403          }
     1404        }
     1405      }
     1406      l=spappend(l,number(e[i]),m[i]);
     1407    }
     1408  }
     1409
    13531410  return(l);
    13541411}
    1355 ///////////////////////////////////////////////////////////////////////////////
    1356 
    1357 static proc sortval(list l)
    1358 {
    1359   def e,m=l[1..2];
    1360   int n=ncols(e);
    1361   poly ee;
    1362   int mm;
    1363   int i,j;
    1364   for(i=n;i>1;i--)
    1365   {
    1366     for(j=i-1;j>=1;j--)
    1367     {
    1368       if(number(e[j])>number(e[i]))
    1369       {
    1370         ee=e[i];
    1371         e[i]=e[j];
    1372         e[j]=ee;
    1373         mm=m[i];
    1374         m[i]=m[j];
    1375         m[j]=mm;
    1376       }
    1377     }
    1378   }
    1379   return(list(e,m));
     1412example
     1413{ "EXAMPLE:"; echo=2;
    13801414}
    13811415///////////////////////////////////////////////////////////////////////////////
    … …  
    13961430  M=jet(hessenberg(M),0);
    13971431  int n=ncols(M);
    1398   number e;
     1432  number e0;
    13991433  intvec v;
    14001434  list l,f;
    1401   int i;
    1402   int j=1;
     1435  int i,j,k;
     1436  j=1;
    14031437  while(j<=n)
    14041438  {
    … …  
    14201454    if(size(v)==1)
    14211455    {
    1422       l=addval(l,M[v,v],1);
     1456      l=spappend(l,number(M[v,v]),1);
    14231457    }
    14241458    else
    … …  
    14271461      for(i=size(f[1]);i>=1;i--)
    14281462      {
    1429         e=number(jet(f[1][i]/var(1),0));
    1430         if(e!=0)
     1463        e0=number(jet(f[1][i]/var(1),0));
     1464        if(e0!=0)
    14311465        {
    1432           l=addval(l,(e*var(1)-f[1][i])/e,f[2][i]);
     1466          l=spappend(l,number((e0*var(1)-f[1][i])/e0),f[2][i]);
    14331467        }
    14341468      }
    14351469    }
    14361470  }
    1437   return(sortval(l));
     1471  return(spnormalize(l[1],l[2]));
    14381472}
    14391473example
    … …  
    14461480///////////////////////////////////////////////////////////////////////////////
    14471481
    1448 proc jordan(matrix M)
     1482proc jordan(matrix M,list #)
    14491483"USAGE:   jordan(M); matrix M
    14501484ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
    … …  
    14701504  M=jet(M,0);
    14711505
    1472   list l=eigenvalues(M);
    1473   def e0,m0=l[1..2];
    1474   kill l;
     1506  if(size(#)==0)
     1507  {
     1508    #=eigenvalues(M);
     1509  }
     1510  def e0,m0=#[1..2];
    14751511
    14761512  int i;
    … …  
    15351571///////////////////////////////////////////////////////////////////////////////
    15361572
    1537 proc jordanbasis(matrix M)
     1573proc jordanbasis(matrix M,list #)
    15381574"USAGE:   jordanbasis(M); matrix M
    15391575ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
    … …  
    15581594  M=jet(M,0);
    15591595
    1560   list l=eigenvalues(M);
    1561   def e,m=l[1..2];
    1562   kill l;
     1596  if(size(#)==0)
     1597  {
     1598    #=eigenvalues(M);
     1599  }
     1600  def e,m=#[1..2];
    15631601
    15641602  int i;
    … …  
    16771715///////////////////////////////////////////////////////////////////////////////
    16781716
    1679 proc jordanform(matrix M)
     1717proc jordanform(matrix M,list #)
    16801718"USAGE:   jordanform(M); matrix M
    16811719ASSUME:  M constant square matrix, eigenvalues of M in coefficient field
    … …  
    16841722"
    16851723{
    1686   list l=jordan(M);
     1724  list l=jordan(M,#);
    16871725  return(jordanmatrix(l[1],l[2]));
    16881726}
Note: See TracChangeset for help on using the changeset viewer.