share/linearalgebra/rtest_linalg.mac

   1 (matrix_element_add : "+", matrix_element_mult : "*",0);
   2 0$
   3
   4 /*--- ctranspose -----*/
   5
   6 (m : matrix([%i,1],[2,%i],[0,2]),0);
   7 0$
   8
   9 ctranspose(m);
  10 matrix([-%i,2,0],[1,-%i,2])$
  11
  12 zeromatrixp(m - ctranspose(ctranspose(m)));
  13 true$
  14
  15 (m : matrix([m,m]),0);
  16 0$
  17
  18 (matrix_element_transpose : 'transpose,0);
  19 0$
  20
  21 zeromatrixp(ctranspose(mat_unblocker(m)) - mat_unblocker(ctranspose(m)));
  22 true$
  23
  24 ctranspose(matrix([%i]));
  25 matrix([-%i])$
  26
  27 (matrix_element_transpose : false,0);
  28 0$
  29
  30 /*--- matrix norm ----*/
  31
  32 (m : matrix([2,3],[6,7]), 0);
  33 0$
  34
  35 mat_norm(m,1);
  36 10$
  37
  38 mat_norm(m,'inf);
  39 13$
  40
  41 mat_norm(m,'frobenius);
  42 sqrt(2^2 + 3^2 + 6^2 + 7^2)$
  43
  44 (m : matrix([2*%i,3],[6,7]), 0);
  45 0$
  46
  47 mat_norm(m,1);
  48 10$
  49
  50 mat_norm(m,'inf);
  51 13$
  52
  53 mat_norm(m,'frobenius);
  54 sqrt(2^2 + 3^2 + 6^2 + 7^2)$
  55
  56
  57 /*--- matrix size ----*/
  58
  59 matrix_size(matrix([1]));
  60 [1,1]$
  61
  62 matrix_size(matrix([false],[false]));
  63 [2,1]$
  64
  65 matrix_size(matrix([false,false]));
  66 [1,2]$
  67
  68 /*--- matrix linalg_rank and nullity -----*/
  69
  70 linalg_rank(matrix([0]));
  71 0$
  72
  73 nullity(matrix([0]));
  74 1$
  75
  76 linalg_rank(matrix([1]));
  77 1$
  78
  79 nullity(matrix([1]));
  80 0$
  81
  82 linalg_rank(matrix([0,0]));
  83 0$
  84
  85 nullity(matrix([0,0]));
  86 2$
  87
  88 linalg_rank(matrix([0,1]));
  89 1$
  90
  91 nullity(matrix([0,1]));
  92 1$
  93
  94 linalg_rank(matrix([0],[0]));
  95 0$
  96
  97 nullity(matrix([0],[0]));
  98 1$
  99
 100 linalg_rank(matrix([1],[0]));
 101 1$
 102
 103 nullity(matrix([1],[0]));
 104 0$
 105
 106 linalg_rank(matrix([0],[1]));
 107 1$
 108
 109 nullity(matrix([0],[1]));
 110 0$
 111
 112 linalg_rank(matrix([1],[1]));
 113 1$
 114
 115 nullity(matrix([1],[1]));
 116 0$
 117
 118 linalg_rank(matrix([0],[%i]));
 119 1$
 120
 121 nullity(matrix([0],[%i]));
 122 0$
 123
 124 linalg_rank(matrix([1],[1]));
 125 1$
 126
 127 nullity(matrix([1],[1]));
 128 0$
 129
 130 linalg_rank(matrix([1],[%i]));
 131 1$
 132
 133 nullity(matrix([1],[%i]));
 134 0$
 135
 136 linalg_rank(matrix([1,2],[0,0]));
 137 1$
 138
 139 nullity(matrix([1,2],[0,0]));
 140 1$
 141
 142 linalg_rank(matrix([1,2],[0,1]));
 143 2$
 144
 145 nullity(matrix([1,2],[0,1]));
 146 0$
 147
 148
 149 /*--- dotproducts ------*/
 150
 151 dotproduct(matrix([1]),matrix([1]));
 152 1$
 153
 154 dotproduct(matrix([%i]),matrix([%i]));
 155 1$
 156
 157 dotproduct(matrix([1],[2]),matrix([-2],[1]));
 158 0$
 159
 160 dotproduct(matrix([1],[%i]),matrix([1],[%i]));
 161 2$
 162
 163 /*--- matrix nullity ----*/
 164
 165 nullity(matrix([0]));
 166 1$
 167
 168 nullity(matrix([1]));
 169 0$
 170
 171 nullity(matrix([0,0]));
 172 2$
 173
 174 nullity(matrix([0,1]));
 175 1$
 176
 177 nullity(matrix([1,0]));
 178 1$
 179
 180 nullity(matrix([1,2],[0,0]));
 181 1$
 182
 183 nullity(matrix([1,2],[0,1]));
 184 0$
 185
 186 /*---orthogonal complements ---*/
 187
 188 (m : matrix([1,2,3]),0);
 189 0$
 190
 191 (s1 : args(columnspace(m)),0);
 192 0$
 193
 194 (s2 : args(nullspace(ctranspose(m))),0);
 195 0$
 196
 197 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 198      cartesian_product(setify(s1),setify(s2))));
 199 true$
 200
 201 /*----*/
 202
 203 (m : matrix([0,0,3]),0);
 204 0$
 205
 206 (s1 : args(columnspace(m)),0);
 207 0$
 208
 209 (s2 : args(nullspace(ctranspose(m))),0);
 210 0$
 211
 212 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 213      cartesian_product(setify(s1),setify(s2))));
 214
 215 true$
 216
 217 /*----*/
 218
 219 (m : matrix([0,0,0]),0);
 220 0$
 221
 222 (s1 : args(columnspace(m)),0);
 223 0$
 224
 225 (s2 : args(nullspace(ctranspose(m))),0);
 226 0$
 227
 228 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 229      cartesian_product(setify(s1),setify(s2))));
 230
 231 true$
 232
 233 /*----*/
 234
 235 (m : matrix([0,0,3],[0,0,1]),0);
 236 0$
 237
 238 (s1 : args(columnspace(m)),0);
 239 0$
 240
 241 (s2 : args(nullspace(ctranspose(m))),0);
 242 0$
 243
 244 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 245      cartesian_product(setify(s1),setify(s2))));
 246
 247 true$
 248
 249 /*----*/
 250
 251 (m : matrix([1,2,3],[%i,2,1]),0);
 252 0$
 253
 254 (s1 : args(columnspace(m)),0);
 255 0$
 256
 257 (s2 : args(nullspace(ctranspose(m))),0);
 258 0$
 259
 260 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 261      cartesian_product(setify(s1),setify(s2))));
 262
 263 true$
 264
 265
 266 /*----*/
 267
 268 (m : matrix([1,2,3],[%i,2,1],[1,2,3],[1,2,3]),0);
 269 0$
 270
 271 (s1 : args(columnspace(m)),0);
 272 0$
 273
 274 (s2 : args(nullspace(ctranspose(m))),0);
 275 0$
 276
 277 every(lambda([e], e=0), map(lambda([e], expand(apply('dotproduct,e))),
 278      cartesian_product(setify(s1),setify(s2))));
 279
 280 true$
 281
 282 /*----*/
 283
 284 (m : matrix([1,2,3],[4,5,6],[7,8,9]),0);
 285 0$
 286
 287 (s1 : args(columnspace(m)),0);
 288 0$
 289
 290 (s2 : args(nullspace(ctranspose(m))),0);
 291 0$
 292
 293 every(lambda([e], e=0), map(lambda([e], apply('dotproduct,e)),
 294      cartesian_product(setify(s1),setify(s2))));
 295
 296 true$
 297
 298
 299 /*-----*/
 300
 301 ptriangularize(matrix([1]),z);
 302 matrix([1])$
 303
 304 ptriangularize(matrix([0]),z);
 305 matrix([0])$
 306
 307 ptriangularize(matrix([z]),z);
 308 matrix([z])$
 309
 310 (m : matrix([1,2],[3,4]) - z*ident(2),0);
 311 0$
 312
 313 abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z))));
 314 1$
 315
 316 (m : matrix([1,2,3],[4,5,6], [7,8,9]) - z*ident(3),0);
 317 0$
 318
 319 abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z))));
 320 1$
 321
 322 (m : matrix([1-z,2,3-z],[4,5-z,6], [7-z,8,9+z]),0);
 323 0$
 324
 325 abs(ratsimp(determinant(m) / determinant(ptriangularize(m,z))));
 326 1$
 327
 328 locate_matrix_entry(matrix([1]),1,1,1,1,lambda([e],e # 0), 'bool);
 329 [1,1]$
 330
 331 locate_matrix_entry(matrix([0]),1,1,1,1,lambda([e],e # 0), 'bool);
 332 false$
 333
 334 locate_matrix_entry(matrix([1]),1,1,1,1,identity, 'min);
 335 [1,1]$
 336
 337 locate_matrix_entry(matrix([1]),1,1,1,1,identity, 'max);
 338 [1,1]$
 339
 340 locate_matrix_entry(matrix([1,1]),1,1,1,2,identity, 'min);
 341 [1,1]$
 342
 343 locate_matrix_entry(matrix([1,1]),1,1,1,2,identity, 'max);
 344 [1,1]$
 345
 346 locate_matrix_entry(matrix([1,2]),1,1,1,2,identity, 'max);
 347 [1,2]$
 348
 349 locate_matrix_entry(matrix([1,2]),1,1,1,2,identity, 'min);
 350 [1,1]$
 351
 352 (m : matrix([1,2,3],[4,5,6],[7,8,9]),0);
 353 0$
 354
 355 locate_matrix_entry(m,1,1,3,3,lambda([e],e+%pi), min);
 356 [1,1]$
 357
 358 locate_matrix_entry(m,1,1,3,3,lambda([e],e+%pi), max);
 359 [3,3]$
 360
 361 locate_matrix_entry(m,1,2,3,3,lambda([e],e+%pi), min);
 362 [1,2]$
 363
 364 locate_matrix_entry(m,1,2,3,3,lambda([e],e+%pi), max);
 365 [3,3]$
 366
 367 columnop(matrix([2,1]),1,2,-2);
 368 matrix([4,1])$
 369
 370 columnop(matrix([2,1]),2,1,x);
 371 matrix([2,1-2*x])$
 372
 373 columnswap(matrix([a,b]),1,2);
 374 matrix([b,a])$
 375
 376 columnswap(matrix([a,b,c]),2,3);
 377 matrix([a,c,b])$
 378
 379 columnswap(columnswap(matrix([a,b],[c,d]),1,2),1,2);
 380 matrix([a,b],[c,d])$
 381
 382 ctranspose(matrix([5+%i*6]));
 383 matrix([5-%i*6])$
 384
 385 ctranspose(matrix([1],[%i]));
 386 matrix([1,-%i])$
 387
 388 (load(unwind_protect), unwind_protect((p[1] : 2014, rowop(matrix([a],[b]),1,2,x)), remarray(p)));
 389 matrix([a-b*x],[b])$
 390
 391 rowop(matrix([a],[b]),1,2,x);
 392 matrix([a-b*x],[b]);
 393
 394 rowswap(rowswap(matrix([a,b],[c,d]),1,2),1,2);
 395 matrix([a,b],[c,d])$
 396
 397 (m : matrix([0,1,1],[0,0,0],[0,0,0]),0);
 398 0$
 399
 400 nullspace(m);
 401 ''(funmake('span, sort([matrix([1],[0],[0]),matrix([0],[-1],[1])])))$
 402
 403 columnspace(m);
 404 span(matrix([1],[0],[0]))$
 405
 406 zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m))));
 407 true$
 408
 409 (m : matrix([0,1,2,3],[0,0,1,0],[0,0,0,0]),0);
 410 0$
 411
 412 zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m))));
 413 true$
 414
 415 (m : matrix([0,0,2,3],[0,0,1,0],[0,0,0,0]),0);
 416 0$
 417
 418 zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m))));
 419 true$
 420
 421 (m : matrix([0,0,%i,3],[0,0,%i,0],[0,0,0,0]),0);
 422 0$
 423
 424 zeromatrixp(ctranspose(apply(addcol,args(columnspace(ctranspose(m))))) . apply(addcol,args(nullspace(m))));
 425 true$
 426
 427 (remvalue(s1,s2,m),0);
 428 0$
 429
 430 hilbert_matrix(2);
 431 matrix([1,1/2],[1/2,1/3])$
 432
 433 kronecker_product(matrix([a]),matrix([b]));
 434 matrix([a*b]);
 435
 436 kronecker_product(matrix([1,2],[3,4]),matrix([a,b],[c,d]));
 437 ''(mat_unblocker(matrix([matrix([a,b],[c,d]), 2*matrix([a,b],[c,d])],[3*matrix([a,b],[c,d]),4*matrix([a,b],[c,d])])))$
 438
 439 every(lambda([e], is(e < 0)), [-1,-2,3]);
 440 false$
 441
 442 every(lambda([e], ratnump(e)), [-1,-2,3]);
 443 true$
 444
 445 every('bfloatp, matrix([1.3b0,7.8b0],[5.7b0, 9.1b0]));
 446 true$
 447
 448 every('bfloatp, matrix([matrix([1.3b0,7.8b0],[5.7b0, 7/8])]));
 449 false$
 450
 451 polytocompanion(a*x^3+b*x^2+c*x+d,x);
 452 matrix([0,0,-d/a],[1,0,-c/a],[0,1,-b/a])$
 453
 454 vandermonde_matrix([a]);
 455 matrix([1])$
 456
 457 vandermonde_matrix([a,b]);
 458 matrix([1,a],[1,b])$
 459
 460 vandermonde_matrix([]);
 461 matrix()$
 462
 463 /* SF bug #1659985 */
 464
 465 vandermonde_matrix([0]);
 466 matrix([1])$
 467
 468 vandermonde_matrix([0,1]);
 469 matrix([1,0],[1,1])$
 470
 471
 472 /*--- copylist and copymatrix ----*/
 473
 474 (m :  [1,[2,3]],0);
 475 0$
 476
 477 (mm : copylist(m),0);
 478 0$
 479
 480 (mm[2][1] : x,0);
 481 0$
 482
 483 is(m[2][1] = 2);
 484 true$
 485
 486 (m :  matrix([1,2],[3,4]),0);
 487 0$
 488
 489 (mm : copymatrix(m),0);
 490 0$
 491
 492 (mm[2,2] : x,0);
 493 0$
 494
 495 is(m[2,2] = 4);
 496 true$
 497
 498 /*-- mat_unblocker & friends --*/
 499
 500 (m : matrix([0,1],[1,0]),0);
 501 0$
 502
 503 (m1 : matrix([m,m],[m,m]),0);
 504 0$
 505
 506 mat_unblocker(m1);
 507 matrix([0,1,0,1],[1,0,1,0],[0,1,0,1],[1,0,1,0])$
 508
 509 mat_fullunblocker(m1);
 510 matrix([0,1,0,1],[1,0,1,0],[0,1,0,1],[1,0,1,0])$
 511
 512 (m2 : matrix([m1,m1],[m1,m1]),0);
 513 0$
 514
 515 mat_unblocker(m2);
 516 ''(matrix([m,m,m,m],[m,m,m,m],[m,m,m,m],[m,m,m,m]))$
 517
 518 mat_fullunblocker(m2);
 519  matrix([0,1,0,1,0,1,0,1],[1,0,1,0,1,0,1,0],[0,1,0,1,0,1,0,1],[1,0,1,0,1,0,1,0],[0,1,0,1,0,1,0,1],[1,0,1,0,1,0,1,0],[0,1,0,1,0,1,0,1],[1,0,1,0,1,0,1,0])$
 520
 521 mat_trace(m1);
 522 0$
 523
 524 mat_trace(identfor(m1));
 525 4$
 526
 527 mat_trace(m1 + identfor(m1));
 528 4$
 529
 530 mat_norm(m1,1);
 531 2$
 532
 533 mat_norm(m2,inf) - mat_norm(mat_fullunblocker(m2),inf);
 534 0$
 535
 536 blockmatrixp(m);
 537 false$
 538
 539 blockmatrixp(m1);
 540 true$
 541
 542 blockmatrixp(2/3);
 543 false$
 544
 545 blockmatrixp(rat(x));
 546 false$
 547
 548 /*---add matrices---*/
 549
 550 (m1 : matrix([8,2],[7,4]),0);
 551 0$
 552
 553 (m2 : matrix([1,2],[3,42]),0);
 554 0$
 555
 556 (m3 : matrix([96,5],[3,42]),0);
 557 0$
 558
 559 addmatrices('max,m1,m2,m3);
 560 matrix([96,5],[7,42])$
 561
 562 addmatrices("+", m1,m2);
 563 matrix([9,4],[10,46])$
 564
 565 toeplitz([1,2,3],[x,y,z]);
 566 matrix([1,y,z],[2,1,y],[3,2,1])$
 567
 568 toeplitz([1,2,3]);
 569 matrix([1,2,3],[2,1,2],[3,2,1])$
 570
 571 toeplitz([1,2,3],[x]);
 572 matrix([1],[2],[3])$
 573
 574 hankel([1,2,3]);
 575 matrix([1,2,3],[2,3,0],[3,0,0])$
 576
 577 hankel([1,2,3],[x]);
 578 matrix([1],[2],[3])$
 579
 580 diag_matrix(1);
 581 matrix([1])$
 582
 583 diag_matrix(1,2);
 584 matrix([1,0],[0,2])$
 585
 586 diag_matrix(diag_matrix(1,2),diag_matrix(3,4));
 587 matrix([matrix([1,0],[0,2]),matrix([0,0],[0,0])],[matrix([0,0],[0,0]),matrix([3,0],[0,4])])$
 588
 589 identfor(diag_matrix(diag_matrix(1,2),diag_matrix(3,4)));
 590 ''(diag_matrix(diag_matrix(1,1),diag_matrix(1,1)))$
 591
 592 zerofor(diag_matrix(diag_matrix(1,2),diag_matrix(3,4)));
 593 ''(diag_matrix(diag_matrix(0,0),diag_matrix(0,0)))$
 594
 595 (m : diag_matrix(diag_matrix(5,6),diag_matrix(8,10)),0);
 596 0$
 597
 598 mat_norm(m,inf) - mat_norm(mat_unblocker(m),inf);
 599 0$
 600
 601 mat_trace(m) - mat_trace(mat_unblocker(m));
 602 0$
 603
 604 hilbert_matrix(3);
 605 matrix([1,1/2,1/3],[1/2,1/3,1/4],[1/3,1/4,1/5])$
 606
 607 vandermonde_matrix([a,b,c]);
 608 matrix([1,a,a^2],[1,b,b^2],[1,c,c^2])$
 609
 610 every('integerp, [1,2,3]);
 611 true$
 612
 613 every('integerp, []);
 614 true$
 615
 616 every('integerp, [1,2,3,x]);
 617 false$
 618
 619 every('integerp, [false]);
 620 false$
 621
 622 every(lambda([x], is (x < 5)), [4,5,6,7]);
 623 false$
 624
 625 errcatch(every('integerp, false));
 626 []$
 627
 628 every(lambda([x], is(x < 1932)), matrix([1,2,3],[8,8,74]));
 629 true$
 630
 631 (remvalue(mm,m1,m2,m3,m),0);
 632 0$
 633
 634 columnspace(matrix([0,1],[0,1]));
 635 span(matrix([1],[1]))$
 636
 637 linalg_rank(matrix([0,1],[0,1]));
 638 1$
 639
 640 sylvester_matrix(0, 0, x);
 641 matrix()$
 642
 643 sylvester_matrix(1, 0, x);
 644 matrix()$
 645
 646 sylvester_matrix(0, 1, x);
 647 matrix()$
 648
 649 sylvester_matrix(1, 1, x);
 650 matrix()$
 651
 652 sylvester_matrix(x, 1, x);
 653 matrix([1])$
 654
 655 sylvester_matrix(1, x, x);
 656 matrix([1])$
 657
 658 sylvester_matrix(x, x, x);
 659 matrix([1,0],[1,0])$
 660
 661 sylvester_matrix(-x, x, x);
 662 matrix([-1,0],[1,0])$
 663
 664 sylvester_matrix(a+b*x, x, x);
 665 matrix([b,a],[1,0])$
 666
 667 sylvester_matrix(a+b*x, x^2, x);
 668 matrix([b,a,0],[0,b,a],[1,0,0])$
 669
 670 /* blockmatrixp -- SF bug # 1905668*/
 671
 672 blockmatrixp(matrix());
 673 false$
 674
 675 blockmatrixp(matrix([]));
 676 false$
 677
 678 blockmatrixp(matrix([],[]));
 679 false$
 680
 681 blockmatrixp(matrix([false]));
 682 false$
 683
 684 blockmatrixp(matrix([matrix()]));
 685 true$
 686
 687 blockmatrixp(matrix([matrix([false])]));
 688 true$
 689
 690 mat_cond(matrix([1,2,3],[4,5,6],[7,8,10]),1);
 691 133$
 692
 693 mat_cond(matrix([1,2,3],[4,5,6],[7,8,10]),inf);
 694 475/3$
 695
 696 /* SF bug #3158: "triangularize gives incorrect result on a matrix containing %i" */
 697
 698 nullspace (matrix([0, 1, 1, 0], [-1, 0, 0, 1], [0, 0, 0, 1], [0, 0, -1, 0]) + %i * ident(4));
 699 span(matrix([2*%i],[2],[0],[0]));
 700