Implementing Gaussian Elimination

We implement Gaussian Elimination in three stages of increasing sophistication,
following the development in Olver and Shakiban, "Applied Linear Algebra," Chapter 1.

The simplest case occurs when we always find the pivot positions to be occupied by non-zero elements. Following Olver and Shakiban, we call this the "regular case."
The next simplest case occurs when there will at least be a suitable nonzero pivot available somewhere, but we may have to swap rows in order to move it into position.
This occurs when the matrix is nonsingular.
In the third and final case, we institute a pivoting policy to improve the numerical performance of the algorithm.

NOTE: In this preliminary version, the algorithms are implemented for square matrices only

Gaussian Elimination - Regular Case

The "Regular Case" refers to a square matrix which can be row reduced to an upper triangular matrix by the following algorithm.
The critical point is that there must be a non-zero element in each pivot position for the algorithm to be successful.
We do not resort to swapping rows. A more general case is handled later.
See Olver and Shakiban, "Applied Linear Algebra," Section 1.3.

gaussianEliminationRegularCase

In[1]:=

gaussianEliminationRegularCase :: Usage = "gaussianEliminationRegularCase[a_?squareMatr ... cesary to swap rows while executing the Gaussian Elimination algorithm on such a matrix." ;

In[2]:=

gaussianEliminationRegularCase[a_ ? squareMatrixQ] := Module[{n = Dimensions[a][[1]] ... ; b[[i]] = b[[i]] - c b[[j]],  {i, j + 1, n}]],  {j, 1, n}] ; b]

squareMatrixQ

In[3]:=

squareMatrixQ :: Usage = "squareMatrixQ[a_?MatrixQ] returns True if a is a square matrix, and otherwise returns False." ;

In[4]:=

squareMatrixQ[a_ ? MatrixQ] := Module[{m = Dimensions[a][[1]], n = Dimensions[a][[2]]}, If[mn, True, False]]

Testing.

In[5]:=

a = ({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}) ; 

gaussianEliminationRegularCase[a] ;

%//MatrixForm

Out[7]//MatrixForm=

( {{1, 2, 3}, {0, -3, -6}, {0, 0, 0}} )

In[8]:=

a = ({{0, 2, 3}, {4, 5, 6}, {7, 8, 9}}) ; 

gaussianEliminationRegularCase[a]

Error: gaussianEliminationRegularCase: A is not regular.

Out[9]=

$Aborted

Gaussian Elimination - Nonsingular Case

The following procedure will row reduce a nonsingular square matrix to an upper triangular matrix.
The naive pivoting policy is merely to move the nearest nonzero candidate into the pivot position.
See Olver and Shakiban, "Applied Linear Algebra," Section 1.4.

gaussianEliminationNonsingularCase

In[10]:=

gaussianEliminationNonsingularCase :: Usage = "gaussianEliminationNonsingularCase[a_?sq ... s an upper triangular matrix which is row equivalent to the nonsingular square matrix a." ;

In[11]:=

gaussianEliminationNonsingularCase[a_ ? squareMatrixQ] := Module[{n = Dimensions[a][ ... ; b[[i]] = b[[i]] - c b[[j]],  {i, j + 1, n}]],  {j, 1, n}] ; b]

allZeroQ

In[12]:=

allZeroQ :: Usage = "allZeroQ[row] returns True if row is a list of zeros, and False otherwise." ;

In[13]:=

allZeroQ[row_] := If[Length[row] 1, First[row] 0, First[row] 0&&allZeroQ[Rest[row]]]

indexOfFirstNonzeroEntry

In[14]:=

indexOfFirstNonzeroEntry :: Usage = "indexOfFirstNonzeroEntry[ls] returns the index of the first nonzero entry in the list ls." ;

In[15]:=

indexOfFirstNonzeroEntry[ls_] := If[ls[[1]] ≠0, 1, 1 + indexOfFirstNonzeroEntry[Rest[ls]]] ;

Testing.

In[16]:=

a = ({{1, 2, 3}, {4, 5, 6}, {7, 8, 2}}) ; 

gaussianEliminationNonsingularCase[a] ;

%//MatrixForm

Out[18]//MatrixForm=

( {{1, 2, 3}, {0, -3, -6}, {0, 0, -7}} )

In[19]:=

a = ({{0, 2, 3}, {0, 0, 6}, {3, 8, 2}}) ; 

gaussianEliminationNonsingularCase[a] ;

%//MatrixForm

Out[21]//MatrixForm=

( {{3, 8, 2}, {0, 2, 3}, {0, 0, 6}} )

In[22]:=

a = ({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}) ; 

gaussianEliminationNonsingularCase[a]

Error: gaussianEliminationNonsingularCase: A is singular.

Out[23]=

$Aborted

Gaussian Elimination - with Partial Pivoting

The following procedure will row reduce a nonsingular square matrix to an upper triangular matrix.
The Partial Pivoting version of Gaussian Elimination uses a more sophisticated pivoting policy to reduce accummulating roundoff errors which may lead to inaccurate solutions. The idea is to always move the largest (in absolute value) available element into the pivot position. Here, the "largest available element" means the largest element in the present column which is in or below the present pivot position. A yet more wide-ranging pivoting policy could even use column operations (which effectively only permute the ordering of the variables) to bring in even larger pivots, and further reduce roundoff errors. Such a policy is called Full Pivoting, but we do not implement it here. A more sophisticated version of partial or full pivoting would use a list of pointers to rows and pointers to columns and simply interchange pointers rather than interchanging whole rows or columns in computer memory.
See Olver and Shakiban, "Applied Linear Algebra," Section 1.7.

gaussianEliminationPartialPivoting

In[24]:=

gaussianEliminationPartialPivoting :: Usage = "gaussianEliminationPartialPivoting[a_?sq ... n an upper triangular matrix which is row equivalent to the nonsingular square matrix a." ;

In[25]:=

gaussianEliminationPartialPivoting[a_ ? squareMatrixQ] := Module[{n = Dimensions[a][ ... ; b[[i]] = b[[i]] - c b[[j]],  {i, j + 1, n}]],  {j, 1, n}] ; b]

allZeroQ

In[26]:=

allZeroQ :: Usage = "allZeroQ[row] returns True if row is a list of zeros, and False otherwise." ;

In[27]:=

allZeroQ[row_] := If[Length[row] 1, First[row] 0, First[row] 0&&allZeroQ[Rest[row]]]

indexOfLargestAvailableEntry

In[28]:=

indexOfLargestAvailableEntry :: Usage = "indexOfLargestAvailableEntry[ls] returns the index of the element in the list ls having the largest absolute value." ;

In[29]:=

indexOfLargestAvailableEntry[ls_] := Module[{max = Max[Abs[ls]]}, If[Abs[ls[[1]]] max, 0, 1 + indexOfLargestAvailableEntry[Rest[ls]]]]

Testing.

In[30]:=

a = ({{1, 2, 3}, {4, 5, 6}, {7, 8, 2}}) ; 

gaussianEliminationPartialPivoting[a] ;

%//MatrixForm

Out[32]//MatrixForm=

( {{7, 8, 2}, {0, 6/7, 19/7}, {0, 0, 7/2}} )

In[33]:=

a = ({{0, 2, 3}, {0, 0, 6}, {3, 8, 2}}) ; 

gaussianEliminationPartialPivoting[a] ;

%//MatrixForm

Out[35]//MatrixForm=

( {{3, 8, 2}, {0, 2, 3}, {0, 0, 6}} )

In[36]:=

a = ({{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}) ; 

gaussianEliminationPartialPivoting[a]

Error: gaussianEliminationPartialPivoting: A is singular.

Out[37]=

$Aborted


Created by Mathematica  (April 10, 2005) Valid XHTML 1.1!