Task

Solve   Ax=b   using Gaussian elimination then backwards substitution.

A   being an   n by n   matrix.

Also,   x and b   are   n by 1   vectors.

To improve accuracy, please use partial pivoting and scaling.



See also

   the Wikipedia entry:   Gaussian elimination


F swap_row(&a, &b, r1, r2)
   I r1 != r2
      swap(&a[r1], &a[r2])
      swap(&b[r1], &b[r2])
 
F gauss_eliminate(&a, &b)
   L(dia) 0 .< a.len
      V (max_row, max) = (dia, a[dia][dia])
      L(row) dia+1 .< a.len
         V tmp = abs(a[row][dia])
         I tmp > max
            (max_row, max) = (row, tmp)
 
      swap_row(&a, &b, dia, max_row)
 
      L(row) dia+1 .< a.len
         V tmp = a[row][dia] / a[dia][dia]
         L(col) dia+1 .< a.len
            a[row][col] -= tmp * a[dia][col]
         a[row][dia] = 0
         b[row] -= tmp * b[dia]
 
   V r = [0.0] * a.len
   L(row) (a.len-1 .. 0).step(-1)
      V tmp = b[row]
      L(j) (a.len-1 .< row).step(-1)
         tmp -= r[j] * a[row][j]
      r[row] = tmp / a[row][row]
   R r
 
V a = [[1.00, 0.00, 0.00,  0.00,  0.00, 0.00],
       [1.00, 0.63, 0.39,  0.25,  0.16, 0.10],
       [1.00, 1.26, 1.58,  1.98,  2.49, 3.13],
       [1.00, 1.88, 3.55,  6.70, 12.62, 23.80],
       [1.00, 2.51, 6.32, 15.88, 39.90, 100.28],
       [1.00, 3.14, 9.87, 31.01, 97.41, 306.02]]
V b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]
 
print(gauss_eliminate(&a, &b))