/* Pure matrix example: Gaussian elimination algorithm. 2008-09-25 AG */ /* The Gaussian elimination algorithm brings a matrix into row echelon form, which makes it possible to solve the original system of linear equations using back substitution. Our version of the algorithm just returns the row echelon form of the matrix, along with the corresponding permutation of the row indices performed during pivoting. */ gauss_elimination x::matrix = p,x when n,m = dim x; p,_,x = foldl step (0..n-1,0,x) (0..m-1) end; /* The elimination step. x is our matrix, i the current row index, j the current column index, and p keeps track of the current permutation of the row indices performed during pivoting. The algorithm returns the updated matrix x, row index i and row permutation p. */ /* Here is a brief rundown of the algorithm: First we find the pivot element in column j of the matrix. (We're doing partial pivoting here, i.e., we only look for the pivot in column j, starting at row i. That's usually good enough to achieve numerical stability.) If the pivot is zero then we're done (the pivot column is already zeroed out). Otherwise, we bring it into the pivot position (swapping row i and the pivot row), divide the pivot row by the pivot, and subtract suitable multiples of the pivot row to eliminate the elements of the pivot column in all subsequent rows. Finally we update i and p accordingly and return the result. */ step (p,i,x) j = if max_x==0 then p,i,x else // updated row permutation and index: transp i max_i p, i+1, {// the top rows of the matrix remain unchanged: x!!(0..i-1,0..m-1); // the pivot row, divided by the pivot element: {x!(i,l)/x!(i,j) | l=0..m-1}; // subtract suitable multiples of the pivot row: {x!(k,l)-x!(k,j)*x!(i,l)/x!(i,j) | k=i+1..n-1; l=0..m-1}} when n,m = dim x; max_i, max_x = pivot i (col x j); x = if max_x>0 then swap x i max_i else x; end with pivot i x = foldl max (0,0) [j,abs (x!j)|j=i..#x-1]; max (i,x) (j,y) = if x