Math 236  •  Spring 2001

About function GJ

Version 1.2:  6 March 2001

My function gaussJordan

This is the function you must modify to get your GJ.

⚠  Before you use gaussGJ as is or as modified, you must have defined scale and swap (or used mine) as well as evaluated all the cells in the section Auxiliary functions, below.

gaussJordan[mat_] := 
Module[
	{a = N[mat], m, n, r, c, p},
	m = Length[a]; n = Last @ Dimensions @ a;
	Do[ (* for each row...let c = index of next pivot column *)
		c = firstNonzeroColumn[Drop[a, r - 1]];
		If[c <= n,         (* then...get index of pivot's row *)
			p = r - 1 + firstNonzeroIndex[Drop[a[[All, c]], r - 1]];
			a = swap[a, r, p];      (* put pivot into position *)       
			a = scale[a, 1./a[[r, c]], r];      (* get leading 1 *)
			a = zeroUpDown[a, r, c]		(* 0s in rest of column *)
			, (* else *)  Null
			]                                        (* end If *)
		, {r, 1, m }
	];                                             (* end Do *)
 a
]

Auxiliary functions (needed for gaussJordan)

General-purpose functions for lists

Usage messages are defined at the end of this group of functions definitions.  To find out the usage of one of these functions, use ? (or ??) in the usual way, for example:  ?firstIndex

Usage messages:

someNonzeroQ::usage  =  "someNonzeroQ[vec] returns True or False \
according to whether the list vec has some nonzero entry or not.";

firstIndex::usage  =  "firstIndex[vec, b] returns the index of \
the first entry in the list vec for which the corresponding entry in the \
Boolean list b is True; if all entries in b are False, the result is \
Length[vec] + 1.";

firstNonzeroIndex::usage  =  "firstNonzeroIndex[vec] returns the \
index of the first nonzero entry in the list vec.";

nonzeroQ::usage  =  "nonzeroQ[x] returns True or False according \
to whether x is not 0 or is 0; nonzeroQ[vec] applies nonzeroQ to each entry \
of the list vec.";

Functions:

nonzeroQ[x_]  :=  N[x]  ≠  0. ;  
    SetAttributes[nonzeroQ,  Listable]

someNonzeroQ[v_List]  :=  Apply[Or,  nonzeroQ[v]]

firstIndex[v_List,  b_List  /;  Not[Apply[Or,  b]]]  :=  
      Length[v]  +  1
firstIndex[v_List,  b_List]  :=  
      First[Flatten[Position[b,  True]]]

firstNonzeroIndex[v_List]  :=  
      firstIndex[v,  nonzeroQ[v]]

Special-purpose functions for manipulating matrices

Usage messages are defined at the end of this group of functions definitions.  To find out the usage of one of these functions, use ? (or ??) in the usual way, for example:  ?firstNonzeroColumn

Usage messages:

nonzeroRowsQ::usage  =  "nonzeroRowsQ[A] returns a Boolean list \
with entry True corresponding to each nonzero row, and entry False for each \
zero row, of the matrix A.";
nonzeroColumnsQ::usage  =  "nonzeroColumnsQ[A] returns a Boolean \
list with entry True corresponding to each nonzero column, and entry False \
for each zero column, of the matrix A.";

firstNonzeroColumn::usage  =  "firstNonzeroColumn[A] returns the \
index of the first nonzero column of the matrix A.";

swap::usage  =  "swap[mat, i, k] swaps the ith and kth rows of \
matrix mat.";

scale::usage  =  "scale[mat, c, i] multiplies the ith row of \
matrix mat by scalar c.";

zeroUpDown::usage  =  "zeroUpDown[A, i, j] applies addrow \
elementary row operations to the matrix A to change to 0 all entries in the \
jth column of A except the ith; A[[i, \
j]] must be 1.";

roundoff::usage  =  "roundoff[mat] replaces entries of matrix mat \
that are less than 10^(-12) in absolute value with 0 exactly.";

Functions:

nonzeroRowsQ[mat_]  :=  Map[someNonzeroQ,  mat]
nonzeroColumnsQ[mat_]  :=  nonzeroRowsQ[Transpose[mat]]

firstNonzeroColumn[mat_]  :=  
      firstIndex[
        Range[Dimensions[
              mat][[2]]],   nonzeroColumnsQ[mat]]

swap[mat_?MatrixQ,  i_Integer,  k_Integer]  :=  
      Module[{A  =  mat},  
        A[[{i, k}]]  =  
          A[[{k, i}]];  A]

scale[mat_?MatrixQ,  c_?NumericQ,  i_Integer]  :=  
      Module[{A  =  mat},  
        A[[i]]  =  
          c A[[i]];  A]

zeroUpDown[mat_,  i_,  
        j_] := Module[{a  =  mat,  is,  
          cs},  
is  =  
          Complement[Range[Length[a]],  {i}];  
        cs  =   a[[is,  j]] ;  
        a[[
            is]]  =   a[[
              is]]  -  
            Outer[Times,  cs,  
              a[[i]]];  
        a ]

roundoff[a_?NumericQ]  :=  
      If[ Abs[a]  <   1.*^-12,  0,  a]
SetAttributes[roundoff,  Listable]

Example

Clear[a]
a  =  {{0, 0, 1, -1, -1, 4}, {2, 4, 2, 4, 2, 4}, {2, 4, 3, 
            3, 3,  4}, {3, 6, 6, 3, 6, 6}};
a  //  MatrixForm
gaussJordan[a]   //  MatrixForm
RowReduce[a]    //  MatrixForm       
gaussJordan[a]  ==  RowReduce[a]

Define your functions

In a separate notebook—not this one—define your putOne and pivotLocate.  Copy my gaussJordan from here to there and edit it suitably to convert it into upir function GJ.

Get the validator package vdGJ.m

At any time before you are ready to validate your function GJ, download my encoded package vdGJ.m in any one of the following three ways:

•  At UMASTR lab only: from the Math_236 directory on drive M: (labelled Math-Stat Files).  An icon for this drive is on the desktop of UMASTR PCs.

• Obtain the file vdGJ.zip from the Files page of the course web site http://www.math.umass.edu/Courses/Math_236.  Then extract the single file vdGJ.m from vdGJ.zip by using WinZip or another "unzip" program.  Be sure that the unzip program is configured so that it does not not change end-of-line characters!

•  If you have an ftp client that does not change end-of-line characters (so not Netscape or Internet Explorer!), get the file vdGJ.m via anonymous ftp from:
	umastr.math.umass.edu
Set the initial remote directory to
 	/pub/lab/Math_236
and use binary mode (not text mode) to download the file vdGJ.m.  For network security reasons, such anonymous ftp access is limited to client compuers that are in the umass.edu domain.   This means that the ftp access must be made directly through a campus network-connected computer or via dial-up connection to OIT PPP ("UMAccess").

If you are using Mathematica on a Mac or under Linux

The latest information I have is that you need to change the end-of-line characters from the DOS convention of CR LF to the convention on your platform (just CR for Mac; just LF for Linux).  For example, under Linux, if you have downloaded vdGJ.zip as a binary file, give the command:
	unzip -a vdGJ.zip
This will extract vdGJ.m and, as it does so, convert the end-of-lines to the *nix LF.

Validate your function GJ

A program of mine is going to generate a set of test matrices and evaluate your function GJ in turn at each of these matrices.  The test matrices it uses for you are not necessarily the same as the test matrices it uses for others.  In fact, some of these matrices involve random numbers, and these random numbers depend in part on the digits in your student ID number.

Do the following four steps:

1.  Be sure that at this Mathematica session you have defined your GJ and all the other objects it references directly or indirectly—your putOne and pivotLocate and any other functions you decided to use; my functions from the section Auxiliary functions, above; your swap and scale or mine (also defined in Auxiliary functions).

2.  Inside the double quotes in the next cell, type your  student ID number and then evaluate the cell.

myID = "";           
(* type your student ID number inside the quotes *)

Check that you correctly defined myID—evaluate the following cell now.  You should see your ID number as result.

myID

3.  Tell Mathematica where to look for my validator package vdGJ.m—edit the following cell so that its first line shows the correct path to the directory to which you downloaded or extracted vdGJ.m:

\(vdGJpath\  = \ "c:/data"; \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      path\ to\ vdGJ . m\ *) 
AppendTo[$Path, \ vdGJpath];\)

4.  Evaluate the preceding cell.  Check that the preceding cell was indeed evaluated—evaluate the following cell.  The result should include at its end the path to the directory holding the file vdGJ.m.

5.  Run the validator  that tests your function—evaluate the following cell now:

Needs["vdGJ`"]            (* run the validator *)

If you have forgotten to define myID correctly, then you can simply edit myID to make it correct—be sure your ID number is enclosed in double quotes—and then run the preceding cell with Needs["vdGJ`"] again.

If, however, at the current Mathematica session you forgot to define GJ or some auxiliary function needed by GJ, or if you need to fix up your GJ to get it to work properly, then you should:
•  save your notebook containing your function definitions;
•  save this notebook AboutGJ.nb with your myID and vdGJpath defined;
•  quit Mathematica and start it again (or, if you know how, leave Mathematica open but quit the kernel and restart the kernel);
•  open your notebook with your function definitions and evaluate them all;
•  open this notebook AboutGJ.nb, evaluate its cells defining my auxiliary functions and your myID and vdGJpath; and
•  run the validator again by evaluating the cell Needs["vdGJ`"] .

What to turn in

Do not turn in a hardcopy of this entire notebook!  Rather, turn in only:

•  A listing of all the objects you defined for GJ, including of course your GJ itself and your putOne and pivotLocate.

•  The validation output that results from evaluating the cell Needs["vdGJ`"].