> read("/afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlgebra/Maple/
> initlib.mpl");

  libname := /afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlg\

        ebra/Maple, "/afs/bp.ncsu.edu/dist/maple551/update",

        "/afs/bp.ncsu.edu/dist/maple551/lib"

> with(refpkg);

       [E_I, E_II, E_III, mydet, myinverse, mysolve, ref, xref]

> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> A := matrix(3,4);

                    A := array(1 .. 3, 1 .. 4, [])

> E1 := E_I(3,2,3);

                               [1    0    0]
                               [           ]
                         E1 := [0    0    1]
                               [           ]
                               [0    1    0]

> E1p := E_I(4,2,4);

                             [1    0    0    0]
                             [                ]
                             [0    0    0    1]
                      E1p := [                ]
                             [0    0    1    0]
                             [                ]
                             [0    1    0    0]

> evalm(E1 &* A);

              [A[1, 1]    A[1, 2]    A[1, 3]    A[1, 4]]
              [                                        ]
              [A[3, 1]    A[3, 2]    A[3, 3]    A[3, 4]]
              [                                        ]
              [A[2, 1]    A[2, 2]    A[2, 3]    A[2, 4]]

> evalm(A &* E1p);

              [A[1, 1]    A[1, 4]    A[1, 3]    A[1, 2]]
              [                                        ]
              [A[2, 1]    A[2, 4]    A[2, 3]    A[2, 2]]
              [                                        ]
              [A[3, 1]    A[3, 4]    A[3, 3]    A[3, 2]]

> E2 := E_II(3,3,alpha);

                             [1    0      0  ]
                             [               ]
                       E2 := [0    1      0  ]
                             [               ]
                             [0    0    alpha]

> evalm(E2 &* A);

        [A[1, 1] , A[1, 2] , A[1, 3] , A[1, 4]]

        [A[2, 1] , A[2, 2] , A[2, 3] , A[2, 4]]

        [alpha A[3, 1] , alpha A[3, 2] , alpha A[3, 3] ,

        alpha A[3, 4]]

> E2p := E_II(4,1,beta);

                            [beta    0    0    0]
                            [                   ]
                            [ 0      1    0    0]
                     E2p := [                   ]
                            [ 0      0    1    0]
                            [                   ]
                            [ 0      0    0    1]

> evalm(A &* E2p);

           [A[1, 1] beta    A[1, 2]    A[1, 3]    A[1, 4]]
           [                                             ]
           [A[2, 1] beta    A[2, 2]    A[2, 3]    A[2, 4]]
           [                                             ]
           [A[3, 1] beta    A[3, 2]    A[3, 3]    A[3, 4]]

> E3 := E_III(3,1,3,alpha);

                             [  1      0    0]
                             [               ]
                       E3 := [  0      1    0]
                             [               ]
                             [alpha    0    1]

> evalm(E3 &* A);

        [A[1, 1] , A[1, 2] , A[1, 3] , A[1, 4]]

        [A[2, 1] , A[2, 2] , A[2, 3] , A[2, 4]]

        [alpha A[1, 1] + A[3, 1] , alpha A[1, 2] + A[3, 2] ,

        alpha A[1, 3] + A[3, 3] , alpha A[1, 4] + A[3, 4]]

> E3p := E_III(4,2,4,beta);

                            [1     0      0    0]
                            [                   ]
                            [0     1      0    0]
                     E3p := [                   ]
                            [0     0      1    0]
                            [                   ]
                            [0    beta    0    1]

> evalm(A &* E3p);

        [A[1, 1] , A[1, 2] + A[1, 4] beta , A[1, 3] , A[1, 4]]
        [                                                    ]
        [A[2, 1] , A[2, 2] + A[2, 4] beta , A[2, 3] , A[2, 4]]
        [                                                    ]
        [A[3, 1] , A[3, 2] + A[3, 4] beta , A[3, 3] , A[3, 4]]

> E3pp := E_III(4,4,2,beta);

                            [1    0    0     0  ]
                            [                   ]
                            [0    1    0    beta]
                    E3pp := [                   ]
                            [0    0    1     0  ]
                            [                   ]
                            [0    0    0     1  ]

> evalm(A &* E3pp);

        [A[1, 1] , A[1, 2] , A[1, 3] , A[1, 2] beta + A[1, 4]]
        [                                                    ]
        [A[2, 1] , A[2, 2] , A[2, 3] , A[2, 2] beta + A[2, 4]]
        [                                                    ]
        [A[3, 1] , A[3, 2] , A[3, 3] , A[3, 2] beta + A[3, 4]]

> print(E1);

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

> inverse(E1);

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

> print(E2);

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

> inverse(E2);

                          [1    0      0  ]
                          [               ]
                          [0    1      0  ]
                          [               ]
                          [            1  ]
                          [0    0    -----]
                          [          alpha]

> print(E3);

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

> inverse(E3);

                          [  1       0    0]
                          [                ]
                          [  0       1    0]
                          [                ]
                          [-alpha    0    1]

> rnd := rand(-10..10);

  rnd := proc()
local t;
global _seed;
    _seed := irem(427419669081*_seed, 999999999989);
    t := _seed;
    irem(t, 21) - 10
end

> B := matrix(4,3);

                    B := array(1 .. 4, 1 .. 3, [])

> for i from 1 to 4 do for j from 1 to 3 do B[i,j]:=rnd(); od;od;
> print(B);

                           [-4     7     8]
                           [              ]
                           [10    -6    -8]
                           [              ]
                           [-5     7     6]
                           [              ]
                           [-6     0     5]

> U := xref(B, 3, 'T');

                            [-4     7       8 ]
                            [                 ]
                            [ 0    23/2    12 ]
                            [                 ]
                       U := [              -50]
                            [ 0     0      ---]
                            [              23 ]
                            [                 ]
                            [ 0     0       0 ]

> print(T);

                      [ 1       0      0     0]
                      [                       ]
                      [5/2      1      0     0]
                      [                       ]
                      [-20                    ]
                      [---     7/46    1     0]
                      [23                     ]
                      [                       ]
                      [        119     91     ]
                      [-4/5    ---     --    1]
                      [        100     50     ]

> evalm(T &* B);

                         [-4     7       8 ]
                         [                 ]
                         [ 0    23/2    12 ]
                         [                 ]
                         [              -50]
                         [ 0     0      ---]
                         [              23 ]
                         [                 ]
                         [ 0     0       0 ]

> xref(B, 3, 'T', 'Tinv');

                         [-4     7       8 ]
                         [                 ]
                         [ 0    23/2    12 ]
                         [                 ]
                         [              -50]
                         [ 0     0      ---]
                         [              23 ]
                         [                 ]
                         [ 0     0       0 ]

> print(Tinv);

                      [ 1       0      0     0]
                      [                       ]
                      [-5/2     1      0     0]
                      [                       ]
                      [        -7             ]
                      [5/4     --      1     0]
                      [        46             ]
                      [                       ]
                      [        -21    -91     ]
                      [3/2     ---    ---    1]
                      [        23     50      ]

>