reduce-historical/r36/xlog/LINALG.LOG

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)