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


void gauss(real matrix a, real matrix b, real scalar det) {
	real scalar i,j,n,s
	real vector js
 
	det = 1
	n = rows(a)
	for (i=1; i<n; i++) {
		maxindex(abs(a[i::n,i]), 1, js=., .)
		j = js[1]+i-1
		if (j!=i) {
			a[(i\j),i..n] = a[(j\i),i..n]
			b[(i\j),.] = b[(j\i),.]
			det = -det
		}
		for (j=i+1; j<=n; j++) {
			s = a[j,i]/a[i,i]
			a[j,i+1..n] = a[j,i+1..n]-s*a[i,i+1..n]
			b[j,.] = b[j,.]-s*b[i,.]
		}
	}
 
	for (i=n; i>=1; i--) {
		for (j=i+1; j<=n; j++) {
			b[i,.] = b[i,.]-a[i,j]*b[j,.]
		}
		b[i,.] = b[i,.]/a[i,i]
		det = det*a[i,i]
	}
}