Artifact 37d558ef66b31ba3ce5d41f775e16245814231295f2f7d4f7069d917d438bc1b:
- Executable file
r36/xlog/LINALG.LOG
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 29418) [annotate] [blame] [check-ins using] [more...]
REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ... if lisp !*rounded then rounded_was_on := t else rounded_was_on := nil; mat1 := mat((1,2,3,4,5),(2,3,4,5,6),(3,4,5,6,7),(4,5,6,7,8),(5,6,7,8,9)); [1 2 3 4 5] [ ] [2 3 4 5 6] [ ] mat1 := [3 4 5 6 7] [ ] [4 5 6 7 8] [ ] [5 6 7 8 9] mat2 := mat((1,1,1,1),(2,2,2,2),(3,3,3,3),(4,4,4,4)); [1 1 1 1] [ ] [2 2 2 2] mat2 := [ ] [3 3 3 3] [ ] [4 4 4 4] mat3 := mat((x),(x),(x),(x)); [x] [ ] [x] mat3 := [ ] [x] [ ] [x] mat4 := mat((3,3),(4,4),(5,5),(6,6)); [3 3] [ ] [4 4] mat4 := [ ] [5 5] [ ] [6 6] mat5 := mat((1,2,1,1),(1,2,3,1),(4,5,1,2),(3,4,5,6)); [1 2 1 1] [ ] [1 2 3 1] mat5 := [ ] [4 5 1 2] [ ] [3 4 5 6] mat6 := mat((i+1,i+2,i+3),(4,5,2),(1,i,0)); [i + 1 i + 2 i + 3] [ ] mat6 := [ 4 5 2 ] [ ] [ 1 i 0 ] mat7 := mat((1,1,0),(1,3,1),(0,1,1)); [1 1 0] [ ] mat7 := [1 3 1] [ ] [0 1 1] mat8 := mat((1,3),(-4,3)); [1 3] mat8 := [ ] [-4 3] mat9 := mat((1,2,3,4),(9,8,7,6)); [1 2 3 4] mat9 := [ ] [9 8 7 6] poly := x^7+x^5+4*x^4+5*x^3+12; 7 5 4 3 poly := x + x + 4*x + 5*x + 12 poly1 := x^2+x*y^3+x*y*z^3+y*x+2+y*3; 2 3 3 poly1 := x + x*y + x*y*z + x*y + 3*y + 2 on errcont; % Basis matrix manipulations. add_columns(mat1,1,2,5*y); [1 5*y + 2 3 4 5] [ ] [2 10*y + 3 4 5 6] [ ] [3 15*y + 4 5 6 7] [ ] [4 5*(4*y + 1) 6 7 8] [ ] [5 25*y + 6 7 8 9] add_rows(mat1,1,2,x); [ 1 2 3 4 5 ] [ ] [x + 2 2*x + 3 3*x + 4 4*x + 5 5*x + 6] [ ] [ 3 4 5 6 7 ] [ ] [ 4 5 6 7 8 ] [ ] [ 5 6 7 8 9 ] add_to_columns(mat1,3,1000); [1 2 1003 4 5] [ ] [2 3 1004 5 6] [ ] [3 4 1005 6 7] [ ] [4 5 1006 7 8] [ ] [5 6 1007 8 9] add_to_columns(mat1,{1,2,3},y); [y + 1 y + 2 y + 3 4 5] [ ] [y + 2 y + 3 y + 4 5 6] [ ] [y + 3 y + 4 y + 5 6 7] [ ] [y + 4 y + 5 y + 6 7 8] [ ] [y + 5 y + 6 y + 7 8 9] add_to_rows(mat1,2,1000); [ 1 2 3 4 5 ] [ ] [1002 1003 1004 1005 1006] [ ] [ 3 4 5 6 7 ] [ ] [ 4 5 6 7 8 ] [ ] [ 5 6 7 8 9 ] add_to_rows(mat1,{1,2,3},x); [x + 1 x + 2 x + 3 x + 4 x + 5] [ ] [x + 2 x + 3 x + 4 x + 5 x + 6] [ ] [x + 3 x + 4 x + 5 x + 6 x + 7] [ ] [ 4 5 6 7 8 ] [ ] [ 5 6 7 8 9 ] augment_columns(mat1,2); [2] [ ] [3] [ ] [4] [ ] [5] [ ] [6] augment_columns(mat1,{1,2,5}); [1 2 5] [ ] [2 3 6] [ ] [3 4 7] [ ] [4 5 8] [ ] [5 6 9] stack_rows(mat1,3); [3 4 5 6 7] stack_rows(mat1,{1,3,5}); [1 2 3 4 5] [ ] [3 4 5 6 7] [ ] [5 6 7 8 9] char_poly(mat1,x); 3 2 x *(x - 25*x - 50) column_dim(mat2); 4 row_dim(mat1); 5 copy_into(mat7,mat1,2,3); [1 2 3 4 5] [ ] [2 3 1 1 0] [ ] [3 4 1 3 1] [ ] [4 5 0 1 1] [ ] [5 6 7 8 9] copy_into(mat7,mat1,5,5); ***** Error in copy_into: the matrix [1 1 0] [ ] [1 3 1] [ ] [0 1 1] does not fit into [1 2 3 4 5] [ ] [2 3 4 5 6] [ ] [3 4 5 6 7] [ ] [4 5 6 7 8] [ ] [5 6 7 8 9] at position 5,5. diagonal(3); [3] % diagonal can take both a list of arguments or just the arguments. diagonal({mat2,mat6}); [1 1 1 1 0 0 0 ] [ ] [2 2 2 2 0 0 0 ] [ ] [3 3 3 3 0 0 0 ] [ ] [4 4 4 4 0 0 0 ] [ ] [0 0 0 0 i + 1 i + 2 i + 3] [ ] [0 0 0 0 4 5 2 ] [ ] [0 0 0 0 1 i 0 ] diagonal(mat1,mat2,mat5); [1 2 3 4 5 0 0 0 0 0 0 0 0] [ ] [2 3 4 5 6 0 0 0 0 0 0 0 0] [ ] [3 4 5 6 7 0 0 0 0 0 0 0 0] [ ] [4 5 6 7 8 0 0 0 0 0 0 0 0] [ ] [5 6 7 8 9 0 0 0 0 0 0 0 0] [ ] [0 0 0 0 0 1 1 1 1 0 0 0 0] [ ] [0 0 0 0 0 2 2 2 2 0 0 0 0] [ ] [0 0 0 0 0 3 3 3 3 0 0 0 0] [ ] [0 0 0 0 0 4 4 4 4 0 0 0 0] [ ] [0 0 0 0 0 0 0 0 0 1 2 1 1] [ ] [0 0 0 0 0 0 0 0 0 1 2 3 1] [ ] [0 0 0 0 0 0 0 0 0 4 5 1 2] [ ] [0 0 0 0 0 0 0 0 0 3 4 5 6] extend(mat1,3,2,x); [1 2 3 4 5 x x] [ ] [2 3 4 5 6 x x] [ ] [3 4 5 6 7 x x] [ ] [4 5 6 7 8 x x] [ ] [5 6 7 8 9 x x] [ ] [x x x x x x x] [ ] [x x x x x x x] [ ] [x x x x x x x] find_companion(mat5,x); 2 x - 2*x - 2 get_columns(mat1,1); { [1] [ ] [2] [ ] [3] [ ] [4] [ ] [5] } get_columns(mat1,{1,2}); { [1] [ ] [2] [ ] [3] [ ] [4] [ ] [5] , [2] [ ] [3] [ ] [4] [ ] [5] [ ] [6] } get_rows(mat1,3); { [3 4 5 6 7] } get_rows(mat1,{1,3}); { [1 2 3 4 5] , [3 4 5 6 7] } hermitian_tp(mat6); [ - i + 1 4 1 ] [ ] [ - i + 2 5 - i] [ ] [ - i + 3 2 0 ] % matrix_augment and matrix_stack can take both a list of arguments % or just the arguments. matrix_augment({mat1,mat2}); ***** Error in matrix_augment: ***** all input matrices must have the same row dimension. matrix_augment(mat4,mat2,mat4); [3 3 1 1 1 1 3 3] [ ] [4 4 2 2 2 2 4 4] [ ] [5 5 3 3 3 3 5 5] [ ] [6 6 4 4 4 4 6 6] matrix_stack(mat1,mat2); ***** Error in matrix_stack: ***** all input matrices must have the same column dimension. matrix_stack({mat6,mat((z,z,z)),mat7}); [i + 1 i + 2 i + 3] [ ] [ 4 5 2 ] [ ] [ 1 i 0 ] [ ] [ z z z ] [ ] [ 1 1 0 ] [ ] [ 1 3 1 ] [ ] [ 0 1 1 ] minor(mat1,2,3); [1 2 4 5] [ ] [3 4 6 7] [ ] [4 5 7 8] [ ] [5 6 8 9] mult_columns(mat1,3,y); [1 2 3*y 4 5] [ ] [2 3 4*y 5 6] [ ] [3 4 5*y 6 7] [ ] [4 5 6*y 7 8] [ ] [5 6 7*y 8 9] mult_columns(mat1,{2,3,4},100); [1 200 300 400 5] [ ] [2 300 400 500 6] [ ] [3 400 500 600 7] [ ] [4 500 600 700 8] [ ] [5 600 700 800 9] mult_rows(mat1,2,x); [ 1 2 3 4 5 ] [ ] [2*x 3*x 4*x 5*x 6*x] [ ] [ 3 4 5 6 7 ] [ ] [ 4 5 6 7 8 ] [ ] [ 5 6 7 8 9 ] mult_rows(mat1,{1,3,5},10); [10 20 30 40 50] [ ] [2 3 4 5 6 ] [ ] [30 40 50 60 70] [ ] [4 5 6 7 8 ] [ ] [50 60 70 80 90] pivot(mat1,3,3); [ - 4 - 2 2 4 ] [------ ------ 0 --- --- ] [ 5 5 5 5 ] [ ] [ - 2 - 1 1 2 ] [------ ------ 0 --- --- ] [ 5 5 5 5 ] [ ] [ 3 4 5 6 7 ] [ ] [ 2 1 - 1 - 2 ] [ --- --- 0 ------ ------] [ 5 5 5 5 ] [ ] [ 4 2 - 2 - 4 ] [ --- --- 0 ------ ------] [ 5 5 5 5 ] rows_pivot(mat1,3,3,{1,5}); [ - 4 - 2 2 4 ] [------ ------ 0 --- --- ] [ 5 5 5 5 ] [ ] [ 2 3 4 5 6 ] [ ] [ 3 4 5 6 7 ] [ ] [ 4 5 6 7 8 ] [ ] [ 4 2 - 2 - 4 ] [ --- --- 0 ------ ------] [ 5 5 5 5 ] remove_columns(mat1,3); [1 2 4 5] [ ] [2 3 5 6] [ ] [3 4 6 7] [ ] [4 5 7 8] [ ] [5 6 8 9] remove_columns(mat1,{2,3,4}); [1 5] [ ] [2 6] [ ] [3 7] [ ] [4 8] [ ] [5 9] remove_rows(mat1,2); [1 2 3 4 5] [ ] [3 4 5 6 7] [ ] [4 5 6 7 8] [ ] [5 6 7 8 9] remove_rows(mat1,{1,3}); [2 3 4 5 6] [ ] [4 5 6 7 8] [ ] [5 6 7 8 9] remove_rows(mat1,{1,2,3,4,5}); ***** Warning in remove_rows: all the rows have been removed. Returning nil. swap_columns(mat1,2,4); [1 4 3 2 5] [ ] [2 5 4 3 6] [ ] [3 6 5 4 7] [ ] [4 7 6 5 8] [ ] [5 8 7 6 9] swap_rows(mat1,1,2); [2 3 4 5 6] [ ] [1 2 3 4 5] [ ] [3 4 5 6 7] [ ] [4 5 6 7 8] [ ] [5 6 7 8 9] swap_entries(mat1,{1,1},{5,5}); [9 2 3 4 5] [ ] [2 3 4 5 6] [ ] [3 4 5 6 7] [ ] [4 5 6 7 8] [ ] [5 6 7 8 1] % Constructors - functions that create matrices. band_matrix(x,5); [x 0 0 0 0] [ ] [0 x 0 0 0] [ ] [0 0 x 0 0] [ ] [0 0 0 x 0] [ ] [0 0 0 0 x] band_matrix({x,y,z},6); [y z 0 0 0 0] [ ] [x y z 0 0 0] [ ] [0 x y z 0 0] [ ] [0 0 x y z 0] [ ] [0 0 0 x y z] [ ] [0 0 0 0 x y] block_matrix(1,2,{mat1,mat2}); ***** Error in block_matrix: row dimensions of ***** matrices into block_matrix are not compatible block_matrix(2,3,{mat2,mat3,mat2,mat3,mat2,mat2}); [1 1 1 1 x 1 1 1 1] [ ] [2 2 2 2 x 2 2 2 2] [ ] [3 3 3 3 x 3 3 3 3] [ ] [4 4 4 4 x 4 4 4 4] [ ] [x 1 1 1 1 1 1 1 1] [ ] [x 2 2 2 2 2 2 2 2] [ ] [x 3 3 3 3 3 3 3 3] [ ] [x 4 4 4 4 4 4 4 4] char_matrix(mat1,x); [x - 1 -2 -3 -4 -5 ] [ ] [ -2 x - 3 -4 -5 -6 ] [ ] [ -3 -4 x - 5 -6 -7 ] [ ] [ -4 -5 -6 x - 7 -8 ] [ ] [ -5 -6 -7 -8 x - 9] cfmat := coeff_matrix({x+y+4*z=10,y+x-z=20,x+y+4}); cfmat := { [4 1 1] [ ] [-1 1 1] [ ] [0 1 1] , [z] [ ] [y] [ ] [x] , [10] [ ] [20] [ ] [-4] } first cfmat * second cfmat; [x + y + 4*z] [ ] [ x + y - z ] [ ] [ x + y ] third cfmat; [10] [ ] [20] [ ] [-4] companion(poly,x); [0 0 0 0 0 0 -12] [ ] [1 0 0 0 0 0 0 ] [ ] [0 1 0 0 0 0 0 ] [ ] [0 0 1 0 0 0 -5 ] [ ] [0 0 0 1 0 0 -4 ] [ ] [0 0 0 0 1 0 -1 ] [ ] [0 0 0 0 0 1 0 ] hessian(poly1,{w,x,y,z}); [0 0 0 0 ] [ ] [ 2 3 2 ] [0 2 3*y + z + 1 3*y*z ] [ ] [ 2 3 2 ] [0 3*y + z + 1 6*x*y 3*x*z ] [ ] [ 2 2 ] [0 3*y*z 3*x*z 6*x*y*z] hilbert(4,1); [ 1 1 1 ] [ 1 --- --- ---] [ 2 3 4 ] [ ] [ 1 1 1 1 ] [--- --- --- ---] [ 2 3 4 5 ] [ ] [ 1 1 1 1 ] [--- --- --- ---] [ 3 4 5 6 ] [ ] [ 1 1 1 1 ] [--- --- --- ---] [ 4 5 6 7 ] hilbert(3,y+x); [ - 1 - 1 - 1 ] [----------- ----------- -----------] [ x + y - 2 x + y - 3 x + y - 4 ] [ ] [ - 1 - 1 - 1 ] [----------- ----------- -----------] [ x + y - 3 x + y - 4 x + y - 5 ] [ ] [ - 1 - 1 - 1 ] [----------- ----------- -----------] [ x + y - 4 x + y - 5 x + y - 6 ] jacobian({x^4,x*y^2,x*y*z^3},{w,x,y,z}); [ 3 ] [0 4*x 0 0 ] [ ] [ 2 ] [0 y 2*x*y 0 ] [ ] [ 3 3 2] [0 y*z x*z 3*x*y*z ] jordan_block(x,5); [x 1 0 0 0] [ ] [0 x 1 0 0] [ ] [0 0 x 1 0] [ ] [0 0 0 x 1] [ ] [0 0 0 0 x] make_identity(11); [1 0 0 0 0 0 0 0 0 0 0] [ ] [0 1 0 0 0 0 0 0 0 0 0] [ ] [0 0 1 0 0 0 0 0 0 0 0] [ ] [0 0 0 1 0 0 0 0 0 0 0] [ ] [0 0 0 0 1 0 0 0 0 0 0] [ ] [0 0 0 0 0 1 0 0 0 0 0] [ ] [0 0 0 0 0 0 1 0 0 0 0] [ ] [0 0 0 0 0 0 0 1 0 0 0] [ ] [0 0 0 0 0 0 0 0 1 0 0] [ ] [0 0 0 0 0 0 0 0 0 1 0] [ ] [0 0 0 0 0 0 0 0 0 0 1] on rounded; % makes things a bit easier to read. random_matrix(3,3,100); [ - 28.8708957912 - 83.4775707772 - 22.5350939803] [ ] [ - 83.4775707772 35.4030294062 5.76473742166 ] [ ] [ - 22.5350939803 5.76473742166 74.9835325987 ] on not_negative; random_matrix(3,3,100); [52.1649704637 63.7238994877 0.886186275566] [ ] [63.7238994877 18.6365765396 77.194379299 ] [ ] [0.886186275566 77.194379299 4.09332974714 ] on only_integer; random_matrix(3,3,100); [49 42 67] [ ] [42 69 84] [ ] [67 84 31] on symmetric; random_matrix(3,3,100); [73 63 85] [ ] [63 5 81] [ ] [85 81 35] off symmetric; on upper_matrix; random_matrix(3,3,100); [4 70 28] [ ] [0 7 58] [ ] [0 0 19] off upper_matrix; on lower_matrix; random_matrix(3,3,100); [7 0 0 ] [ ] [46 49 0 ] [ ] [42 70 65] off lower_matrix; on imaginary; off not_negative; random_matrix(3,3,100); [ 38*i - 14 - 72*i + 19 51*i - 30 ] [ ] [ - 99*i + 72 94*i - 59 3*i - 46 ] [ ] [ 47*i - 54 - 9*i - 73*i - 28] off rounded; % toeplitz and vandermonde can take both a list of arguments or just % the arguments. toeplitz({1,2,3,4,5}); [1 2 3 4 5] [ ] [2 1 2 3 4] [ ] [3 2 1 2 3] [ ] [4 3 2 1 2] [ ] [5 4 3 2 1] toeplitz(x,y,z); [x y z] [ ] [y x y] [ ] [z y x] vandermonde({1,2,3,4,5}); [1 1 1 1 1 ] [ ] [1 2 4 8 16 ] [ ] [1 3 9 27 81 ] [ ] [1 4 16 64 256] [ ] [1 5 25 125 625] vandermonde(x,y,z); [ 2] [1 x x ] [ ] [ 2] [1 y y ] [ ] [ 2] [1 z z ] % kronecker_product a1 := mat((1,2),(3,4),(5,6)); [1 2] [ ] a1 := [3 4] [ ] [5 6] a2 := mat((1,x,1),(2,2,2),(3,3,3)); [1 x 1] [ ] a2 := [2 2 2] [ ] [3 3 3] kronecker_product(a1,a2); [1 x 1 2 2*x 2 ] [ ] [2 2 2 4 4 4 ] [ ] [3 3 3 6 6 6 ] [ ] [3 3*x 3 4 4*x 4 ] [ ] [6 6 6 8 8 8 ] [ ] [9 9 9 12 12 12] [ ] [5 5*x 5 6 6*x 6 ] [ ] [10 10 10 12 12 12] [ ] [15 15 15 18 18 18] clear a1,a2; % High level algorithms. on rounded; % makes output easier to read. ch := cholesky(mat7); ch := { [1 0 0 ] [ ] [1 1.41421356237 0 ] [ ] [0 0.707106781187 0.707106781187] , [1 1 0 ] [ ] [0 1.41421356237 0.707106781187] [ ] [0 0 0.707106781187] } tp first ch - second ch; [0 0 0] [ ] [0 0 0] [ ] [0 0 0] tmp := first ch * second ch; [1 1 0] [ ] tmp := [1 3.0 1] [ ] [0 1 1] tmp - mat7; [0 0 0] [ ] [0 0 0] [ ] [0 0 0] off rounded; gram_schmidt({1,0,0},{1,1,0},{1,1,1}); {{1,0,0},{0,1,0},{0,0,1}} gram_schmidt({1,2},{3,4}); 1 2 2*sqrt(5) - sqrt(5) {{---------,---------},{-----------,------------}} sqrt(5) sqrt(5) 5 5 on rounded; % again, makes large quotients a bit more readable. % The algorithm used for lu_decom sometimes swaps the rows of the input % matrix so that (given matrix A, lu_decom(A) = {L,U,vec}), we find L*U % does not equal A but a row equivalent of it. The call convert(A,vec) % will return this row equivalent (ie: L*U = convert(A,vec)). lu := lu_decom(mat5); lu := { [4 0 0 0 ] [ ] [1 0.75 0 0 ] [ ] [1 0.75 2.0 0 ] [ ] [3 0.25 4.0 4.33333333333] , [1 1.25 0.25 0.5 ] [ ] [0 1 1 0.666666666667] [ ] [0 0 1 0 ] [ ] [0 0 0 1 ] , [3,3,3,4]} mat5; [1 2 1 1] [ ] [1 2 3 1] [ ] [4 5 1 2] [ ] [3 4 5 6] tmp := first lu * second lu; [4 5.0 1 2.0] [ ] [1 2.0 1 1 ] tmp := [ ] [1 2.0 3.0 1 ] [ ] [3 4.0 5.0 6.0] tmp1 := convert(mat5,third lu); [4 5 1 2] [ ] [1 2 1 1] tmp1 := [ ] [1 2 3 1] [ ] [3 4 5 6] tmp - tmp1; [0 0 0 0] [ ] [0 0 0 0] [ ] [0 0 0 0] [ ] [0 0 0 0] % and the complex case... lu1 := lu_decom(mat6); lu1 := { [ 1 0 0 ] [ ] [ 4 - 4*i + 5 0 ] [ ] [i + 1 3 0.414634146341*i + 2.26829268293] , [1 i 0 ] [ ] [0 1 0.19512195122*i + 0.243902439024] [ ] [0 0 1 ] , [3,2,3]} mat6; [i + 1 i + 2 i + 3] [ ] [ 4 5 2 ] [ ] [ 1 i 0 ] tmp := first lu1 * second lu1; [ 1 i 0 ] [ ] tmp := [ 4 5 2.0 ] [ ] [i + 1 i + 2 i + 3.0] tmp1 := convert(mat6,third lu1); [ 1 i 0 ] [ ] tmp1 := [ 4 5 2 ] [ ] [i + 1 i + 2 i + 3] tmp - tmp1; [0 0 0] [ ] [0 0 0] [ ] [0 0 0] mat9inv := pseudo_inverse(mat9); [ - 0.199999999996 0.100000000013 ] [ ] [ - 0.0499999999988 0.0500000000037 ] mat9inv := [ ] [ 0.0999999999982 - 5.57816640101e-12] [ ] [ 0.249999999995 - 0.0500000000148 ] mat9 * mat9inv; [ 0.999999999982 - 0.0000000000557817125824] [ ] [5.54201129432e-12 1.00000000002 ] simplex(min,2*x1+14*x2+36*x3,{-2*x1+x2+4*x3>=5,-x1-2*x2-3*x3<=2}); {45.0,{x1=0,x2=0,x3=1.25}} simplex(max,10000 x1 + 1000 x2 + 100 x3 + 10 x4 + x5,{ x1 <= 1, 20 x1 + x2 <= 100, 200 x1 + 20 x2 + x3 <= 10000, 2000 x1 + 200 x2 + 20 x3 + x4 <= 1000000, 20000 x1 + 2000 x2 + 200 x3 + 20 x4 + x5 <= 100000000}); {1.0e+8,{x1=0,x2=0,x3=0,x4=0,x5=100000000.0}} simplex(max, 5 x1 + 4 x2 + 3 x3, { 2 x1 + 3 x2 + x3 <= 5, 4 x1 + x2 + 2 x3 <= 11, 3 x1 + 4 x2 + 2 x3 <= 8 }); {13.0,{x1=2.0,x2=0,x3=1.0}} simplex(min,3 x1 + 5 x2,{ x1 + 2 x2 >= 2, 22 x1 + x2 >= 3}); {5.04651162791,{x1=0.093023255813953,x2=0.95348837209302}} simplex(max,10x+5y+5.5z,{5x+3z<=200,0.2x+0.1y+0.5z<=12,0.1x+0.2y+0.3z<=9, 30x+10y+50z<=1500}); {525.0,{x=40.0,y=25.0,z=0}} %example of extra variables (>=0) being added. simplex(min,x-y,{x>=-3}); *** Warning: variable y not defined in input. Has been defined as >=0. ***** Error in simplex: The problem is unbounded. % unfeasible as simplex algorithm implies all x>=0. simplex(min,x,{x<=-100}); ***** Error in simplex: Problem has no feasible solution. % three error examples. simplex(maxx,x,{x>=5}); ***** Error in simplex(first argument): must be either max or min. simplex(max,x,x>=5); ***** Error in simplex(third argument}: must be a list. simplex(max,x,{x<=y}); ***** Error in simplex(third argument): ***** right hand side of each inequality must be a number simplex(max, 346 X11 + 346 X12 + 248 X21 + 248 X22 + 399 X31 + 399 X32 + 200 Y11 + 200 Y12 + 75 Y21 + 75 Y22 + 2.35 Z1 + 3.5 Z2, { 4 X11 + 4 X12 + 2 X21 + 2 X22 + X31 + X32 + 250 Y11 + 250 Y12 + 125 Y21 + 125 Y22 <= 25000, X11 + X12 + X21 + X22 + X31 + X32 + 2 Y11 + 2 Y12 + Y21 + Y22 <= 300, 20 X11 + 15 X12 + 30 Y11 + 20 Y21 + Z1 <= 1500, 40 X12 + 35 X22 + 50 X32 + 15 Y12 + 10 Y22 + Z2 = 5000, X31 = 0, Y11 + Y12 <= 50, Y21 + Y22 <= 100 }); {99250.0, {y21=0, y22=0, x31=0, x11=75.0, z1=0, x21=225.0, z2=5000.0, x32=0, x22=0, x12=0, y12=0, y11=0}} % from Marc van Dongen. Finding the first feasible solution for the % solution of systems of linear diophantine inequalities. simplex(max,0,{ 3*X259+4*X261+3*X262+2*X263+X269+2*X270+3*X271+4*X272+5*X273+X229=2, 7*X259+11*X261+8*X262+5*X263+3*X269+6*X270+9*X271+12*X272+15*X273+X229=4, 2*X259+5*X261+4*X262+3*X263+3*X268+4*X269+5*X270+6*X271+7*X272+8*X273=1, X262+2*X263+5*X268+4*X269+3*X270+2*X271+X272+2*X229=1, X259+X262+2*X263+4*X268+3*X269+2*X270+X271-X273+3*X229=2, X259+2*X261+2*X262+2*X263+3*X268+3*X269+3*X270+3*X271+3*X272+3*X273+X229=1, X259+X261+X262+X263+X268+X269+X270+X271+X272+X273+X229=1}); {0, {x229=0.5, x259=0.5, x261=0, x262=0, x263=0, x268=0, x269=0, x270=0, x271=0, x272=0, x273=0}} svd_ans := svd(mat8); svd_ans := { [ 0.289784137735 0.957092029805] [ ] [ - 0.957092029805 0.289784137735] , [5.1491628629 0 ] [ ] [ 0 2.9130948854] , [ - 0.687215403194 0.726453707825 ] [ ] [ - 0.726453707825 - 0.687215403194] } tmp := tp first svd_ans * second svd_ans * third svd_ans; [ 0.99999998509 2.9999999859 ] tmp := [ ] [ - 4.00000004924 2.99999995342] tmp - mat8; [ - 0.0000000149096008872 - 0.0000000141042817425] [ ] [ - 0.0000000492430629606 - 0.0000000465832750152] mat9inv := pseudo_inverse(mat9); [ - 0.199999999996 0.100000000013 ] [ ] [ - 0.0499999999988 0.0500000000037 ] mat9inv := [ ] [ 0.0999999999982 - 5.57816640101e-12] [ ] [ 0.249999999995 - 0.0500000000148 ] mat9 * mat9inv; [ 0.999999999982 - 0.0000000000557817125824] [ ] [5.54201129432e-12 1.00000000002 ] % Predicates. matrixp(mat1); t matrixp(poly); squarep(mat2); t squarep(mat3); symmetricp(mat1); t symmetricp(mat3); if not rounded_was_on then off rounded; END; (TIME: linalg 3359 3469)