Changeset e5ae343 in git


Ignore:
Timestamp:
Feb 16, 2002, 11:58:39 AM (22 years ago)
Author:
Mathias Schulze <mschulze@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b52fc4b2495505785981d640dcf7eb3e456778ef')
Children:
adf5bd69bf9afcfb41fbd46edfe3f594d586b715
Parents:
d9d2417dd93f0480a9af34fb65f29dd5ac08309c
Message:
*mschulze: moved spnf to gaussman.lib, eigenvalues to kernel


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/linalg.lib

    rd9d241 re5ae343  
    11//GMG last modified: 04/25/2000
    22//////////////////////////////////////////////////////////////////////////////
    3 version="$Id: linalg.lib,v 1.22 2002-01-22 15:27:42 mschulze Exp $";
     3version="$Id: linalg.lib,v 1.23 2002-02-16 10:58:39 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  spnf(a[,m][,V]);     spectrum normal form of (a,m,V)
    2927 eigenvals(M);        eigenvalues and multiplicities of M
    3028 jordan(M);           Jordan data of M
     
    12371235///////////////////////////////////////////////////////////////////////////////
    12381236
    1239 static proc swap(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 
    1262 static 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 
    1281 static 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]-m*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 
    1300 proc hessenberg(matrix M)
    1301 "USAGE:   hessenberg(M); matrix M
    1302 ASSUME:  M constant square matrix
    1303 RETURN:  matrix H;  Hessenberg form of M
    1304 EXAMPLE: example hessenberg; shows examples
    1305 "
    1306 {
    1307   int n=ncols(M);
    1308   int i,j;
    1309   for(int k=1;k<n-1;k++)
    1310   {
    1311     j=k+1;
    1312     while(j<n&&jet(M[j,k],0)==0)
    1313     {
    1314       j++;
    1315     }
    1316     if(jet(M[j,k],0)!=0)
    1317     {
    1318       M=swap(M,j,k+1);
    1319       for(i=j+1;i<=n;i++)
    1320       {
    1321         M=rowelim(M,i,k+1,k);
    1322       }
    1323     }
    1324   }
    1325   return(M);
    1326 }
    1327 example
    1328 { "EXAMPLE:"; echo=2;
    1329   ring R=0,x,dp;
    1330   matrix M[3][3]=3,2,1,0,2,1,0,0,3;
    1331   print(M);
    1332   print(hessenberg(M));
    1333 }
    1334 ///////////////////////////////////////////////////////////////////////////////
    1335 
    1336 proc spnf(ideal e,list #)
    1337 "USAGE:   spnf(e[,m][,V]); ideal e, intvec m, list V
    1338 ASSUME:  ncols(e)==size(m)==size(V); typeof(V[i])=="int"
    1339 RETURN:
    1340 @format
    1341 list sp;  spectrum normal form of (e,m,V)
    1342   ideal sp[1];  numbers in e in increasing order
    1343   intvec sp[2];
    1344     int sp[2][i];  multiplicity of number sp[1][i] in e
    1345   list sp[3];
    1346     module sp[3][i];  module associated to number sp[1][i]
    1347 @end format
    1348 EXAMPLE: example spnf; shows examples
    1349 "
    1350 {
    1351   int n=ncols(e);
    1352   intvec m;
    1353   module v;
    1354   list V;
    1355   int i,j;
    1356   while(i<size(#))
    1357   {
    1358     i++;
    1359     if(typeof(#[i])=="intvec")
    1360     {
    1361       m=#[i];
    1362     }
    1363     if(typeof(#[i])=="module")
    1364     {
    1365       v=#[i];
    1366       for(j=n;j>=1;j--)
    1367       {
    1368         V[j]=module(v[j]);
    1369       }
    1370     }
    1371     if(typeof(#[i])=="list")
    1372     {
    1373       V=#[i];
    1374     }
    1375   }
    1376   if(m==0)
    1377   {
    1378     for(i=n;i>=1;i--)
    1379     {
    1380       m[i]=1;
    1381     }
    1382   }
    1383 
    1384   int k;
    1385   ideal e0;
    1386   intvec m0;
    1387   list V0;
    1388   number e1;
    1389   int m1;
    1390   for(i=n;i>=1;i--)
    1391   {
    1392     if(m[i]!=0)
    1393     {
    1394       for(j=i-1;j>=1;j--)
    1395       {
    1396         if(m[j]!=0)
    1397         {
    1398           if(number(e[i])>number(e[j]))
    1399           {
    1400             e1=number(e[i]);
    1401             e[i]=e[j];
    1402             e[j]=e1;
    1403             m1=m[i];
    1404             m[i]=m[j];
    1405             m[j]=m1;
    1406             if(size(V)>0)
    1407             {
    1408               v=V[i];
    1409               V[i]=V[j];
    1410               V[j]=v;
    1411             }
    1412           }
    1413           if(number(e[i])==number(e[j]))
    1414           {
    1415             m[i]=m[i]+m[j];
    1416             m[j]=0;
    1417             if(size(V)>0)
    1418             {
    1419               V[i]=V[i]+V[j];
    1420             }
    1421           }
    1422         }
    1423       }
    1424       k++;
    1425       e0[k]=e[i];
    1426       m0[k]=m[i];
    1427       if(size(V)>0)
    1428       {
    1429         V0[k]=V[i];
    1430       }
    1431     }
    1432   }
    1433 
    1434   if(size(V0)>0)
    1435   {
    1436     n=size(V0);
    1437     module U=std(V0[n]);
    1438     for(i=n-1;i>=1;i--)
    1439     {
    1440       V0[i]=simplify(reduce(V0[i],U),1);
    1441       if(i>=2)
    1442       {
    1443         U=std(U+V0[i]);
    1444       }
    1445     }
    1446   }
    1447 
    1448   list l;
    1449   if(k>0)
    1450   {
    1451     l=e0,m0;
    1452     if(size(V0)>0)
    1453     {
    1454       l[3]=V0;
    1455     }
    1456   }
    1457   return(l);
    1458 }
    1459 example
    1460 { "EXAMPLE:"; echo=2;
    1461 }
    1462 ///////////////////////////////////////////////////////////////////////////////
    1463 
    14641237proc eigenvals(matrix M)
    14651238"USAGE:   eigenvals(M); matrix M
     
    14761249"
    14771250{
    1478   M=jet(hessenberg(M),0);
    1479   int n=ncols(M);
    1480   int k;
    1481   ideal e;
    1482   intvec m;
    1483   number e0;
    1484   intvec v;
    1485   list l;
    1486   int i,j;
    1487   j=1;
    1488   while(j<=n)
    1489   {
    1490     v=j;
    1491     j++;
    1492     if(j<=n)
    1493     {
    1494       while(j<n&&M[j,j-1]!=0)
    1495       {
    1496         v=v,j;
    1497         j++;
    1498       }
    1499       if(M[j,j-1]!=0)
    1500       {
    1501         v=v,j;
    1502         j++;
    1503       }
    1504     }
    1505     if(size(v)==1)
    1506     {
    1507       k++;
    1508       e[k]=M[v,v];
    1509       m[k]=1;
    1510     }
    1511     else
    1512     {
    1513       l=factorize(det(submat(M,v,v)-var(1)));
    1514       for(i=size(l[1]);i>=1;i--)
    1515       {
    1516         e0=number(jet(l[1][i]/var(1),0));
    1517         if(e0!=0)
    1518         {
    1519           k++;
    1520           e[k]=(e0*var(1)-l[1][i])/e0;
    1521           m[k]=l[2][i];
    1522         }
    1523       }
    1524     }
    1525   }
    1526   return(spnf(e,m));
     1251  def R=basering;
     1252  ring P=0,t,dp;
     1253  map zero=R,0;
     1254  matrix M=zero(M);
     1255  list l=system("eigenvalues",M);
     1256  setring R;
     1257  map zero=P,0;
     1258  return(zero(l));
    15271259}
    15281260example
Note: See TracChangeset for help on using the changeset viewer.