Math 236  •  26 February 2001

Elementary row operations:
swap, scale, addrow (and roundoff)

You may use this notebook in any of several ways:
• as a source for the functions swap, scale, addrow, and roundoff used to carry out elementary row operations;
• as solutions to Homework Set 2, Problem 1 (except that the requisite testing is not shown here); and
• as a tutorial on several aspects to programming in Mathematica.

Feel free to use any of the definitions of swap and scale here in other homework problems in place of your own functions and in place of the ones available from the rowops package.

Definitions of swap and scale

Here are definitions of swap and scale along with usage messages for them:

Clear[swap,  scale]

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

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

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

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

Explanation of definitions

Here's the explanation of what happened in swap in the expression A[[{i, k}]] = A[[{k, i}]].  For any list v—including a list of lists such as a matrix—an expression of the form v[[{i, k}]] produces the list consisting of the ith and kth members of v, in that order.  Then an assignment of the form v[[{i, k}]] = v[[{k, i}]] interchanges the ith and kth members of v.

If you did not understand how to write swap the way shown above, you could introduce an extra local variable temp, say, that holds one of the rows being swapped while you put the other one in its place:

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

The general form of a Module expression is 
	Module[{vars}, body]
where the first argument {vars} is the list of temporary (local) variables and the second argument body is an expression that could be a compound expression.  Indeed, in both swap and scale the expression body is a compound expression formed from two expressions separated by a semicolon.

As written, swap and scale initialize the local variable A within the first argument {A = mat}.  You could, instead, do that initialization separately, as the first statement in the compound expression that forms the body of the Module.  For example:
	swap[mat_?MatrixQ, i_Integer, k_Integer] := Module[{A},
		A = mat; A[[{i,k}]] = A[[{k,i}]]; A]
However, it is shorter and probably clearer to initialize A within the temporary variable list, as shown in the original definitions.

In scale, you do not need any multiplication symbol between c and A[[i]] in the expression c A[[i]].

The function Return is not needed to return the value

In both swap and scale, the second argument of Module—the body of the Module—is a compound expression consisting of two expressions separated by a semicolon ( ; ).  For example, the body of the scale Module is A[[i]]=c A[[i]]; A.  The result of an evaluated Module expression is always the value of the final expression in the body, which here is the (new) value of A.

Accordingly, there is absolutely no need to use Return[A] as the second of the two expressions in the body of the Modules.  The only situation when Return is used is when you need to "jump out" of something prematurely, for example, an If expression testing for an error that is nested deeply in an expression.

⚠  Wrong to return result in MatrixForm

It would be dead wrong to write your scale, say, in a form such as the following:

wrongScale[mat_?MatrixQ,  c_?NumericQ,  i_Integer] := \(Module[{A\  = \ mat}, \ 
\ \ \ \ 	A[[
              i]]\  = \ 
            c\ A[[i]]; \ 
          MatrixForm[A]]\)

The reason that it would be wrong is that if A is a matrix—a particular form of a list of lists—then MatrixForm[A] is not a matrix.  Rather, it is a special display form of a matrix.  So if you returned such a display form as the result of your function, then you would be unable to operate further upon it.  Try the following, for example:

Clear[A]
A  =  { {1,  2,  3},  {4,  5,  
            6} };
wrongScale[A,  10,  1]
5 wrongScale[A,  10,  1]

Type-checking of input

The input type-checking represented by ?MatrixQ and Integer in the pattern names mat_?MatrixQ, i_Integer, and k_Integer is useful but not essential.  The patterns i_Integer and k_Integer match any expression whose Head is Integer; in other words, they match any expression whose value is an integer.  The pattern mat_?MatrixQ matches any expression whose value is a matrix—a list of lists where each of those inner lists includes no deeper-level list and where all of those inner lists have the same length.

Fancier type-checking of input

The checking that i and k are integers is really insufficient if you want to "idiot-proof" the functions swap and scale.  What you really want to check is that these integers are among the row numbers of the matrix.  To check that, the easiest way is through use of Condition (abbreviation /; ), as will be done below.

For convenience, first create a utility function rowIndexQ to check whether an integer is a legal row index of a matrix:

Clear[rowIndexQ]
rowIndexQ[mat_?MatrixQ,  j_Integer]  :=  
      MemberQ[Range[Length[mat]],  j]

For example:

Clear[A]
A  =  
      Partition[Table[Random[Integer,  {-10,  50}],  {12}],  4];  
    A  //  MatrixForm

rowIndexQ[A,  2]
rowIndexQ[A,  5]

Finally, incorporate into the definition of swap the test for whether the arguments i and k are legitimate row indices of the matrix by means of Condition, using the abbreviation (that is, special input form) /;:

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

swap[A,  1,  2]  //  MatrixForm

swap[A,  1,  5]

The Condition (expressed through the /; form) just used in the revised definition of swap can also be placed at the end of the right side of the entire epression:

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

Exercise .   Use

Exercise.  Similarly use Condition ( /; ) in the defintion of scale so as to check that the integer argument is a legitimate index.

⚠  Why ReplaceAll ( /. ) with Rule ( -> ) do not  work correctly

It is very tempting to try something like the following method for defining swap:

badScale[mat_?MatrixQ,   c_?NumericQ,  
        i_Integer]  :=   mat  /.   mat[[
            i]]  ->   c mat[[
              i]]

And in some—but not all—situations this will seem to work.  For example:

Clear[A]
A  =  Partition[Range[12],  4];  A  //  MatrixForm

badScale[A,  10,  1]

But consider the following example:

A  =  { {1,  2,  3,  4},  {5,  6,  7,  8},  {1,  2,  3,  
          4} };  A  //  MatrixForm

badScale[A,  10,  1]

What went wrong?  Nothing, really: Mathematica did exactly what it was asked to do!  The trouble is in what it means to use ReplaceAll ( /. ) with a Rule ( -> ).  ReplaceAll tries to use the rule mat[[i]] -> c mat[[i]] everywhere within the left argument mat to /..  It looks for the pattern math[[i]] everywhere within the left argument mat.  And in the preceding example, since for i = 1 the pattern mat[[i]] is the constant vector {1, 2, 3, 4}, this pattern matches both the first and the third rows of mat.  Accordingly, both the first and third rows of mat get scaled by 10.

Exercise.  Explain in detail why the following definition for swap would be wrong.  
	badSwap[mat_, i_, k_] := mat /. {mat[[i]] ->  mat[[k]], mat[[k]] -> mat[[i]]}
Demonstrate that badSwap gives an unwanted result for a suitable matrix.

Definition of addrow

The function addrow, along with a usage message for it, is given below.

Clear[addrow]
addrow::usage  =  "addrow[mat, c, k, i] replaces the ith row of \
matrix mat by the sum of that row and scalar c times the kth row.";
addrow[mat_List,  c_?NumericQ,  k_Integer,  i_Integer]  :=  
      Module[{A  =  mat},  
        A[[i]]  =  A[[i]]  +  c A[[k]]; A]

Exercise.  Include further input type-checking in this version of addrow.

Definition of roundoff

Clear[roundoff]
roundoff::usage  =  "roundoff[mat] replaces entries of matrix mat \
that are less than 10^(-12) in absolute value with 0 exactly.";
roundoff[a_?NumericQ]  :=  
      If[ Abs[a]  <   1.*^-12,  0,  a]
SetAttributes[roundoff,  Listable]

Why make roundoff Listable?

The reason for the SetAttributes expression above is that, without it, roundoff would balk at evaluating with a list as argument.  [The difficulty is that, whereas Abs is Listable, neither Less ( < ) nor If is Listable.]  To see that this is so, try the following:

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

scalarRoundoff[{10. ^-15,  100. }]

Another way around this difficulty, without explicitly giving roundoff the Listable attribute, would be to use Map each time that you needed to roundoff a list of numbers:

Map[scalarRoundoff,  {10. ^-15,  100. }]

But such a use of Map would not directly work if the argument were a list of lists—in particular, if it were a matrix:

Map[
      scalarRoundoff,  { {1,  2,  3},  {10. ^-14,  10. ^-15,  
          10. ^-16} }]

In such a case, you could use Map with a third argument, specifying the level to which scalarRoundoff is to be mapped.  For a matrix, or any list of lists, level {2} indicates the individual entries of the bottm-level lists.

Map[
      scalarRoundoff,  { {1,  2,  3},  {10. ^-14,  10. ^-15,  
          10. ^-16} },  {2}]

You can see now why it's so much simpler here to take care of all numerical lists, no matter how deeply nested, by using SetAttributes to give roundoff the attribute Listable.