Gauss-Jordan Elimination

The first method we will look at to solve a system of equations is Gauss-Jordan elimination. As we stated previously, the general solution process is to reduce the A matrix to a form such that the system of equations can be solved directly. With Gauss-Jordan elimination, the A matrix is reduced to the identity matrix. One advantage of Gauss-Jordan is that it will also give you the inverse of the A matrix. Gauss-Jordan, when pivoted, is a very stable algorithm. One disadvantage is that it requires about three times the number of operations of Gaussian elimination or LU decomposition and thus is slower than those methods . For this reason, Gauss-Jordan elimination is less frequently used than Gaussian elimination or LU decomposition.

In addition to solving the system of equations from Eq. (19.2), Gauss-Jordan elimination can be used to simultaneously solve for the inverse of the A matrix starting from Eq. (19.6).

Equation 19.6

In Eq. (19.6), A ^“ ¹ is the inverse of the A matrix and I is the identity matrix.

Equation 19.7

graphics/19equ07.gif

What Gauss-Jordan elimination does is to convert the original A matrix in both Eq. (19.2) and (19.6) to the identity matrix by adding multiples of one row to another. When this is done, we are left with the equations in Eq. (19.8).

Equation 19.8

The original b vector of Eq. (19.2) is overwritten to contain the vector of unknowns. The identity matrix on the right-hand side of Eq. (19.6) becomes the inverse A matrix. You don't have to allocate memory to store the identity (and later inverse A ) matrix on the right-hand side of Eq. (19.6). The A matrix is reduced column by column, and while this is going on the inverse A matrix is being formed column by column. The Gauss-Jordan algorithm is written such that the columns of the A ^“ ¹ matrix replace the columns of the A matrix one by one. When the solution process is complete, the original A matrix now contains its inverse.

The Gauss-Jordan method is unstable unless pivoting is used. When you pivot the original A matrix (and perform the same operations on the b vector), the resulting A ^“ ¹ matrix will be scrambled. If the A matrix were partially pivoted, the inverse matrix can be unscrambled by swapping its columns in the opposite order that you swapped the rows of the original A matrix.

We will implement the Gauss-Jordan elimination solution process in a method named gaussJordan() . There are three parts to the method. First, the partialPivot() method is called to perform partial pivoting on the A matrix. Next, the pivoted A matrix is reduced to the identity matrix. The same operations are performed on the b vector. Finally, the A matrix, which now contains the inverse A matrix, is unscrambled. The code listing for the gaussJordan() method is shown here.

 public static void gaussJordan(double a[][], double b[]) {   int i,j,k,m;   double temp;   int numRows = a.length;   int numCols = a[0].length;   int index[][] = new int[numRows][2];   //  Perform an implicit partial pivoting of the   //  a[][] array and b[]  array.   partialPivot(a, b, index);   //  Perform the elimination row by row. First dividing   //  the current row and b element by a[i][i]   for(i=0; i<numRows; ++i) {     temp = a[i][i];     for(j=0; j<numCols; ++j) {       a[i][j] /= temp;     }     b[i] /= temp;     a[i][i] = 1.0/temp;     //  Reduce the other rows by subtracting a multiple     //  of the current row from them. Don't reduce the     //  current row. As each column of the a[][] matrix     //  is reduced its elements are replaced with the     //  inverse a[][] matrix.     for(k=0; k<numRows; ++k) {       if (k != i) {         temp = a[k][i];         for(j=0; j<numCols; ++j) {           a[k][j] -= temp*a[i][j];         }         b[k] -= temp*b[i];         a[k][i] = -temp*a[i][i];       }     }   }   //  Unscramble the inverse a[][] matrix.   //  The columns are swapped in the opposite order   //  that the rows were during the pivoting.   for(j=numCols-1; j>=0; j) {     k = index[j][0];     m = index[j][1];     if (k != m) {       for(i=0; i<numRows; ++i) {         temp = a[i][m];         a[i][m] = a[i][k];         a[i][k] = temp;       }     }   }   return; }

The gaussJordan() method is tested in the "Testing the EqnSolver Class Methods" section later in this chapter.