[pure-lang-svn] SF.net SVN: pure-lang:[851] pure/trunk

[<ag...@us...>]

Revision: 851
          http://pure-lang.svn.sourceforge.net/pure-lang/?rev=851&view=rev
Author:   agraef
Date:     2008-09-25 00:55:11 +0000 (Thu, 25 Sep 2008)

Log Message:
-----------
Add Gaussian elimination example.

Modified Paths:
--------------
    pure/trunk/NEWS

Added Paths:
-----------
    pure/trunk/examples/gauss.pure

Modified: pure/trunk/NEWS
===================================================================
--- pure/trunk/NEWS	2008-09-25 00:54:30 UTC (rev 850)
+++ pure/trunk/NEWS	2008-09-25 00:55:11 UTC (rev 851)
@@ -78,7 +78,7 @@
 
 - #x			total number of elements
 - dim x			number of rows and columns (as a pair)
-- x'			transposed matrix
+- x'			transpose of a matrix
 - x!i			ith matrix element in row-major order (zero-based)
 - x!(i,j)		jth element in ith row of the matrix (zero-based)
 - x!!is, x!!(is,js)	slicing (is and js are lists of machine ints)

Added: pure/trunk/examples/gauss.pure
===================================================================
--- pure/trunk/examples/gauss.pure	                        (rev 0)
+++ pure/trunk/examples/gauss.pure	2008-09-25 00:55:11 UTC (rev 851)
@@ -0,0 +1,69 @@
+
+/* 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.) 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 picot row by the pivot, and subtract suitable multiples of the
+   pivot row to eliminate 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
+    // 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:
+     {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}}
+  else p,i,x
+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<y then j,y else i,x;
+end;
+
+/* Swap rows i and j of the matrix x. This operation is used in the partial
+   pivoting step. (This is not really needed because we could just index the
+   rows indirectly through the row permutation readily available in our
+   implementation of the algorithm, but we omitted this here for clarity.) */
+
+swap x i j = x!!(transp i j (0..n-1),0..m-1) when n,m = dim x end;
+
+/* A little helper function to apply a transposition to a permutation. */
+
+transp i j p = [p!tr k | k=0..#p-1]
+with tr k = if k==i then j else if k==j then i else k end;
+
+/* For convenience, print a double matrix in "short" format a la Octave. */
+
+using system;
+__show__ x::matrix
+= strcat [printd j (x!(i,j))|i=0..n-1; j=0..m-1] + "\n"
+with printd 0 = sprintf "\n%10.5f"; printd _ = sprintf "%10.5f" end
+when n,m = dim x end if dmatrixp x;
+
+/* Example: */
+
+let x = dmatrix {2,1,-1,8; -3,-1,2,-11; -2,1,2,-3};
+x; gauss_elimination x;


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.