http://www.math.hut.fi/teaching/y3/opemate/G-J.txt
----------------------------------------------------------------------------

> # Gauss-Jordan, GREE esim. kts. s. 13 yleinen selostus G-J:sta
> with(linalg):
--------------------------------------------------------------------------------
> A:=matrix([[1,-1,1,0,0,2],[1,1,2,4,1,2],[0,2,1,4,3,1],[1,-3,0,-4,-2,1]]);

                               [ 1  -1  1   0   0  2 ]
                               [                     ]
                               [ 1   1  2   4   1  2 ]
                          A := [                     ]
                               [ 0   2  1   4   3  1 ]
                               [                     ]
                               [ 1  -3  0  -4  -2  1 ]
--------------------------------------------------------------------------------
> b:=vector([0,3,2,0]);

                               b := [ 0, 3, 2, 0 ]
--------------------------------------------------------------------------------
> Ab:=augment(A,b);

                              [ 1  -1  1   0   0  2  0 ]
                              [                        ]
                              [ 1   1  2   4   1  2  3 ]
                        Ab := [                        ]
                              [ 0   2  1   4   3  1  2 ]
                              [                        ]
                              [ 1  -3  0  -4  -2  1  0 ]
--------------------------------------------------------------------------------
> # 1. sarake:
> Ab:=addrow(Ab,1,2,-1);

                              [ 1  -1  1   0   0  2  0 ]
                              [                        ]
                              [ 0   2  1   4   1  0  3 ]
                        Ab := [                        ]
                              [ 0   2  1   4   3  1  2 ]
                              [                        ]
                              [ 1  -3  0  -4  -2  1  0 ]
--------------------------------------------------------------------------------
> Ab:=addrow(Ab,1,4,-1);

                             [ 1  -1   1   0   0   2  0 ]
                             [                          ]
                             [ 0   2   1   4   1   0  3 ]
                       Ab := [                          ]
                             [ 0   2   1   4   3   1  2 ]
                             [                          ]
                             [ 0  -2  -1  -4  -2  -1  0 ]
--------------------------------------------------------------------------------
> #2. sarake:
> Ab:=addrow(Ab,2,3,-1);

                             [ 1  -1   1   0   0   2   0 ]
                             [                           ]
                             [ 0   2   1   4   1   0   3 ]
                       Ab := [                           ]
                             [ 0   0   0   0   2   1  -1 ]
                             [                           ]
                             [ 0  -2  -1  -4  -2  -1   0 ]
--------------------------------------------------------------------------------
> Ab:=addrow(Ab,2,4,1);

                              [ 1  -1  1  0   0   2   0 ]
                              [                         ]
                              [ 0   2  1  4   1   0   3 ]
                        Ab := [                         ]
                              [ 0   0  0  0   2   1  -1 ]
                              [                         ]
                              [ 0   0  0  0  -1  -1   3 ]
--------------------------------------------------------------------------------
> # 5. sarake:
> Ab:=addrow(Ab,3,4,1/2);

                             [ 1  -1  1  0  0    2    0  ]
                             [                           ]
                             [ 0   2  1  4  1    0    3  ]
                       Ab := [                           ]
                             [ 0   0  0  0  2    1    -1 ]
                             [                           ]
                             [ 0   0  0  0  0  -1/2  5/2 ]
--------------------------------------------------------------------------------
> # Takaisinsijoitus:
> # Normeeraus:
> Ab:=mulrow(Ab,4,-2);

                               [ 1  -1  1  0  0  2   0 ]
                               [                       ]
                               [ 0   2  1  4  1  0   3 ]
                         Ab := [                       ]
                               [ 0   0  0  0  2  1  -1 ]
                               [                       ]
                               [ 0   0  0  0  0  1  -5 ]
--------------------------------------------------------------------------------
> Ab:=mulrow(Ab,3,1/2);

                             [ 1  -1  1  0  0   2     0  ]
                             [                           ]
                             [ 0   2  1  4  1   0     3  ]
                       Ab := [                           ]
                             [ 0   0  0  0  1  1/2  -1/2 ]
                             [                           ]
                             [ 0   0  0  0  0   1    -5  ]
--------------------------------------------------------------------------------
> Ab:=mulrow(Ab,2,1/2);

                           [ 1  -1   1   0   0    2     0  ]
                           [                               ]
                           [ 0   1  1/2  2  1/2   0    3/2 ]
                     Ab := [                               ]
                           [ 0   0   0   0   1   1/2  -1/2 ]
                           [                               ]
                           [ 0   0   0   0   0    1    -5  ]
--------------------------------------------------------------------------------
> #6. sarake:
> Ab:=addrow(Ab,4,3,-1/2);

                            [ 1  -1   1   0   0   2   0  ]
                            [                            ]
                            [ 0   1  1/2  2  1/2  0  3/2 ]
                      Ab := [                            ]
                            [ 0   0   0   0   1   0   2  ]
                            [                            ]
                            [ 0   0   0   0   0   1   -5 ]
--------------------------------------------------------------------------------
> Ab:=addrow(Ab,4,1,-2);

                            [ 1  -1   1   0   0   0   10 ]
                            [                            ]
                            [ 0   1  1/2  2  1/2  0  3/2 ]
                      Ab := [                            ]
                            [ 0   0   0   0   1   0   2  ]
                            [                            ]
                            [ 0   0   0   0   0   1   -5 ]
--------------------------------------------------------------------------------
> # 5. sarake:
> Ab:=addrow(Ab,3,2,-1/2);

                             [ 1  -1   1   0  0  0   10 ]
                             [                          ]
                             [ 0   1  1/2  2  0  0  1/2 ]
                       Ab := [                          ]
                             [ 0   0   0   0  1  0   2  ]
                             [                          ]
                             [ 0   0   0   0  0  1   -5 ]
--------------------------------------------------------------------------------
> # 2. sarake
> Ab:=addrow(Ab,2,1,1);

                             [ 1  0  3/2  2  0  0  21/2 ]
                             [                          ]
                             [ 0  1  1/2  2  0  0   1/2 ]
                       Ab := [                          ]
                             [ 0  0   0   0  1  0    2  ]
                             [                          ]
                             [ 0  0   0   0  0  1   -5  ]
--------------------------------------------------------------------------------
> # Nyt on rref eli redusoitu porrasmuoto valmis. Tarkistetaan viela
> 
> rref(augment(A,b));    ## rref:n synonyymi:  gaussjord
                         ## joissakin mikroluokan koneissa rref tökkii,
                         ## käytä silloin gaussjord:ia

                          [ 1  0  3/2  2  0  0  21/2 ]
                          [                          ]
                          [ 0  1  1/2  2  0  0   1/2 ]
                          [                          ]
                          [ 0  0   0   0  1  0    2  ]
                          [                          ]
                          [ 0  0   0   0  0  1   -5  ]
--------------------------------------------------------------------------------
> # kylla vaan!!


######################################################################
# Lopuksi permutoidaan sarakkeet niin, että "johtavat sarakkeet"
# tulevat ensimmäisiksi, jolloin saadaan rakennne, jossa on
#           [I(4x4),P(4x2),b(4x1)]
# eli:
#
# [4x4-yksikkömatriisi, 2 "parametrisaraketta", oikea puoli (b-sarake)]
--------------------------------------------------------------------------------
> submatrix(Ab,[1,2,3,4],[1,2,5,6,3,4,7]);

                          [ 1  0  0  0  3/2  2  21/2 ]
                          [                          ]
                          [ 0  1  0  0  1/2  2   1/2 ]
                          [                          ]
                          [ 0  0  1  0   0   0    2  ]
                          [                          ]
                          [ 0  0  0  1   0   0   -5  ]
--------------------------------------------------------------------------------
>