@@ -1,1687 +1,1687 @@ -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) +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)