Task

Invert matrix   A   using Gauss-Jordan method.

A   being an   n × n   matrix.




V Eps = 1e-10
 
F transformToRref(&mat)
   V lead = 0
 
   L(r) 0 .< mat.len
      I lead >= mat[0].len
         R
 
      V i = r
      L mat[i][lead] == 0
         i++
         I i == mat.len
            i = r
            lead++
            I lead == mat[0].len
               R
      swap(&mat[i], &mat[r])
 
      V d = mat[r][lead]
      I abs(d) > Eps
         L(&item) mat[r]
            item /= d
 
      L(i) 0 .< mat.len
         I i != r
            V m = mat[i][lead]
            L(c) 0 .< mat[0].len
               mat[i][c] -= mat[r][c] * m
 
      lead++
 
F inverse(mat)
   V augmat = [[0.0] * (2 * mat.len)] * mat.len
   L(i) 0 .< mat.len
      L(j) 0 .< mat.len
         augmat[i][j] = mat[i][j]
      augmat[i][mat.len + i] = 1
 
   transformToRref(&augmat)
 
   V result = [[0.0] * mat.len] * mat.len
   L(i) 0 .< mat.len
      L(j) 0 .< mat.len
         I augmat[i][j] != Float(i == j)
            X ValueError(‘matrix is singular’)
         result[i][j] = augmat[i][mat.len + j]
   R result
 
F print_mat(mat)
   L(row) mat
      V line = ‘’
      L(val) row
         I !line.empty
            line ‘’= ‘ ’
         line ‘’= ‘#3.5’.format(val)
      print(line)
 
F runTest(mat)
   print(‘Matrix:’)
   print_mat(mat)
   print()
   print(‘Inverse:’)
   print_mat(inverse(mat))
   print()
   print()
 
V m1 = [[Float(1), 2, 3],
        [Float(4), 1, 6],
        [Float(7), 8, 9]]
 
V m2 = [[Float( 2), -1,  0],
        [Float(-1),  2, -1],
        [Float( 0), -1,  2]]
 
V m3 = [[Float(-1), -2,  3,  2],
        [Float(-4), -1,  6,  2],
        [Float( 7), -8,  9,  1],
        [Float( 1), -2,  1,  3]]
 
runTest(m1)
runTest(m2)
runTest(m3)