> 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]

> ?xref
> 
# Examples:
> infolevel['refpkg']:=10:
> with(refpkg, xref);

                                [xref]

> A := matrix([ [1, 2, 1, 3], [3, -1, -3, -1], [2, 3, 1, 4] ]);

                           [1     2     1     3]
                           [                   ]
                      A := [3    -1    -3    -1]
                           [                   ]
                           [2     3     1     4]

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

                               [1     2     1      3]
                               [                    ]
                 i=, 1, j=, 1, [0    -7    -6    -10]
                               [                    ]
                               [0    -1    -1     -2]


                              [ 1    0    0]
                              [            ]
                          L=, [-3    1    0]
                              [            ]
                              [-2    0    1]


                                 [1    0    0]
                                 [           ]
                        L^(-1)=, [3    1    0]
                                 [           ]
                                 [2    0    1]


                             [1     2     1       3  ]
                             [                       ]
               i=, 2, j=, 2, [0    -7     -6     -10 ]
                             [                       ]
                             [0     0    -1/7    -4/7]


                           [  1       0      0]
                           [                  ]
                       L=, [ -3       1      0]
                           [                  ]
                           [-11/7    -1/7    1]


                                [1     0     0]
                                [             ]
                       L^(-1)=, [3     1     0]
                                [             ]
                                [2    1/7    1]


                         [1     2     1       3  ]
                         [                       ]
                    U := [0    -7     -6     -10 ]
                         [                       ]
                         [0     0    -1/7    -4/7]

> print(T);

                         [  1       0      0]
                         [                  ]
                         [ -3       1      0]
                         [                  ]
                         [-11/7    -1/7    1]

> evalm(T &* A);

                      [1     2     1       3  ]
                      [                       ]
                      [0    -7     -6     -10 ]
                      [                       ]
                      [0     0    -1/7    -4/7]

> evalm(L &* U);

                        [1     2     1     3]
                        [                   ]
                        [3    -1    -3    -1]
                        [                   ]
                        [2     3     1     4]

> B := matrix(4,3,[1,2,3, 0,0,0, 0,0,4, 0,1,5]);

                               [1    2    3]
                               [           ]
                               [0    0    0]
                          B := [           ]
                               [0    0    4]
                               [           ]
                               [0    1    5]

> U2 := xref(B, 3, 'T2', 'L2', 'P');

                                   [1    2    3]
                                   [           ]
                                   [0    0    0]
                     i=, 1, j=, 1, [           ]
                                   [0    0    4]
                                   [           ]
                                   [0    1    5]


                            [1    0    0    0]
                            [                ]
                            [0    1    0    0]
                        L=, [                ]
                            [0    0    1    0]
                            [                ]
                            [0    0    0    1]


                              [1    0    0    0]
                              [                ]
                              [0    1    0    0]
                     L^(-1)=, [                ]
                              [0    0    1    0]
                              [                ]
                              [0    0    0    1]


                         [1    3    2]
                         [           ]      [1    0    0]
                         [0    0    0]      [           ]
             i=, 3, AA=, [           ], P=, [0    0    1]
                         [0    4    0]      [           ]
                         [           ]      [0    1    0]
                         [0    5    1]


                                   [1    3    2]
                                   [           ]
                                   [0    0    0]
                     i=, 3, j=, 2, [           ]
                                   [0    4    0]
                                   [           ]
                                   [0    0    1]


                          [1    0     0      0]
                          [                   ]
                          [0    1     0      0]
                      L=, [                   ]
                          [0    0     1      0]
                          [                   ]
                          [0    0    -5/4    1]


                             [1    0     0     0]
                             [                  ]
                             [0    1     0     0]
                    L^(-1)=, [                  ]
                             [0    0     1     0]
                             [                  ]
                             [0    0    5/4    1]


                               [1    3    2]
                               [           ]
                               [0    0    0]
                         U2 := [           ]
                               [0    4    0]
                               [           ]
                               [0    0    1]

> print(L2);

                         [1    0     0     0]
                         [                  ]
                         [0    1     0     0]
                         [                  ]
                         [0    0     1     0]
                         [                  ]
                         [0    0    5/4    1]

> print(P);

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

> evalm(L2 &* U2 &* P);

                            [1    2    3]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    4]
                            [           ]
                            [0    1    5]

# 
> printlevel:=3;

                           printlevel := 3

> C := matrix(3,3);

                    C := array(1 .. 3, 1 .. 3, [])

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

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

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

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

> Cinv := myinverse(C);

                   Elem step: exchange rows , 1, 1


                 Elem step: mult row , 1,  by , -1/4


           Elem step: add , 10,  times row , 1,  to row , 2


           Elem step: add , 1,  times row , 1,  to row , 3


                 [1     2    1/4 ]  [-1/4    0    0]
                 [               ]  [              ]
                 [0    22    11/2], [-5/2    1    0]
                 [               ]  [              ]
                 [0    -7    33/4]  [-1/4    0    1]


                   Elem step: exchange rows , 2, 2


                 Elem step: mult row , 2,  by , 1/22


           Elem step: add , -2,  times row , 2,  to row , 1


           Elem step: add , 7,  times row , 2,  to row , 3


                                   [-1      -1      ]
                                   [--      --     0]
                                   [44      11      ]
                 [1    0    -1/4]  [                ]
                 [              ]  [-5              ]
                 [0    1    1/4 ], [--     1/22    0]
                 [              ]  [44              ]
                 [0    0     10 ]  [                ]
                                   [-23             ]
                                   [---    7/22    1]
                                   [22              ]


                   Elem step: exchange rows , 3, 3


                 Elem step: mult row , 3,  by , 1/10


          Elem step: add , 1/4,  times row , 3,  to row , 1


          Elem step: add , -1/4,  times row , 3,  to row , 2


                               [-43     -73         ]
                               [---     ---     1/40]
                               [880     880         ]
                [1    0    0]  [                    ]
                [           ]  [-7               -1 ]
                [0    1    0], [--     3/80      -- ]
                [           ]  [80               40 ]
                [0    0    1]  [                    ]
                               [-23                 ]
                               [---    7/220    1/10]
                               [220                 ]


                            [-43     -73         ]
                            [---     ---     1/40]
                            [880     880         ]
                            [                    ]
                            [-7               -1 ]
                    Cinv := [--     3/80      -- ]
                            [80               40 ]
                            [                    ]
                            [-23                 ]
                            [---    7/220    1/10]
                            [220                 ]

> evalm(C &* Cinv);

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

>