Matrix Inversion Using Gauss-Jordan Elimination

The inverse of a square matrix C is defined by the matrix equation

CC⁻¹ = C⁻¹C = I

where C⁻¹ is the inverse of C, and I is the identity matrix, for example

    |1 0 0|
I = |0 1 0|
    |0 0 1|   when C is a 3x3 matrix.

The matrix C⁻¹ exists if and only if C is nonsingular, i.e., if the determinant of C<>0.

One approach to finding C⁻¹, is to consider the matrix equation

CC⁻¹ = I (1)

and apply a series of elementary row operations to C on the left hand side of the equation and simultaneously to I on the right hand side, with the objective of transforming the C matrix into I. Ultimately, the left hand side of the matrix equation (1) then becomes IC⁻¹, or simply C⁻¹.

These same operations must transform I on the right hand side of equation (1) into C⁻¹ in order to preserve the equality. Thus equation (1) eventually becomes.

I C⁻¹ = C⁻¹ (2)

This is the Gauss Jordan elimination method.

The Program

The program given here uses the ideas above, and in addition employs several techniques for maximizing pivot elements in pivot rows and columns. In order to save storage space, the input matrix C is overwritten by the inverted matrix C⁻¹ (i.e., the original matrix C is destroyed). The program has no error exits, and if a singular matrix is entered, erroneous results are returned. Within the program however a “feasibility index” is computed, and an extremely small feasibility index is usually indicative of a singularity in which case, the matrix has no true inverse.

To run the program, the user must respond to the prompt “INPUT DIMENSION N” by entering N, the dimension of the square matrix C to be inverted, and to the prompt “INPUT C-MATRIX”, by entering the C matrix by rows, i.e., a₁₁, a₁₂, a₁₃… a₁ₙ, a₂₁, a₂₂, a₂₃… aₙₙ.

Example

Invert the matrix

    |2 1 1|
C = |1 2 1|
    |1 1 2|

When A is substituted into equation (1) we get

|2 1 1|       |1 0 0|
|1 2 1| C⁻¹ = |0 1 0|
|1 1 2|       |0 0 1|

The program will prompt: “INPUT DIMENSION N.” Since we have a 3×3 matrix, enter “3”. The next prompt is “INPUT C-MATRIX.” Now the C-matrix is entered as 2, 1, 1, 1, 2, 1, 1, 1, 2.

The final results are printed out as .75000, -.25000, -.25000, -.25000, .75000, -.25 000, -.25000, -.25000, .75000.

(Editor’s note): Ken’s future articles will feature: Numerical integration of continuous functions using the Trapezoid Rule and Simpson’s Rule.

Statistics package, including mean, mode, median and standard deviation computation from point input.

Jacobi Method for simultaneous computation of eigenvectors and eigenvalues.

100 REM MATINV
120 REM THIS PROGRAM INVERTS
130 REM INPUT MATRIX C
200 DIM C(15,15)
250 DIM J(40)
300 PRINT "INPUT DIMENSION N"
400 INPUT N
500 PRINT "INPUT C-MATRIX"
600 FOR I=1 TO N
700 FOR J=1 TO N
800 INPUT C(I,J)
900 NEXT J
1000 NEXT I
1100 LET PD=1
1200 FOR L=1 TO N
1300 LET DD=0
1400 FOR K=1 TO N
1500 LET DD=DD+C(L,K)*C(L,K)
1600 NEXT K
1700 LET DD=SQR (DD)
1800 LET PD=PD*DD
1900 NEXT L
2000 LET DETM=1
2100 FOR L=1 TO N
2200 LET J(L+20)=L
2300 NEXT L
2400 FOR L=1 TO N
2500 LET CC=0
2600 LET M=L
2700 FOR K=L TO N
2800 IF ((ABS (CC)-ABS (C(L,K))>=0)) THEN GOTO 3100
2900 LET M=K
3000 LET CC=C(L,K)
3100 NEXT K
3200 IF (L=M) THEN GOTO 4100
3300 LET K=J(M+20)
3400 LET J(M+20)=J(L+20)
3500 LET J(L+20)=K
3600 FOR K=1 TO N
3700 LET S=C(K,L)
3800 LET C(K,L)=C(K,M)
3900 LET C(K,M)=S
4000 NEXT K
4100 LET C(L,L)=1
4200 LET DETM=DETM*CC
4300 FOR M=1 TO N
4400 LET C(L,M)=C(L,M)/CC
4500 NEXT M
4600 FOR M=1 TO N
4700 IF L=M THEN GOTO 5400
4800 LET CC=C(M,L)
4900 IF (CC=0) THEN GOTO 5400
5000 LET C(M,L)=0
5100 FOR K=1 TO N
5200 LET C(M,K)=C(M,K)-CC*C(L,K)
5300 NEXT K
5400 NEXT M
5410 NEXT L
5500 FOR L=1 TO N
5600 IF (J(L+20)=L) THEN GOTO 6800
5700 LET M=L
5800 LET M=M+1
5900 IF (J(M+20)=L) THEN GOTO 6100
6000 IF (N>M) THEN GOTO 5800
6100 LET J(M+20)=J(L+20)
6200 FOR K=1 TO N
6300 LET CC=C(L,K)
6400 LET C(L,K)=C(M,K)
6500 LET C(M,K)=CC
6600 NEXT K
6700 LET J(L+20)=L
6800 NEXT L
6900 LET DETM=ABS (DETM)
7000 LET DTNRM=DETM/PD
7100 PRINT "THE FEASIBILITY INDX=",DTNRM
7200 FOR I=1 TO N
7300 FOR J=1 TO N
7400 PRINT C(I,J)
7500 NEXT J
7600 PRINT " "
7700 NEXT I
7800 STOP