Changeset a0c62d in git for Singular/LIB/linalg.lib
 Timestamp:
 Aug 25, 2001, 4:56:41 PM (22 years ago)
 Branches:
 (u'spielwiese', '8e0ad00ce244dfd0756200662572aef8402f13d5')
 Children:
 e48054420f3b08092bdc82ec27f35dadd698c988
 Parents:
 0ff6b53f29ede7e671edc683ddee8ffff811f123
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

Singular/LIB/linalg.lib
r0ff6b5 ra0c62d 1 1 //GMG last modified: 04/25/2000 2 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1.1 8 20010825 10:51:07mschulze Exp $";3 version="$Id: linalg.lib,v 1.19 20010825 14:56:41 mschulze Exp $"; 4 4 category="Linear Algebra"; 5 5 info=" … … 9 9 10 10 PROCEDURES: 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); GramSchmidt orthogonalization 16 diag_test(A); test whether A can be diagnolized 17 busadj(A); coefficients of Adj(E*tA) and coefficients of det(E*tA) 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,uppertriang 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]=ABBA 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); GramSchmidt orthogonalization 16 diag_test(A); test whether A can be diagnolized 17 busadj(A); coefficients of Adj(E*tA) and coefficients of det(E*tA) 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,uppertriang 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]=ABBA 34 35 "; 35 36 … … 1334 1335 /////////////////////////////////////////////////////////////////////////////// 1335 1336 1336 static proc addval(list l,polye,int m)1337 static proc spappend(list l,number e,int m) 1337 1338 { 1338 1339 if(size(l)==0) 1339 1340 { 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 1353 proc spnormalize(ideal e,list #) 1354 "USAGE: spnormalize(a,w[,m]); 1355 RETURN: 1356 @format 1357 list 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 1362 EXAMPLE: 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) 1346 1370 { 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=i1;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 1353 1410 return(l); 1354 1411 } 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=i1;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)); 1412 example 1413 { "EXAMPLE:"; echo=2; 1380 1414 } 1381 1415 /////////////////////////////////////////////////////////////////////////////// … … 1396 1430 M=jet(hessenberg(M),0); 1397 1431 int n=ncols(M); 1398 number e ;1432 number e0; 1399 1433 intvec v; 1400 1434 list l,f; 1401 int i ;1402 intj=1;1435 int i,j,k; 1436 j=1; 1403 1437 while(j<=n) 1404 1438 { … … 1420 1454 if(size(v)==1) 1421 1455 { 1422 l= addval(l,M[v,v],1);1456 l=spappend(l,number(M[v,v]),1); 1423 1457 } 1424 1458 else … … 1427 1461 for(i=size(f[1]);i>=1;i) 1428 1462 { 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) 1431 1465 { 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]); 1433 1467 } 1434 1468 } 1435 1469 } 1436 1470 } 1437 return(s ortval(l));1471 return(spnormalize(l[1],l[2])); 1438 1472 } 1439 1473 example … … 1446 1480 /////////////////////////////////////////////////////////////////////////////// 1447 1481 1448 proc jordan(matrix M )1482 proc jordan(matrix M,list #) 1449 1483 "USAGE: jordan(M); matrix M 1450 1484 ASSUME: M constant square matrix, eigenvalues of M in coefficient field … … 1470 1504 M=jet(M,0); 1471 1505 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]; 1475 1511 1476 1512 int i; … … 1535 1571 /////////////////////////////////////////////////////////////////////////////// 1536 1572 1537 proc jordanbasis(matrix M )1573 proc jordanbasis(matrix M,list #) 1538 1574 "USAGE: jordanbasis(M); matrix M 1539 1575 ASSUME: M constant square matrix, eigenvalues of M in coefficient field … … 1558 1594 M=jet(M,0); 1559 1595 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]; 1563 1601 1564 1602 int i; … … 1677 1715 /////////////////////////////////////////////////////////////////////////////// 1678 1716 1679 proc jordanform(matrix M )1717 proc jordanform(matrix M,list #) 1680 1718 "USAGE: jordanform(M); matrix M 1681 1719 ASSUME: M constant square matrix, eigenvalues of M in coefficient field … … 1684 1722 " 1685 1723 { 1686 list l=jordan(M );1724 list l=jordan(M,#); 1687 1725 return(jordanmatrix(l[1],l[2])); 1688 1726 }
Note: See TracChangeset
for help on using the changeset viewer.