Changeset 6188357 in git
- Timestamp:
- Jan 3, 2001, 2:24:49 AM (23 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- 9cc5d37df563ce9febe8a449f39cf9986d705e1e
- Parents:
- cda1ad905ad654f1520fe50c63f3628c521558b0
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/linalg.lib
rcda1ad r6188357 1 1 //GMG last modified: 04/25/2000 2 2 ////////////////////////////////////////////////////////////////////////////// 3 version="$Id: linalg.lib,v 1. 7 2000-12-22 14:12:32greuel Exp $";3 version="$Id: linalg.lib,v 1.8 2001-01-03 01:24:49 greuel Exp $"; 4 4 category="Linear Algebra"; 5 5 info=" 6 6 LIBRARY: linalg.lib Algorithmic Linear Algebra 7 AUTHOR: Ivor Saynisch (ivs@math.tu-cottbus.de) 7 AUTHORS: Ivor Saynisch (ivs@math.tu-cottbus.de) 8 @* Mathias Schulze (mschulze@mathematik.uni-kl.de) 8 9 9 10 PROCEDURES: … … 17 18 charpoly(A,v); characteristic polynomial of A ( using busadj(A) ) 18 19 adjoint(A); adjoint of A ( using busadj(A) ) 19 det_B(A); determinant of A . ( usesbusadj(A) )20 det_B(A); determinant of A ( using busadj(A) ) 20 21 gaussred(A); gaussian reduction: P*A=U*S, S a row reduced form of A 21 22 gaussred_pivot(A); gaussian reduction: P*A=U*S, uses row pivoting 22 23 gauss_nf(A); gaussian normal form of A 23 24 mat_rk(A); rank of constant matrix A 24 U_D_O(A); P*A=U*D*O, P,D,U,O =permutaion,diag,lower-,upper-triang25 U_D_O(A); P*A=U*D*O, P,D,U,O=permutaion,diag,lower-,upper-triang 25 26 pos_def(A,i); test symmetric matrix for positive definiteness 27 jordan(M[,opt]); eigenvalues, Jordan block sizes, transformation matrix 28 jordanmatrix(l); Jordan matrix with eigenvalues, Jordan block sizes 29 jordanform(M); Jordan normal form of constant square matrix M 26 30 "; 27 31 28 32 LIB "matrix.lib"; 29 33 LIB "ring.lib"; 30 34 LIB "elim.lib"; 31 35 ////////////////////////////////////////////////////////////////////////////// 32 36 // help functions … … 39 43 } 40 44 41 static 42 proc const_mat(matrix A) 45 static proc const_mat(matrix A) 43 46 "RETURN: 1 (0) if A is (is not) a constant matrix" 44 47 { … … 99 102 ///////////////////////////////////////////////////////////////////////////// 100 103 101 proc inverse(matrix A )104 proc inverse(matrix A, list #) 102 105 "USAGE: inverse(A); A = constant, square matrix 103 106 RETURN: matrix, the inverse matrix of A 104 NOTE: parameters and minpoly are allowed , uses std applied to A enlarged105 bei unitmatrix107 NOTE: parameters and minpoly are allowed. The procedure applies std with 108 ordering (c,dp) to A-enlarged-by-unit-matrix 106 109 EXAMPLE: example inverse; shows an example" 107 110 { … … 113 116 114 117 if (ncols(A)!=n){ 115 " input is not a square matrix";118 "// ** input is not a square matrix"; 116 119 matrix invmat; 117 120 return(invmat); 118 121 } 119 122 if(!const_mat(A)){ 120 " input is not a constant matrix";123 "// ** input is not a constant matrix"; 121 124 matrix invmat; 122 125 return(invmat); … … 130 133 matrix B=transpose(concat(A,unitmat(n))); 131 134 matrix D=transpose(std(B)); 135 136 //option(intStrategy); 137 //option(nointStrategy); 138 //if (#[1]==1) {matrix D=transpose(interred(B));} 139 132 140 matrix D1=submat(D,1..n,1..n); 133 141 matrix D2=submat(D,1..n,(n+1)..2*n); … … 146 154 for (i=1;i<=n;i=i+1){ 147 155 if (deg(B[n+i,i])<0){ 148 " matrix is not invertible";156 "// ** matrix is not invertible"; 149 157 setring BR; 150 158 kill @R;kill @@R; … … 181 189 182 190 if (ncols(A)!=n){ 183 " input is not a square matrix";;191 "// ** input is not a square matrix";; 184 192 return(A); 185 193 } 186 194 187 195 if(!const_mat(A)){ 188 " input is not a constant matrix";196 "// ** input is not a constant matrix"; 189 197 return(A); 190 198 } 191 199 192 200 if(deg(std(A-transpose(A))[1])!=-1){ 193 " input is not a symmetric matrix";201 "// ** input is not a symmetric matrix"; 194 202 return(A); 195 203 } … … 222 230 223 231 if(!const_mat(A)){ 224 " input is not a constant matrix";232 "// ** input is not a constant matrix"; 225 233 matrix B; 226 234 return(B); … … 232 240 for(j=1;j<i;j=j+1) { 233 241 k=inner_product(B[j],B[j]); 234 if (k==0) { "Error: vector with zero length"; return(matrix(B)); }242 if (k==0) { "Error: vector of length zero"; return(matrix(B)); } 235 243 B[i]=B[i]-(inner_product(B[i],B[j])/k)*B[j]; 236 244 } … … 250 258 proc diag_test(matrix A) 251 259 "USAGE: diag_test(A); A = const square matrix 252 NOTE: Der Test ist nur fuer zerfallende Matrizen aussagefaehig. 253 Parameterringe werden nicht unterstuetzt (benutzt factorize,gcd). 254 RETURN: int, 1 falls A diagonalisierbar, 0 nicht diagonalisierbar 255 -1 keine Aussage moeglich, da A nicht zerfallend. 260 RETURN: int, 1 if A is diagonalisable, 0 if not 261 -1 no statement possible, since A does not split. 262 NOTE: The test works only for split matrices, i.e if eigenvalues of A 263 are in the ground field. 264 Does not work with parameters (uses factorize,gcd). 256 265 EXAMPLE: example diag_test; shows an example" 257 266 { … … 305 314 for(i=1;i<=nf;i++){ 306 315 if(lead(Xfactors[1,i])>=@t^2){ 307 //" Die Matrix ist nicht zerfallend";316 //" matrix does not split"; 308 317 setring BR; 309 318 return(-1); … … 322 331 323 332 if(0!=reduce(b,g)){ 324 //" Die Matrix ist nicht diagonalisierbar!";333 //" matrix not diagonalizable"; 325 334 setring BR; 326 335 return(0); … … 346 355 ////////////////////////////////////////////////////////////////////////////// 347 356 proc busadj(matrix A) 348 "USAGE: busadj(A); A = square matrix 349 RETURN: list L; 350 L[1] contains the (n+1) coefficients of the characteristic polynomial 351 X of A, i.e. 357 "USAGE: busadj(A); A = square matrix (of size nxn) 358 RETURN: list L: 359 @format 360 L[1] contains the (n+1) coefficients of the characteristic 361 polynomial X of A, i.e. 352 362 X = L[1][1]+..+L[1][k]*t^(k-1)+..+(L[1][n+1])*t^n 353 363 L[2] contains the n (nxn)-matrices Hk which are the coefficients of 354 the busadjoint bA =adjoint(E*t-A) of A, i.e.364 the busadjoint bA = adjoint(E*t-A) of A, i.e. 355 365 bA = (Hn-1)*t^(n-1)+...+Hk*t^k+...+H0, ( Hk=L[2][k+1] ) 366 @end format 356 367 EXAMPLE: example busadj; shows an example" 357 368 { … … 400 411 ////////////////////////////////////////////////////////////////////////////// 401 412 proc charpoly(matrix A, list #) 402 "USAGE: charpoly(A,[v]); A = square matrix, v = name of a ringvariable 413 "USAGE: charpoly(A[,v]); A square matrix, v string, name of a variable 414 RETURN: poly, the characteristic polynomial det(E*v-A) 415 (default: v=name of last variable) 416 NOTE: A must be independent of the variable v. The computation uses det. 417 If printlevel>0, det(E*v-A) is diplayed recursively. 418 EXAMPLE: example charpoly; shows an example" 419 { 420 int n = nrows(A); 421 int z = nvars(basering); 422 int i,j; 423 string v; 424 poly X; 425 if(ncols(A) != n) 426 { 427 "// input is not a square matrix"; 428 return(X); 429 } 430 //---------------------- test for correct variable ------------------------- 431 if( size(#)==0 ){ 432 #[1] = varstr(z); 433 } 434 if( typeof(#[1]) == "string") { v = #[1]; } 435 else 436 { 437 "// 2nd argument must be a name of a variable not contained in the matrix"; 438 return(X); 439 } 440 j=-1; 441 for(i=1; i<=z; i++) 442 { 443 if(varstr(i)==v){j=i;} 444 } 445 if(j==-1) 446 { 447 "// "+v+" is not a variable in the basering"; 448 return(X); 449 } 450 if ( size(select1(module(A),j)) != 0 ) 451 { 452 "// matrix must not contain the variable "+v; 453 "// change to a ring with an extra variable, if necessary." 454 return(X); 455 } 456 //-------------------------- compute charpoly ------------------------------ 457 X = det(var(j)*unitmat(n)-A); 458 if( printlevel-voice+2 >0) { showrecursive(X,var(j));} 459 return(X); 460 } 461 example 462 { "EXAMPLE"; echo=2; 463 ring r=0,(x,t),dp; 464 matrix A[3][3]=1,x2,x,x2,6,4,x,4,1; 465 print(A); 466 charpoly(A,"t"); 467 } 468 469 ////////////////////////////////////////////////////////////////////////////// 470 proc charpoly_B(matrix A, list #) 471 "USAGE: charpoly_B(A[,v]); A square matrix, v string, name of a variable 472 RETURN: poly, the characteristic polynomial det(E*v-A) 473 (default: v=name of last variable) 403 474 NOTE: A must be constant in the variable v. The computation uses busadj(A). 404 RETURN: poly, the characteristic polynomial det(E*v-A) 405 EXAMPLE: example charpoly; shows an example" 475 EXAMPLE: example charpoly_B; shows an example" 406 476 { 407 477 int i,j; … … 420 490 //test for correct variable 421 491 if( size(#)==0 ){ 422 #[1] = varstr( 1);492 #[1] = varstr(nvars(BR)); 423 493 } 424 494 if( typeof(#[1]) == "string"){ … … 453 523 for(i=1; i<=n; i++){ 454 524 if(deg(lead(A)[i])>=1){ 455 " matrix must not contain the variable "+v;525 "// matrix must not contain the variable "+v; 456 526 kill @R; 457 527 setring BR; … … 475 545 matrix A[3][3]=1,x2,x,x2,6,4,x,4,1; 476 546 print(A); 477 charpoly (A,"t");547 charpoly_B(A,"t"); 478 548 } 479 549 … … 481 551 proc adjoint(matrix A) 482 552 "USAGE: adjoint(A); A = square matrix 553 RETURN: adjoint matrix of A, i.e. Adj*A=det(A)*E 483 554 NOTE: computation uses busadj(A) 484 RETURN: adjoint485 555 EXAMPLE: example adjoint; shows an example" 486 556 { … … 511 581 ////////////////////////////////////////////////////////////////////////////// 512 582 proc inverse_B(matrix A) 513 "USAGE: inverse_B(A); A = square matrix 514 RETURN: list Inv with Inv[1]=matrix I and Inv[2]=poly p 515 and I*A = unitmat(n)*p; 516 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not 517 EXAMPLE: example inverse_B; shows an example" 583 "USAGE: inverse_B(A); A = square matrix 584 RETURN: list Inv with 585 - Inv[1] = matrix I and 586 - Inv[2] = poly p 587 such that I*A = unitmat(n)*p; 588 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not; 589 the computation uses busadj. 590 EXAMPLE: example inverse_B; shows an example" 518 591 { 519 592 int i; … … 545 618 { "EXAMPLE"; echo=2; 546 619 ring r=0,(x,y),lp; 547 matrix A[3][3]=x,y,1,1,x2,y,x +y,6,7;620 matrix A[3][3]=x,y,1,1,x2,y,x,6,0; 548 621 print(A); 549 622 list Inv=inverse_B(A); … … 580 653 "USAGE: inverse_L(A); A = square matrix 581 654 RETURN: list Inv representing a left inverse of A, i.e 582 Inv[1] = matrix I and Inv[2]=poly p such that 583 I*A = unitmat(n)*p; 584 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not 655 - Inv[1] = matrix I and 656 - Inv[2] = poly p 657 such tha I*A = unitmat(n)*p; 658 NOTE: p=1 if 1/det(A) is computable and p=det(A) if not; 659 the computation uses lift 585 660 EXAMPLE: example inverse_L; shows an example" 586 661 { … … 621 696 { "EXAMPLE"; echo=2; 622 697 ring r=0,(x,y),lp; 623 matrix A[3][3]=x,y,1,1,x2,y,x +y,6,7;698 matrix A[3][3]=x,y,1,1,x2,y,x,6,0; 624 699 print(A); 625 700 list Inv=inverse_L(A); … … 631 706 ////////////////////////////////////////////////////////////////////////////// 632 707 proc gaussred(matrix A) 633 "USAGE: gaussred(A); A any constant matrix 634 RETURN: list Z: Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A) 635 gives a row reduced matrix S, a permutation matrix P and a 636 normalized lower triangular matrix U, with P*A=U*S 637 EXAMPLE: example gaussred; shows an example" 708 "USAGE: gaussred(A); A any constant matrix 709 RETURN: list Z: Z[1]=P , Z[2]=U , Z[3]=S , Z[4]=rank(A) 710 gives a row reduced matrix S, a permutation matrix P and a 711 normalized lower triangular matrix U, with P*A=U*S 712 NOTE: This procedure is designed for teaching purposes mainly. 713 The straight forward implementation in the interpreted library 714 is not very efficient (no standard basis computation). 715 EXAMPLE: example gaussred; shows an example" 638 716 { 639 717 int i,j,l,k,jp,rang; … … 812 890 proc gauss_nf(matrix A) 813 891 "USAGE: gauss_nf(A); A any constant matrix 814 RETURN: matrix; gauss normal form of A 892 RETURN: matrix; gauss normal form of A (uses gaussred) 815 893 EXAMPLE: example gauss_nf; shows an example" 816 894 { … … 838 916 list Z; 839 917 if(!const_mat(A)){ 840 "// input is not a constant matrix";918 "// input is not a constant matrix"; 841 919 return(-1); 842 920 } … … 857 935 gives a permutation matrix P, 858 936 a normalized lower triangular matrix U , 859 a diagonal matrix D, and a normalized upper triangular matrix O 937 a diagonal matrix D, and 938 a normalized upper triangular matrix O 860 939 with P*A=U*D*O 861 NOTE: Z[1]=-1 means that A is not regular 940 NOTE: Z[1]=-1 means that A is not regular (proc uses gaussred) 862 941 EXAMPLE: example U_D_O; shows an example" 863 942 { … … 920 999 proc pos_def(matrix A) 921 1000 "USAGE: pos_def(A); A = constant, symmetric square matrix 922 RETURN: int; 1 if A is positive definit , 0 if not, -1 if unknown 1001 RETURN: int: 1002 1 if A is positive definit , 1003 0 if not, 1004 -1 if unknown 923 1005 EXAMPLE: example pos_def; shows an example" 924 1006 { … … 981 1063 proc linsolve(matrix A, matrix b) 982 1064 "USAGE: linsolve(A,b); A a constant nxm-matrix, b a constant nx1-matrix 983 RETURN: a solution of inhomogeneous linear system A*X = b 1065 RETURN: a 1xm matrix X, solution of inhomogeneous linear system A*X = b 1066 return the 0-matrix if system is not solvable 1067 NOTE: uses gaussred 984 1068 EXAMPLE: example linsolve; shows an example" 985 1069 { … … 1089 1173 ////////////////////////////////////////////////////////////////////////////// 1090 1174 1175 /////////////////////////////////////////////////////////////////////////////// 1176 // PROCEDURES for Jordan Normal form 1177 // AUTHOR: Mathias Schulze, email: mschulze@mathematik.uni-kl.de 1178 /////////////////////////////////////////////////////////////////////////////// 1179 1180 static proc countblocks(matrix M) 1181 { 1182 int b,r,r0; 1183 1184 int i=1; 1185 while(i<=nrows(M)) 1186 { 1187 b++; 1188 r=nrows(M[i]); 1189 r0=r; 1190 1191 dbprint(printlevel-voice+2,"//searching for block "+string(b)+"..."); 1192 while(i<r0&&i<nrows(M)) 1193 { 1194 i++; 1195 if(i<=nrows(M)) 1196 { 1197 r=nrows(M[i]); 1198 if(r>r0) 1199 { 1200 r0=r; 1201 } 1202 } 1203 } 1204 dbprint(printlevel-voice+2,"//...block "+string(b)+" found"); 1205 1206 i++; 1207 } 1208 1209 return(b); 1210 } 1211 /////////////////////////////////////////////////////////////////////////////// 1212 1213 static proc getblock(matrix M,intvec v) 1214 { 1215 matrix M0[size(v)][size(v)]=M[v,v]; 1216 return(M0); 1217 } 1218 /////////////////////////////////////////////////////////////////////////////// 1219 1220 proc jordan(matrix M,list #) 1221 "USAGE: jordan(M[,opt]); M constant square matrix, opt integer 1222 ASSUME: The eigenvalues of M are in the coefficient field. 1223 RETURN: The procedure returns a list jd with 3 entries: 1224 @format 1225 - jd[1] ideal, eigenvalues of M, 1226 - jd[2] list of intvecs, jd[2][i][j] size of j-th Jordan block with 1227 eigenvalue jd[1][i], and 1228 - jd[3] a matrix, jd[3]^(-1)*M*jd[3] in Jordan normal form. 1229 Depending on opt, only certain entries of jd are computed. 1230 If opt=-1, jd[1] is computed, 1231 if opt= 0, jd[1] and jd[2] are computed, 1232 if opt= 1, jd[1], jd[2], and jd[3] are computed, and, 1233 if opt= 2, jd[1] and jd[3] are computed. 1234 By default, opt=0. 1235 @end format 1236 NOTE: A non constant polynomial matrix M is replaced by its constant part. 1237 DISPLAY: The procedure displays comments if printlevel>=1. 1238 EXAMPLE: example jordan; shows an example. 1239 " 1240 { 1241 int n=nrows(M); 1242 if(n==0) 1243 { 1244 print("//empty matrix"); 1245 return(list()); 1246 } 1247 if(n!=ncols(M)) 1248 { 1249 print("//no square matrix"); 1250 return(list()); 1251 } 1252 1253 M=jet(M,0); 1254 1255 dbprint(printlevel-voice+2,"//counting blocks of matrix..."); 1256 int i=countblocks(M); 1257 dbprint(printlevel-voice+2,"//...blocks of matrix counted"); 1258 if(i==1) 1259 { 1260 dbprint(printlevel-voice+2,"//matrix has 1 block"); 1261 } 1262 else 1263 { 1264 dbprint(printlevel-voice+2,"//matrix has "+string(i)+" blocks"); 1265 } 1266 1267 dbprint(printlevel-voice+2,"//counting blocks of transposed matrix..."); 1268 int j=countblocks(transpose(M)); 1269 dbprint(printlevel-voice+2,"//...blocks of transposed matrix counted"); 1270 if(j==1) 1271 { 1272 dbprint(printlevel-voice+2,"//transposed matrix has 1 block"); 1273 } 1274 else 1275 { 1276 dbprint(printlevel-voice+2,"//transposed matrix has "+string(j)+" blocks"); 1277 } 1278 1279 if(i<j) 1280 { 1281 dbprint(printlevel-voice+2,"//transposing matrix..."); 1282 M=transpose(M); 1283 dbprint(printlevel-voice+2,"//...matrix transposed"); 1284 } 1285 1286 list fd; 1287 matrix M0; 1288 poly cp; 1289 ideal eM,eM0; 1290 intvec mM,mM0; 1291 intvec u; 1292 int b,r,r0; 1293 1294 i=1; 1295 while(i<=nrows(M)) 1296 { 1297 b++; 1298 u=i; 1299 r=nrows(M[i]); 1300 r0=r; 1301 1302 dbprint(printlevel-voice+2,"//searching for block "+string(b)+"..."); 1303 while(i<r0&&i<nrows(M)) 1304 { 1305 i++; 1306 if(i<=nrows(M)) 1307 { 1308 u=u,i; 1309 r=nrows(M[i]); 1310 if(r>r0) 1311 { 1312 r0=r; 1313 } 1314 } 1315 } 1316 dbprint(printlevel-voice+2,"//...block "+string(b)+" found"); 1317 1318 if(size(u)==1) 1319 { 1320 dbprint(printlevel-voice+2,"//1x1-block:"); 1321 dbprint(printlevel-voice+2,M[u[1]][u[1]]); 1322 1323 if(mM[1]==0) 1324 { 1325 eM=M[u[1]][u[1]]; 1326 mM=1; 1327 } 1328 else 1329 { 1330 eM=eM,ideal(M[u[1]][u[1]]); 1331 mM=mM,1; 1332 } 1333 } 1334 else 1335 { 1336 dbprint(printlevel-voice+2, 1337 "//"+string(size(u))+"x"+string(size(u))+"-block:"); 1338 M0=getblock(M,u); 1339 dbprint(printlevel-voice+2,M0); 1340 1341 dbprint(printlevel-voice+2,"//characteristic polynomial:"); 1342 cp=det(module(M0-var(1)*freemodule(size(u)))); 1343 dbprint(printlevel-voice+2,cp); 1344 1345 dbprint(printlevel-voice+2,"//factorizing characteristic polynomial..."); 1346 fd=factorize(cp,2); 1347 dbprint(printlevel-voice+2,"//...characteristic polynomial factorized"); 1348 1349 dbprint(printlevel-voice+2,"//computing eigenvalues..."); 1350 eM0,mM0=fd[1..2]; 1351 if(1<var(1)) 1352 { 1353 for(j=ncols(eM0);j>=1;j--) 1354 { 1355 if(deg(eM0[j])>1) 1356 { 1357 print("//eigenvalues not in the coefficient field"); 1358 return(list()); 1359 } 1360 if(eM0[j][1]==0) 1361 { 1362 eM0[j]=0; 1363 } 1364 else 1365 { 1366 eM0[j]=-eM0[j][2]/(eM0[j][1]/var(1)); 1367 } 1368 } 1369 } 1370 else 1371 { 1372 for(j=ncols(eM0);j>=1;j--) 1373 { 1374 if(deg(eM0[j])>1) 1375 { 1376 print("//eigenvalues not in the coefficient field"); 1377 return(list()); 1378 } 1379 if(eM0[j][2]==0) 1380 { 1381 eM0[j]=0; 1382 } 1383 else 1384 { 1385 eM0[j]=-eM0[j][1]/(eM0[j][2]/var(1)); 1386 } 1387 } 1388 } 1389 dbprint(printlevel-voice+2,"//...eigenvalues computed"); 1390 1391 if(mM[1]==0) 1392 { 1393 eM=eM0; 1394 mM=mM0; 1395 } 1396 else 1397 { 1398 eM=eM,eM0; 1399 mM=mM,mM0; 1400 } 1401 } 1402 1403 i++; 1404 } 1405 1406 dbprint(printlevel-voice+2,"//sorting eigenvalues..."); 1407 poly e; 1408 int m; 1409 for(i=ncols(eM);i>=2;i--) 1410 { 1411 for(j=i-1;j>=1;j--) 1412 { 1413 if(eM[i]<eM[j]) 1414 { 1415 e=eM[i]; 1416 eM[i]=eM[j]; 1417 eM[j]=e; 1418 m=mM[i]; 1419 mM[i]=mM[j]; 1420 mM[j]=m; 1421 } 1422 } 1423 } 1424 dbprint(printlevel-voice+2,"//...eigenvalues sorted"); 1425 1426 dbprint(printlevel-voice+2,"//removing multiple eigenvalues..."); 1427 i=1; 1428 j=2; 1429 while(j<=ncols(eM)) 1430 { 1431 if(eM[i]==eM[j]) 1432 { 1433 mM[i]=mM[i]+mM[j]; 1434 } 1435 else 1436 { 1437 i++; 1438 eM[i]=eM[j]; 1439 mM[i]=mM[j]; 1440 } 1441 j++; 1442 } 1443 eM=eM[1..i]; 1444 mM=mM[1..i]; 1445 dbprint(printlevel-voice+2,"//...multiple eigenvalues removed"); 1446 1447 dbprint(printlevel-voice+2,"//eigenvalues:"); 1448 dbprint(printlevel-voice+2,eM); 1449 dbprint(printlevel-voice+2,"//multiplicities:"); 1450 dbprint(printlevel-voice+2,mM); 1451 1452 int opt=0; 1453 if(size(#)>0) 1454 { 1455 if(typeof(#[1])=="int") 1456 { 1457 opt=#[1]; 1458 } 1459 } 1460 if(opt<0) 1461 { 1462 return(list(eM)); 1463 } 1464 int k,l; 1465 matrix I=freemodule(n); 1466 matrix Mi,Ni; 1467 module sNi; 1468 list K; 1469 if(opt>=1) 1470 { 1471 module V,K1,K2; 1472 matrix v[n][1]; 1473 } 1474 if(opt<=1) 1475 { 1476 list bM; 1477 intvec bMi; 1478 } 1479 1480 for(i=ncols(eM);i>=1;i--) 1481 { 1482 Mi=M-eM[i]*I; 1483 1484 dbprint(printlevel-voice+2, 1485 "//computing kernels of powers of matrix minus eigenvalue " 1486 +string(eM[i])); 1487 K=list(module()); 1488 for(Ni,sNi=Mi,0;size(sNi)<mM[i];Ni=Ni*Mi) 1489 { 1490 sNi=syz(Ni); 1491 K=K+list(sNi); 1492 } 1493 dbprint(printlevel-voice+2,"//...kernels computed"); 1494 1495 if(opt<=1) 1496 { 1497 dbprint(printlevel-voice+2, 1498 "//computing Jordan block sizes for eigenvalue " 1499 +string(eM[i])+"..."); 1500 bMi=0; 1501 bMi[size(K[2])]=0; 1502 for(j=size(K);j>=2;j--) 1503 { 1504 for(k=size(bMi);k>size(bMi)+size(K[j-1])-size(K[j]);k--) 1505 { 1506 bMi[k]=bMi[k]+1; 1507 } 1508 } 1509 bM=list(bMi)+bM; 1510 dbprint(printlevel-voice+2,"//...Jordan block sizes computed"); 1511 } 1512 1513 if(opt>=1) 1514 { 1515 dbprint(printlevel-voice+2, 1516 "//computing Jordan basis vectors for eigenvalue " 1517 +string(eM[i])+"..."); 1518 if(size(K)>1) 1519 { 1520 for(j,K1=2,0;j<=size(K)-1;j++) 1521 { 1522 K2=K[j]; 1523 K[j]=interred(reduce(K[j],std(K1+module(Mi*K[j+1])))); 1524 K1=K2; 1525 } 1526 K[j]=interred(reduce(K[j],std(K1))); 1527 } 1528 for(j=size(K);j>=2;j--) 1529 { 1530 for(k=size(K[j]);k>=1;k--) 1531 { 1532 v=K[j][k]; 1533 for(l=j;l>=1;l--) 1534 { 1535 V=module(v)+V; 1536 v=Mi*v; 1537 } 1538 } 1539 } 1540 dbprint(printlevel-voice+2,"//...Jordan basis vectors computed"); 1541 } 1542 } 1543 1544 list jd=eM; 1545 if(opt<=1) 1546 { 1547 jd[2]=bM; 1548 } 1549 if(opt>=1) 1550 { 1551 jd[3]=V; 1552 } 1553 return(jd); 1554 } 1555 example 1556 { "EXAMPLE:"; echo=2; 1557 ring R=0,x,dp; 1558 matrix M[3][3]=3,2,1,0,2,1,0,0,3; 1559 print(M); 1560 jordan(M); 1561 } 1562 /////////////////////////////////////////////////////////////////////////////// 1563 1564 proc jordanmatrix(list jd) 1565 "USAGE: jordanmatrix(jd); jd list of ideal and list of intvecs 1566 RETURN: The procedure returns the Jordan matrix J, satisfying: 1567 @* - eigenvalues of J are given by jd[1] 1568 @* - j-th Jordan block of J with eigenvalue jd[1][i] has size jd[2][i][j] 1569 DISPLAY: The procedure displays comments if printlevel>=1. 1570 EXAMPLE: example jordanmatrix; shows an example. 1571 " 1572 { 1573 if(size(jd)<2) 1574 { 1575 print("//not enough entries in argument list"); 1576 matrix J[1][0]; 1577 return(J); 1578 } 1579 def eJ,bJ=jd[1..2]; 1580 if(typeof(eJ)!="ideal") 1581 { 1582 print("//first entry in argument list not an ideal"); 1583 matrix J[1][0]; 1584 return(J); 1585 } 1586 if(typeof(bJ)!="list") 1587 { 1588 print("//second entry in argument list not a list"); 1589 matrix J[1][0]; 1590 return(J); 1591 } 1592 if(size(eJ)<size(bJ)) 1593 { 1594 int s=size(eJ); 1595 } 1596 else 1597 { 1598 int s=size(bJ); 1599 } 1600 1601 int i,j,k,n; 1602 for(i=s;i>=1;i--) 1603 { 1604 if(typeof(bJ[i])!="intvec") 1605 { 1606 print("//second entry in argument list not a list of intvecs"); 1607 matrix J[1][0]; 1608 return(J); 1609 } 1610 else 1611 { 1612 for(j=size(bJ[i]);j>=1;j--) 1613 { 1614 k=bJ[i][j]; 1615 if(k>0) 1616 { 1617 n=n+k; 1618 } 1619 } 1620 } 1621 } 1622 1623 int l; 1624 matrix J[n][n]; 1625 for(i,l=1,1;i<=s;i++) 1626 { 1627 for(j=1;j<=size(bJ[i]);j++) 1628 { 1629 k=bJ[i][j]; 1630 if(k>0) 1631 { 1632 while(k>=2) 1633 { 1634 J[l,l]=eJ[i]; 1635 J[l,l+1]=1; 1636 k,l=k-1,l+1; 1637 } 1638 J[l,l]=eJ[i]; 1639 l++; 1640 } 1641 } 1642 } 1643 1644 return(J); 1645 } 1646 example 1647 { "EXAMPLE:"; echo=2; 1648 ring R=0,x,dp; 1649 list l; 1650 l[1]=ideal(2,3); 1651 l[2]=list(intvec(1),intvec(2)); 1652 print(jordanmatrix(l)); 1653 } 1654 /////////////////////////////////////////////////////////////////////////////// 1655 1656 proc jordanform(matrix M) 1657 "USAGE: jordanform(M); M constant square matrix 1658 ASSUME: The eigenvalues of M are in the coefficient field. 1659 RETURN: matrix, the Jordan normal form of M. 1660 NOTE: A non constant polynomial matrix M is replaced by its constant part. 1661 DISPLAY: The procedure displays more comments for higher printlevel. 1662 EXAMPLE: example jordanform; shows an example. 1663 " 1664 { 1665 return(jordanmatrix(jordan(M))); 1666 } 1667 example 1668 { "EXAMPLE:"; echo=2; 1669 ring R=0,x,dp; 1670 matrix M[3][3]=3,2,1,0,2,1,0,0,3; 1671 print(M); 1672 print(jordanform(M)); 1673 } 1674 ///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset
for help on using the changeset viewer.