```inv ← {vals} ##.gauss_jordan mat            ⍝ Gauss-Jordan elimination.

Gauss-Jordan elimination is a classic algorithm, implemented the D-style.

NB: This  function  is  included  only  for interest as APL provides both matrix
inverse and matrix division as a primitive function: ⌹.

Further, the results of this D-function differ very slightly from those of APL's
primitive  function. This may be because APL's C-code maintains 80-bit floating-
point  accuracy during (much, if not all of) its calculation, whereas the inter-
mediate  arrays of this function use only 64 bits. Of course, the function could
be  adapted  to use, for example, rational →rats← arithmetic for perfect precis-
ion (see test script: ##.scripts.gauss_jordan).

Matrix  right  argument  [mat] represents the coefficients of a system of linear
equations and,  if present, left argument [vals] represents the equations' value
vector (right-hand side).

For example:

x +   y + 2z = ¯3
¯1x + ¯2y + 3z = 14
3x + ¯7y + 4z = ¯3

which may be written:

A x = v
¯ ¯   ¯
where:
┌ 1  1 2┐   v = [¯3 14 ¯3]
A = │¯1 ¯1 3│   ¯
¯   │ 3 ¯7 4│
└       ┘

Then the solution vector x is:
¯
x = v ⌹ A               ⍝ using dyadic ⌹.
¯   ¯   ¯
or
x = (⌹ A) × v           ⍝ using monadic ⌹ (where × is matrix product +.×).
¯      ¯  ¯ ¯                                    ¯
Gauss-Jordan uses three simple transformations, none of which changes the solut-
ion of the equations. For example, starting with:

[0] Initial equations:

x +   y + 2z = ¯3
¯1x + ¯2y + 3z = 14
3x + ¯7y + 4z = ¯3

[1] Exchange rows:

¯1x + ¯2y + 3z = 14     ⍝ swap first and second rows.
x +   y + 2z = ¯3
3x + ¯7y + 4z = ¯3

[2] Multiply each term in a row by the same factor:

¯1x + ¯2y + 3z = 14
3x +  3y + 6z = ¯9     ⍝ multiply second row by 3.
3x + ¯7y + 4z = ¯3

[3] Subtract one row from another:

¯1x +  ¯2y +  3z = ¯2
3x +   3y +  6z = ¯9
0x + ¯10y + ¯2z =  6   ⍝ subtract second row from third.

The algorithm applies these three steps repeatedly, until the coefficient matrix
is transformed into an identity matrix, giving:

1x + 0y + 0z = ¯6
0x + 1y + 0z = ¯1
0x + 0y + 1z =  2

or:

x    ·    · = ¯6
·    y    · = ¯1
·    ·    z =  2

------------
We  find  the explicit matrix inverse by starting with an identity matrix on the
right:

[0]                                         ⍝ initial matrix:
1x +  1y + 2z = 1 0 0
¯1x + ¯2y + 3z = 0 1 0
3x + ¯7y + 4z = 0 0 1
[1]                                         ⍝ exchange first and third rows:
3x + ¯7y + 4z = 0 0 1
¯1x + ¯2y + 3z = 0 1 0
1x +  1y + 2z = 1 0 0
[2]
1x + ¯2.333y + 1.333z = 0 0 0.333      ⍝ multiply first row by ÷3:
¯1x +     ¯2y +     3z = 0 1 0
1x +      1y +     2z = 1 0 0
[3]                                         ⍝ subtract first row from others:
1x + ¯2.333y + 1.333z = 0 0  0.333
0x + ¯4.333y + 4.333z = 0 1  0.333
0x +  3.333y + 0.667z = 1 0 ¯0.333
...                                         ⍝ and so on, until:
1x +  0y + 0z = 0.25 ¯0.346  0.135
0x +  1y + 0z = 0.25 ¯0.038 ¯0.096
0x +  0y + 1z = 0.25  0.192 ¯0.019

Choosing a pivot value
----------------------
In  a  finite-precision implementation, such as IEEE floating point, subtracting
relatively  large  but  commensurate  numbers  introduces  a  significant error.
Gauss-Jordan  reduces this as much as possible by choosing the "pivot" value for
transformation [2] to be the item with the largest absolute magnitude within the
remaining rows of the column under consideration.

before                                                      after
1  0  · · · ·                                               1  0  · · · ·
0  1  · · · ·     In this example, working on the third     0  1  · · · ·
0  0  3 · · ·     column,  the pivot (max abs) value is     0  0 ¯4 · · ·
0  0 ¯2 · · ·     ¯4,  resulting in the exchange of the     0  0 ¯2 · · ·
0  0 ¯4 · · ·     third and fifth rows.                     0  0  3 · · ·
0  0  1 · · ·                                               0  0  1 · · ·

Notice that values in the first two rows of the third column are ignored.

Put simply,  the pivot value for column ⍺ is the item from rows ⍺ downwards with
the largest absolute value: ⌈/|⍺↓⍵[;⍺].

Illustration
------------
Here is an illustration of the steps in the Gauss-Jordan elimination of:

┌────────┐
│4 8 4  0│
│1 4 7  2│
│1 5 4 ¯3│
│1 3 0 ¯2│
└────────┘

Append an identity matrix:

┌────────┬───────┐
│4 8 4  0│1 0 0 0│ append identity matrix.
│1 4 7  2│0 1 0 0│
│1 5 4 ¯3│0 0 1 0│
│1 3 0 ¯2│0 0 0 1│
└────────┴───────┘

Eliminate off-diagonal values from first column:
┌┬───────┬───────┐
│4 8 4  0│1 0 0 0│ first col; pivot value is 4.        (4)· · ·
│1 4 7  2│0 1 0 0│                                      1 · · ·
│1 5 4 ¯3│0 0 1 0│                                      1 · · ·
│1 3 0 ¯2│0 0 0 1│                                      1 · · ·
└┴───────┴───────┘
┌────────┬──────────┐
├1─2─1──0┼0.25─0─0─0┤ divide first row by pivot value.
│1 4 7  2│0    1 0 0│
│1 5 4 ¯3│0    0 1 0│
│1 3 0 ¯2│0    0 0 1│
└────────┴──────────┘
┌─────────┬───────────┐
│1 2  1  0│ 0.25 0 0 0│ subtract first row from remaining rows.
│0 2  6  2│¯0.25 1 0 0│ to leave 0s in all off-diagonal rows
│0 3  3 ¯3│¯0.25 0 1 0│ of first column.
│0 1 ¯1 ¯2│¯0.25 0 0 1│
└─────────┴───────────┘

Eliminate off-diagonal values from second column:
┌──┬──────┬───────────┐
│1 2  1  0│ 0.25 0 0 0│                                 · · · ·
│0 2  6  2│¯0.25 1 0 0│                                 · 2 · ·
│0 3  3 ¯3│¯0.25 0 1 0│ second col; pivot value is 3.   ·(3)· ·
│0 1 ¯1 ¯2│¯0.25 0 0 1│                                 · 1 · ·
└──┴──────┴───────────┘
┌─────────┬───────────┐
│1 2  1  0│ 0.25 0 0 0│                                 · · · ·
├0─3──3─¯3┼¯0.25─0─1─0┤ swap second row with            ·(3)· ·
├0─2──6──2┼¯0.25─1─0─0┤ pivot value row.                · 2 · ·
│0 1 ¯1 ¯2│¯0.25 0 0 1│                                 · 1 · ·
└─────────┴───────────┘
┌─────────┬───────────────────┐
│1 2  1  0│ 0.25    0 0      0│
├0─1──1─¯1┼¯0.08333─0─0.3333─0┤ divide second row by pivot value.
│0 2  6  2│¯0.25    1 0      0│
│0 1 ¯1 ¯2│¯0.25    0 0      1│
└─────────┴───────────────────┘
┌─────────┬────────────────────┐
│1 0 ¯1  2│ 0.4167  0 ¯0.6667 0│ subtract multiples of second row from
│0 1  1 ¯1│¯0.08333 0  0.3333 0│ remaining rows to leave 0s in all
│0 0  4  4│¯0.08333 1 ¯0.6667 0│ off-diagonal rows of second column.
│0 0 ¯2 ¯1│¯0.1667  0 ¯0.3333 1│
└─────────┴────────────────────┘

Eliminate off-diagonal values from third column:
┌─────┬───┬────────────────────┐
│1 0 ¯1  2│ 0.4167  0 ¯0.6667 0│                                · · · ·
│0 1  1 ¯1│¯0.08333 0  0.3333 0│                                · · · ·
│0 0  4  4│¯0.08333 1 ¯0.6667 0│ third col; pivot value is 4.   · ·(4)·
│0 0 ¯2 ¯1│¯0.1667  0 ¯0.3333 1│                                · ·¯2 ·
└─────┴───┴────────────────────┘
┌─────────┬───────────────────────┐
│1 0 ¯1  2│ 0.4167  0    ¯0.6667 0│
│0 1  1 ¯1│¯0.08333 0     0.3333 0│
├0─0──1──1┼¯0.02083─0.25─¯0.1667─0┤ divide third row by pivot value.
│0 0 ¯2 ¯1│¯0.1667  0    ¯0.3333 1│
└─────────┴───────────────────────┘
┌────────┬────────────────────────┐
│1 0 0  3│ 0.3958   0.25 ¯0.8333 0│ subtract multiples of third row from
│0 1 0 ¯2│¯0.0625  ¯0.25  0.5    0│ remaining rows to leave 0s in all
│0 0 1  1│¯0.02083  0.25 ¯0.1667 0│ off-diagonal rows of third column.
│0 0 0  1│¯0.2083   0.5  ¯0.6667 1│
└────────┴────────────────────────┘

Eliminate off-diagonal values from fourth column:
┌───────┬┬────────────────────────┐
│1 0 0  3│ 0.3958   0.25 ¯0.8333 0│                                 · · · ·
│0 1 0 ¯2│¯0.0625  ¯0.25  0.5    0│                                 · · · ·
│0 0 1  1│¯0.02083  0.25 ¯0.1667 0│                                 · · · ·
│0 0 0  1│¯0.2083   0.5  ¯0.6667 1│ fourth col; pivot value is 1.   · · ·(1)
└───────┴┴────────────────────────┘
┌───────┬────────────────────────┐
│1 0 0 0│ 1.021  ¯1.25  1.167  ¯3│ subtract multiples of fourth row
│0 1 0 0│¯0.4792  0.75 ¯0.8333  2│ from remaining rows to leave 0s in all
│0 0 1 0│ 0.1875 ¯0.25  0.5    ¯1│ off-diagonal rows of fourth column.
│0 0 0 1│¯0.2083  0.5  ¯0.6667  1│
└───────┴────────────────────────┘

Drop leading identity matrix for inverse of original matrix:
┌────────────────────────┐
│ 1.021  ¯1.25  1.167  ¯3│
│¯0.4792  0.75 ¯0.8333  2│
│ 0.1875 ¯0.25  0.5    ¯1│
│¯0.2083  0.5  ¯0.6667  1│
└────────────────────────┘

(muse:

As supplied,  the function is offensive to the functional programming purist
because it uses destructive assignment:

mat←⍵                               ⍝ name matrix for updating.
mat[⍺ p;]←mat[p ⍺;]                 ⍝ exchange ⍺th and pth rows.
mat[⍺;]÷←mat[⍺;⍺]                   ⍝ reduce col diagonal to 1.

For fun, we could recast it in a (non-destructive) purer form:

gauss_jordan←{⎕ML ⎕IO←0                   ⍝ Gauss-Jordan elimination.

elim←{                                ⍝ elimination of row/col ⍺.
pivt←{⍵⍳⌈/⍵}|⍵[;⍺]×∨\mask         ⍝ position of pivot row.
mat1←⍺ pivt swap ⍵                ⍝ exchanged ⍺th and pth rows.
mat2←mat1÷[0]⍵[pivt;⍺]*mask       ⍝ ⍺th row divided by pivot.
mat2-(mat2[;⍺]×~mask)∘.×mat2[⍺;]  ⍝ remaining rows reduced by pvt.
}

swap←{                                ⍝ matrix ⍵ with rows ⍺ exchanged.
i←⍳⊃⍴⍵                            ⍝ index vector for rows.
⍵[(n\n/i)+y\⌽y/i;]                ⍝ merge with reversed rows.
}

0::⎕SIGNAL ⎕EN                        ⍝ pass back error to caller.

⍺←=/↑⍳⍴⍵                              ⍝ id matrix for monadic case.

(⍴⍺)⍴(0 1×⍴⍵)↓↑elim/(⌽⍳⌊/⍴⍵),⊂⍵,⍺     ⍝ elimination/ ··· 2 1 0 (⍵,⍺)

}

Further,  we could remove square-bracket-indexing and represent our matrices
as vectors-of-row-vectors:

gauss_jordan←{⎕ML ⎕IO←0                 ⍝ Gauss-Jordan elimination.

elim←{                              ⍝ elimination of row/col ⍺.
pivt←{⍵⍳⌈/⍵}|(⍺⊃¨⍵)×∨\mask      ⍝ position of pivot row.
mat1←⍺ pivt swap ⍵              ⍝ exchanged ⍺th and pivot'th rows.
mat2←mat1÷(pivt ⍺⊃⍵)*mask       ⍝ ⍺th row divided by pivot.
mat2-((⍺⊃¨mat2)×~mask)×⊂⍺⊃mat2  ⍝ remaining rows reduced by pivot.
}

swap←{                              ⍝ swap rows ⍺ in matrix ⍵.
i←⍳⍴⍵                           ⍝ index vector for rows.
((n\n/i)+y\⌽y/i)⊃¨⊂⍵            ⍝ merge with reversed rows.
}

0::⎕SIGNAL ⎕EN                      ⍝ pass back error to caller.

⍺←=/↑⍳⍴⍵                            ⍝ id matrix for monadic case.

(⍴⍺)⍴(0 1×⍴⍵)↓↑↑elim/(⌽⍳⌊/⍴⍵),⊂↓⍵,⍺ ⍝ elimination/ ··· 2 1 0 (⍵,⍺)
}

The  clinically  obsessive would now be in a position to refine the function
still further by replacing local names with additional inner D-functions and
operators.

Pursuing  this  approach  to  the edge of reason, we wind up with a function
that  is a single expression of functions and operators, with no assignments
or  guards.  Here  is  the  code for the monadic case (the dyadic equivalent
replaces  (=/↑⍳⍴⍵)  with  ⍺ in the first line). The expression is strung out
over a number of lines with some white space:

{
·   (=/↑⍳⍴⍵){
·   ·   (⍴⍺)⍴(0 1×⍴⍵)↓↑↑{
·   ·   ·   ⍺(
·   ·   ·   ·   (⍺=⍳⍴⍵){
·   ·   ·   ·   ·   ⍺(
·   ·   ·   ·   ·   ·   ⍺⍺{
·   ·   ·   ·   ·   ·   ·   ⍵-(
·   ·   ·   ·   ·   ·   ·   ·   (⍺⊃¨⍵)×~⍺⍺
·   ·   ·   ·   ·   ·   ·   )×⊂⍺⊃⍵
·   ·   ·   ·   ·   ·   }
·   ·   ·   ·   ·   )⍺(
·   ·   ·   ·   ·   ·   (
·   ·   ·   ·   ·   ·   ·   {⍵⍳⌈/⍵}|(⍺⊃¨⍵)×∨\⍺⍺
·   ·   ·   ·   ·   ·   ){
·   ·   ·   ·   ·   ·   ·   (
·   ·   ·   ·   ·   ·   ·   ·   ⍺ ⍺⍺{
·   ·   ·   ·   ·   ·   ·   ·   ·   ⍺(
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (⍳⍴⍵){
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ↑⍵{
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (⍵\⍵/⍵⍵)+⍺\⌽⍺/⍵⍵
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   )⊃¨⊂⍺⍺
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   }⍺⍺/1 0=⊂⍺⍺∊⍺
·   ·   ·   ·   ·   ·   ·   ·   ·   ·   }
·   ·   ·   ·   ·   ·   ·   ·   ·   )⍵
·   ·   ·   ·   ·   ·   ·   ·   }⍵
·   ·   ·   ·   ·   ·   ·   )÷(⍺⍺ ⍺⊃⍵)*⍵⍵
·   ·   ·   ·   ·   ·   }⍺⍺
·   ·   ·   ·   ·   )⍵
·   ·   ·   ·   }
·   ·   ·   )⍵
·   ·   }/(⌽⍳⌊/⍴⍵),⊂↓⍵,⍺
·   }⍵
}

We  can check that the coding is correct by pasting the lines into character
variable gj_src, say, and then executing:

(⎕io ⎕ml) ⎕pp ← 0 4                             ⍝ condition environment.

(⍎gj_src~⎕tc,'·') 3 3⍴ 1 2 3, ¯3 1 5, 2 4 ¯1    ⍝ apply function.
0.4286 ¯0.2857 ¯0.1429
¯0.1429  0.1429  0.2857
0.2857  0      ¯0.1429
)

Examples:

⎕←mat←4 4⍴ 4 8 4 0, 1 4 7 2, 1 5 4 ¯3, 1 3 0 ¯2
4 8 4  0
1 4 7  2
1 5 4 ¯3
1 3 0 ¯2

⎕pp←4                           ⍝ display numbers to 4 sig-figs.

gauss_jordan mat                ⍝ matrix inverse.
1.021  ¯1.25  1.167  ¯3
¯0.4792  0.75 ¯0.8333  2
0.1875 ¯0.25  0.5    ¯1
¯0.2083  0.5  ¯0.6667  1

1 2 3 4 gauss_jordan mat        ⍝ matrix divide.
¯9.979 6.521 ¯2.813 2.792

(⌹mat)≡gauss_jordan mat         ⍝ check against primitive matrix inverse.
1

(1 2 3 4⌹mat)≡1 2 3 4 gauss_jordan mat      ⍝ ditto matrix divide.
1

⍝ Using an order-⍵ Hilbert matrix, we can see the slight variation between
⍝ primitive ⌹ and the D-function, as errors accumulate.

hil←{÷1+∘.+⍨(⍳⍵)-⎕IO}           ⍝ order-⍵ Hilbert matrix.

⍝ The results are identical up to hil 5:

↑{(⌹⍵)(gauss_jordan ⍵)}∘hil¨ 0 to 5
┌───────────────────────────────────┬───────────────────────────────────┐
│0                                  │0                                  │
├───────────────────────────────────┼───────────────────────────────────┤
│1                                  │1                                  │
├───────────────────────────────────┼───────────────────────────────────┤
│ 4 ¯6                              │ 4 ¯6                              │
│¯6 12                              │¯6 12                              │
├───────────────────────────────────┼───────────────────────────────────┤
│  9  ¯36   30                      │  9  ¯36   30                      │
│¯36  192 ¯180                      │¯36  192 ¯180                      │
│ 30 ¯180  180                      │ 30 ¯180  180                      │
├───────────────────────────────────┼───────────────────────────────────┤
│  16  ¯120   240  ¯140             │  16  ¯120   240  ¯140             │
│¯120  1200 ¯2700  1680             │¯120  1200 ¯2700  1680             │
│ 240 ¯2700  6480 ¯4200             │ 240 ¯2700  6480 ¯4200             │
│¯140  1680 ¯4200  2800             │¯140  1680 ¯4200  2800             │
├───────────────────────────────────┼───────────────────────────────────┤
│   25   ¯300    1050   ¯1400    630│   25   ¯300    1050   ¯1400    630│
│ ¯300   4800  ¯18900   26880 ¯12600│ ¯300   4800  ¯18900   26880 ¯12600│
│ 1050 ¯18900   79380 ¯117600  56700│ 1050 ¯18900   79380 ¯117600  56700│
│¯1400  26880 ¯117600  179200 ¯88200│¯1400  26880 ¯117600  179200 ¯88200│
│  630 ¯12600   56700  ¯88200  44100│  630 ¯12600   56700  ¯88200  44100│
└───────────────────────────────────┴───────────────────────────────────┘

⍝ However, at hil 6, owing to error accumulation,
⍝ small differences begin to appear:

⍪ {(⌹⍵)(gauss_jordan ⍵)}hil 6
┌────────────────────────────────────────────────────────────────────────────────────────────┐
│   36            ¯630.0000001     3360          ¯7560.000001     7560.000001    ¯2772.000001│
│ ¯630.0000001   14700           ¯88200.00001   211680         ¯220500           83160.00002 │
│ 3360          ¯88200.00001     564480.0001  ¯1411200         1512000         ¯582120.0001  │
│¯7560.000001   211680         ¯1411200        3628800.001    ¯3969000.001     1552320       │
│ 7560.000001  ¯220500          1512000       ¯3969000.001     4410000.001    ¯1746360       │
│¯2772.000001    83160.00002    ¯582120.0001   1552320        ¯1746360          698544.0001  │
├────────────────────────────────────────────────────────────────────────────────────────────┤
│   36           ¯630           3360          ¯7560     7560.000001    ¯2772                 │
│ ¯630          14700         ¯88200.00001   211680  ¯220500           83160.00001           │
│ 3360         ¯88200.00001   564480       ¯1411200  1512000         ¯582120                 │
│¯7560         211680       ¯1411200        3628800 ¯3969000         1552320                 │
│ 7560.000001 ¯220500        1512000       ¯3969000  4410000        ¯1746360                 │
│¯2772          83160.00001  ¯582120        1552320 ¯1746360          698544.0001            │
└────────────────────────────────────────────────────────────────────────────────────────────┘