reduce-historical/r37/packages/linalg/linalg.rlg

Fri Feb 12 00:06:04 MET 1999
REDUCE 3.7, 15-Jan-99 ...

1: 1: 
2: 2: 2: 2: 2: 2: 2: 2: 2: 
3: 3: 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);


[ - 8.11911717343   - 75.7167729277    30.62058083   ]
[                                                    ]
[ - 50.0325962624   47.1655452861     35.8674263384  ]
[                                                    ]
[ - 49.3715543826   - 97.5563670864   - 18.8861862756]


on not_negative;


random_matrix(3,3,100);


[43.8999853223  33.7140980286   33.75065406 ]
[                                           ]
[49.7333355117  98.9642944905  58.5331568816]
[                                           ]
[39.9146060895  67.7954727837  24.8684367642]


on only_integer;


random_matrix(3,3,100);


[16  77  49]
[          ]
[28  84  51]
[          ]
[84  56  57]


on symmetric;


random_matrix(3,3,100);


[89  74  91]
[          ]
[74  95  41]
[          ]
[91  41  87]


off symmetric;


on upper_matrix;


random_matrix(3,3,100);


[41  3   8 ]
[          ]
[0   31  80]
[          ]
[0   0   12]


off upper_matrix;


on lower_matrix;


random_matrix(3,3,100);


[69  0   0 ]
[          ]
[34  87  0 ]
[          ]
[78  72  13]


off lower_matrix;


on imaginary;


off not_negative;


random_matrix(3,3,100);


[ - 95*i - 72   - 57*i + 59  52*i + 46]
[                                     ]
[ - 40*i - 54      70*i      39*i + 28]
[                                     ]
[ - 40*i + 45   28*i - 81    9*i + 74 ]


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.5409010713e-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});


{100000000,{x1=0,x2=0,x3=0,x4=0,x5=1.0e+08}}


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.0930233,x2=0.953488}}


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.0000000149096013313   - 0.0000000141042812984]
[                                                  ]
[ - 0.0000000492430638488   - 0.0000000465832745711]



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.5409010713e-12         1.00000000002       ]



% triang_adjoint(in_mat) calculates the
% triangularizing adjoint of in_mat

triang_adjoint(mat1);


[1   0  0   0  0]
[               ]
[-2  1  0   0  0]
[               ]
[-1  2  -1  0  0]
[               ]
[0   0  0   0  0]
[               ]
[0   0  0   0  0]


triang_adjoint(mat2);


[1   0  0  0]
[           ]
[-2  1  0  0]
[           ]
[0   0  0  0]
[           ]
[0   0  0  0]


triang_adjoint(mat5);


[1    0   0   0]
[              ]
[-1   1   0   0]
[              ]
[-3   3   0   0]
[              ]
[10  -12  -4  6]


triang_adjoint(mat6);


[   1       0      0  ]
[                     ]
[  -4     i + 1    0  ]
[                     ]
[4*i - 5    3    i - 3]


triang_adjoint(mat7);


[1    0    0]
[           ]
[-1   1    0]
[           ]
[1    - 1  2]


triang_adjoint(mat8);


[1  0]
[    ]
[4  1]


triang_adjoint(mat9);


***** Error in triang_adjoint: input matrix should be square. 


% testing triang_adjoint with random matrices

% the range of the integers is in one case from
% -1000 to 1000. in the other case it is from
% -1 to 1 so that the deteminant of the i-th
% submatrix equals very often to zero.

% random matrix contains arbitrary real values
off only_integer;


tmp:=random_matrix(5,5,1000);


tmp := mat(( - 558.996086656*i + 1.71931812953,76.9987188735*i + 1.19004104683,

             - 739.283479439*i - 1.32534106204,742.101952123*i + 1.35926854848,

            680.515777254*i + 1.56403177895),

           ( - 689.196170962*i + 1.49289170118,

             - 232.584493916*i - 1.38227180395,280.109305836*i + 1.38865247861,

            298.151479065*i - 1.19035182389, - 602.312143386*i - 1.82183796879),

           ( - 739.195910955*i - 1.45944960213,859.293884826*i + 1.70488070379,

            359.856032683*i - 1.2966991869, - 105.409833087*i - 1.9360631701,

            234.350529301*i - 1.15598520849),

           (155.629059348*i + 1.09264385739, - 16.1559469072*i - 1.9425176505,

            725.11578405*i - 1.05760723025,783.020942195*i - 1.28625265346,

             - 544.129360355*i + 1.74790906085),

           (373.562370318*i - 1.95218354686, - 722.109349973*i + 1.56309793677,

             - 746.664959169*i - 1.9915755693,186.154794517*i - 1.09842189916,

            435.90998986*i - 1.46175649496))


triang_adjoint tmp;


mat((1,0,0,0,0),

    (689.196170962*i - 1.49289170118, - 558.996086656*i + 1.71931812953,0,0,0),

    ( - 1253.37955588*i + 7.64148089854e+5, - 1516.42713845*i - 4.23429448803e+5

     ,1078.01877642*i - 1.830851973e+5,0,0),

      102791325687.0*i + 7.3752778526e+8
    (------------------------------------,
              i - 169.834887206

       - 3.66748178757e+10*i - 6.62162769101e+6
     -------------------------------------------,
                  i - 169.834887206

     9.85957444629e+7*i - 1.01033337718e+6,

      - 7.49414742893e+8*i - 2.25311577415e+6,0),

       - 547052849318.0*i + 4.06181988045e+13
    (-----------------------------------------,
                 i - 112.974983172

       - 141265342333.0*i + 4.13350560163e+12
     -----------------------------------------,
                 i - 112.974983172

      845804392649.0*i - 9.62488227345e+13
     --------------------------------------,
               i - 112.974983172

      876106032577.0*i - 2.66464796763e+13
     --------------------------------------,
               i - 112.974983172

      1.47617976407e+12*i - 1.66771384004e+14
     -----------------------------------------))
                 i - 169.834887206



tmp:=random_matrix(1,1,1000);


tmp := [ - 463.860434427*i + 1.35500571348]


triang_adjoint tmp;


[1]



% random matrix contains complex real values
on imaginary;


tmp:=random_matrix(5,5,1000);


tmp := mat((107.345792105*i - 1.98704739339,188.868545358*i + 1.22417796742,

             - 630.485915434*i + 1.32195292724,

             - 542.510039297*i - 1.94318764036,359.860945563*i - 1.69174206177),

           ( - 469.501213476*i - 1.17375946319, - 62.2197820375*i - 1.4051578261

            , - 98.6604380996*i + 1.64691610034,

             - 216.296595937*i + 1.56809020199,797.19877393*i - 1.31894550853),

           (2.07054207792*i + 1.3891068942,393.038868455*i - 1.60894498437,

             - 215.390393738*i - 1.00068640594,

             - 195.674928032*i + 1.22123114986,211.921323796*i - 1.42499533273),

           ( - 750.357435524*i - 1.19871674827,

             - 792.333836712*i - 1.63151974094, - 494.87049225*i + 1.99554801527

            ,638.173945822*i + 1.23793954612,111.418959978*i - 1.96273029328),

           ( - 255.359922267*i + 1.99035939892,

             - 575.376389757*i - 1.03533681609,463.961589382*i - 1.86476410547,

            83.8856338571*i + 1.10369785887, - 129.597812786*i - 1.4917934624))


triang_adjoint tmp;


mat((1,0,0,0,0),

    (469.501213476*i + 1.17375946319,107.345792105*i - 1.98704739339,0,0,0),

    (383.407897912*i + 1.84407237435e+5,1218.59364331*i + 41798.5118562,

     769.235159465*i - 81990.7504399,0,0),

       - 1.411092405e+10*i - 1.91497165215e+8
    (-----------------------------------------,
                 i - 106.587367245

       - 2.06157034475e+10*i + 1.09218575639e+8
     -------------------------------------------,
                  i - 106.587367245

      - 2.4008888901e+8*i + 13175.2571592,

      - 1.02728261373e+8*i + 9.22309484944e+5,0),

       - 203213290519.0*i - 3.07405185302e+12
    (-----------------------------------------,
                 i - 184.764270765

      311149245317.0*i + 2.05618234856e+13
     --------------------------------------,
               i - 184.764270765

      212889617996.0*i - 4.13210409411e+13
     --------------------------------------,
               i - 184.764270765

       - 7.79955805661e+10*i - 5.10418442965e+12
     --------------------------------------------,
                  i - 184.764270765

      7.62835257557e+10*i - 1.40944700076e+13
     -----------------------------------------))
                 i - 106.587367245



tmp:=random_matrix(1,1,1000);


tmp := [276.278111177*i + 1.74724262616]


triang_adjoint tmp;


[1]


off imaginary;



% random matrix contains rounded real values
on rounded;


tmp:=random_matrix(5,5,1000);


tmp := mat(( - 293.224093687, - 99.023221037, - 819.400851656,796.020234589,

            593.862087611),

           ( - 137.84203019,354.3234619, - 852.314261681, - 217.485901759,

            256.139775139),

           (324.37828726, - 56.5718498235, - 118.293003834,108.279501424,

            23.2385400299),

           ( - 976.556156754,684.207160793,146.328625386,502.848132905,

            312.766816689),

           (211.783458501,166.556239469,175.715904944,251.57997022,280.123720131

            ))


triang_adjoint tmp;


mat((1,0,0,0,0),

    (137.84203019, - 293.224093687,0,0,0),

    ( - 1.07136859076e+5, - 48709.2122316, - 1.17545737812e+5,0,0),

    (1.27980020917e+8, - 1.64635380167e+8,4.76863677307e+8,1.43208428244e+8,0),

    (5.82963241185e+10,3.9383738062e+10, - 437637051137.0, - 111757830528.0,

     261327212376.0))



tmp:=random_matrix(1,1,1000);


tmp := [406.584701921]


triang_adjoint tmp;


[1]


off rounded;



% random matrix contains only integer values
on only_integer;


tmp:=random_matrix(7,7,1000);


       [969   210    8    244   -887  -39   -916]
       [                                        ]
       [-774  296   -475  -694  -909  560    89 ]
       [                                        ]
       [-390  -559  -551  -567  241   -306  -655]
       [                                        ]
tmp := [-478  809   181   -987  -144  929   -886]
       [                                        ]
       [188   267   -778  660   374   590    30 ]
       [                                        ]
       [ 73   971   -946  -43   -215  386   -365]
       [                                        ]
       [-792  -852  558   -797  343   219   110 ]


triang_adjoint tmp;


mat((1,0,0,0,0,0,0),

    (774,969,0,0,0,0,0),

    (548106,459771,449364,0,0,0,0),

    (-108937808,399369604,-497500435,-461605941,0,0,0),

    (-386678984240,-1001551613816,454549593485,637690866447,433944480084,0,0),

    (-604165739229705,-320961967400919,-165015285307395,-1008712187269380,

     -1670995725485274,1433408878792557,0),

    (-181830640557070260,295390292387079435,816541226477288004,

     850494616785589377,458176543109779557,-1709784109660828152,

     -1475366833406131953))



tmp:=random_matrix(7,7,1);


       [0  0  0  0  0  0  0]
       [                   ]
       [0  0  0  0  0  0  0]
       [                   ]
       [0  0  0  0  0  0  0]
       [                   ]
tmp := [0  0  0  0  0  0  0]
       [                   ]
       [0  0  0  0  0  0  0]
       [                   ]
       [0  0  0  0  0  0  0]
       [                   ]
       [0  0  0  0  0  0  0]


triang_adjoint tmp;


[1  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]
[                   ]
[0  0  0  0  0  0  0]



% random matrix contains only complex integer
% values

on imaginary;


tmp:=random_matrix(5,5,1000);


tmp := mat((12*(38*i + 83),3*(153*i - 305),2*(141*i + 427), - 553*i + 617,

            3*(83*i + 157)),

           (164*i - 635, - 133*i + 991, - 373*i + 210,965*i - 608,2*(358*i - 55)

            ),

           ( - 230*i + 227,32*i + 339,2*(485*i - 219),707*i - 767, - 985*i - 51)

            ,

           ( - 609*i + 647, - 441*i + 187,930*i - 349,250*i - 211,274*i - 451),

           ( - 374*i - 135,793*i + 592, - 81*i - 1,89*i + 26,3*( - 40*i + 201)))


triang_adjoint tmp;


mat((1,0,0,0,0),

    ( - 164*i + 635,12*(38*i + 83),0,0,0),

    (293397*i - 414880,9*(14243*i - 47243),3*(253651*i + 180645),0,0),

       - 253324472288717*i + 71265413812547
    (---------------------------------------,
                253651*i + 180645

      2*( - 220885726602145*i - 1441709355714)
     ------------------------------------------, - 1436348339*i + 1393250309,
                 253651*i + 180645

     511458435*i - 1454012933,0),

      13983048003979950612955437881*i - 71498490838832832842693585028
    (-----------------------------------------------------------------,
                  65634686423804933*i - 9174596297286164

      89295323223054915316808489269*i - 37624299403809895760446255007
     -----------------------------------------------------------------,
                  65634686423804933*i - 9174596297286164

      2*( - 71881165390656818494884812727*i - 25318671134083617432051412624)
     ------------------------------------------------------------------------,
                      65634686423804933*i - 9174596297286164

      134577377248105484011524135103*i + 3495516738012600790097438251
     -----------------------------------------------------------------,
                  65634686423804933*i - 9174596297286164

      6*(65634686423804933*i - 9174596297286164)
     --------------------------------------------))
                  253651*i + 180645



tmp:=random_matrix(5,5,2);


       [i - 1   i      i       0       - (i + 1)]
       [                                        ]
       [  0     i     -1     - i + 1    i + 1   ]
       [                                        ]
tmp := [ -1     0      0     - i + 1      -1    ]
       [                                        ]
       [ -1     - i   - i      - i      i + 1   ]
       [                                        ]
       [i - 1   0    i + 1     -1         0     ]


triang_adjoint tmp;


[      1             0             0             0             0     ]
[                                                                    ]
[      0           i - 1           0             0             0     ]
[                                                                    ]
[  - (i + 1)       i + 1                                             ]
[------------     -------      - (i + 1)         0             0     ]
[   i - 1          i - 1                                             ]
[                                                                    ]
[  - (i + 1)                  2*(2*i + 1)       - 2*i                ]
[------------        0       -------------    --------         0     ]
[     i                          i - 1         i - 1                 ]
[                                                                    ]
[ 2*(3*i - 4)    2*(i + 2)    5*(3*i + 1)     - 7*i + 1    2*(i + 2) ]
[-------------  -----------  -------------  ------------  -----------]
[   4*i + 3        i - 1        4*i + 3       4*i + 3        i - 1   ]



% 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;

4: 4: 4: 4: 4: 4: 4: 4: 4: 

Time for test: 1460 ms, plus GC time: 60 ms

5: 5: 
Quitting
Fri Feb 12 00:06:12 MET 1999