Changeset e5ae343 in git
 Timestamp:
 Feb 16, 2002, 11:58:39 AM (22 years ago)
 Branches:
 (u'fiekerDuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b52fc4b2495505785981d640dcf7eb3e456778ef')
 Children:
 adf5bd69bf9afcfb41fbd46edfe3f594d586b715
 Parents:
 d9d2417dd93f0480a9af34fb65f29dd5ac08309c
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/linalg.lib
rd9d241 re5ae343 1 1 //GMG last modified: 04/25/2000 2 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1.2 2 20020122 15:27:42mschulze Exp $";3 version="$Id: linalg.lib,v 1.23 20020216 10:58:39 mschulze Exp $"; 4 4 category="Linear Algebra"; 5 5 info=" … … 25 25 U_D_O(A); P*A=U*D*O, P,D,U,O=permutaion,diag,lower,uppertriang 26 26 pos_def(A,i); test symmetric matrix for positive definiteness 27 hessenberg(M); Hessenberg form of M28 spnf(a[,m][,V]); spectrum normal form of (a,m,V)29 27 eigenvals(M); eigenvalues and multiplicities of M 30 28 jordan(M); Jordan data of M … … 1237 1235 /////////////////////////////////////////////////////////////////////////////// 1238 1236 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)==0jet(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)==0jet(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 M1302 ASSUME: M constant square matrix1303 RETURN: matrix H; Hessenberg form of M1304 EXAMPLE: example hessenberg; shows examples1305 "1306 {1307 int n=ncols(M);1308 int i,j;1309 for(int k=1;k<n1;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 example1328 { "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 V1338 ASSUME: ncols(e)==size(m)==size(V); typeof(V[i])=="int"1339 RETURN:1340 @format1341 list sp; spectrum normal form of (e,m,V)1342 ideal sp[1]; numbers in e in increasing order1343 intvec sp[2];1344 int sp[2][i]; multiplicity of number sp[1][i] in e1345 list sp[3];1346 module sp[3][i]; module associated to number sp[1][i]1347 @end format1348 EXAMPLE: example spnf; shows examples1349 "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=i1;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=n1;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 example1460 { "EXAMPLE:"; echo=2;1461 }1462 ///////////////////////////////////////////////////////////////////////////////1463 1464 1237 proc eigenvals(matrix M) 1465 1238 "USAGE: eigenvals(M); matrix M … … 1476 1249 " 1477 1250 { 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,j1]!=0) 1495 { 1496 v=v,j; 1497 j++; 1498 } 1499 if(M[j,j1]!=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)); 1527 1259 } 1528 1260 example
Note: See TracChangeset
for help on using the changeset viewer.