36
Chapter 2. Solution of Linear Algebraic Equations
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
visit website or call 1-800-872-7423 (North America only),or send email to (outside North America).
Coleman, T.F., and Van Loan, C. 1988,
Handbook for Matrix Computations
(Philadelphia: S.I.A.M.).
Forsythe, G.E., and Moler, C.B. 1967,
Computer Solution of Linear Algebraic Systems
(Engle-
wood Cliffs, NJ: Prentice-Hall).
Wilkinson, J.H., and Reinsch, C. 1971,
Linear Algebra
,vol.IIof
Handbook for Automatic Com-
putation
(New York: Springer-Verlag).
Westlake, J.R. 1968,
A Handbook of Numerical Matrix Inversion and Solution of Linear Equations
(New York: Wiley).
Johnson, L.W., and Riess, R.D. 1982,
Numerical Analysis
, 2nd ed. (Reading, MA: Addison-
Wesley), Chapter 2.
Ralston, A., and Rabinowitz, P. 1978,
A First Course in Numerical Analysis
, 2nd ed. (New York:
McGraw-Hill), Chapter 9.

extend the equations to the case of N × N matrices, with M sets of right-hand
side vectors, in completely analogous fashion. The routine implemented below
is, of course, general.
2.1 Gauss-Jordan Elimination
37
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
visit website or call 1-800-872-7423 (North America only),or send email to (outside North America).
Elimination on Column-Augmented Matrices
Consider the linear matrix equation


a
11
a
12
a
13
a
14
a
21
a
22
a
23
a
24

41





x
12
x
22
x
32
x
42





x
13
x
23
x
33
x
43




y
43
y
44




=




b
11
b
21
b
31
b
41





b
12
b
22

(2.1.1)
Here the raised dot (·) signifies matrix multiplication, while the operator  just
signifies column augmentation, that is, removing the abutting parentheses and
making a wider matrix out of the operands of the  operator.
It should not take you long to write out equation (2.1.1) and to see that it simply
states that x
ij
is the ith component (i =1,2,3,4) of the vector solution of the jth
right-hand side (j =1,2,3), the one whose coefficients are b
ij
,i =1,2,3,4;and
that the matrix of unknown coefficients y
ij
is the inverse matrix of a
ij
.Inother
words, the matrix solution of
[A] · [x
1
 x
2
 x
3
 Y]=[b
1
b
2
b
3
1](2.1.2)

• Interchanging any two columns of A gives the same solution set only
if we simultaneously interchange corresponding rows of the x’s and of
Y. In other words, this interchange scrambles the order of the rows in
the solution. If we do this, we will need to unscramble the solution by
restoring the rows to their original order.
Gauss-Jordan elimination uses one or more of the above operations to reduce
the matrix A to the identity matrix. When this is accomplished, the right-hand side
becomes the solution set, as one sees instantly from (2.1.2).
38
Chapter 2. Solution of Linear Algebraic Equations
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
visit website or call 1-800-872-7423 (North America only),or send email to (outside North America).
Pivoting
In “Gauss-Jordan elimination with no pivoting,” only the second operation in
the above list is used. The first row is divided by the element a
11
(this being a
trivial linear combination of the first row with any other row — zero coefficient for
the other row). Then the right amount of the first row is subtracted from each other
row to make all the remaining a
i1
’s zero. The first column of A now agrees with
the identity matrix. We move to the second column and divide the second row by
a
22
, then subtract the right amount of the second row from rows 1, 3, and 4, so as to
make their entries in the second column zero. The second column is now reduced

the third linear equation in our original set and multiply it by a factor of a million, it
is almost guaranteed that it will contribute the first pivot; yet the underlying solution
of the equations is not changed by this multiplication! One therefore sometimes sees
routines which choose as pivot that element which would have been largest if the
original equations had all been scaled to have their largest coefficient normalized to
unity. This is called implicit pivoting. There is some extra bookkeeping to keep track
of the scale factors by which the rows would have been multiplied. (The routines in
§2.3 include implicit pivoting, but the routine in this section does not.)
Finally, let us consider the storage requirements of the method. With a little
reflection you will see that at every stage of the algorithm,either an element of A is
2.1 Gauss-Jordan Elimination
39
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
visit website or call 1-800-872-7423 (North America only),or send email to (outside North America).
predictably a one or zero (if it is already in a part of the matrix that has been reduced
to identity form) or else the exactly corresponding element of the matrix that started
as 1 is predictably a one or zero (if its mate in A has not been reduced to the identity
form). Therefore the matrix 1 does not have to exist as separate storage: The matrix
inverse of A is gradually built up in A as the original A is destroyed. Likewise,
the solution vectors x can gradually replace the right-hand side vectors b and share
the same storage, since after each column in A is reduced, the corresponding row
entry in the b’s is never again used.
Here is the routine for Gauss-Jordan elimination with full pivoting:
#include <math.h>
#include "nrutil.h"
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
void gaussj(float **a, int n, float **b, int m)

icol=k;
}
} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
}
++(ipiv[icol]);
We now have the pivot element, so we interchange rows, if needed, to put the pivot
element on the diagonal. The columns are not physically interchanged, only relabeled:
indxc[i], the column of the ith pivot element, is the ith column that is reduced, while
indxr[i] is the row in which that pivot element was originally located. If indxr[i] =
indxc[i] there is an implied column interchange. With this form of bookkeeping, the
solution b’s will end up in the correct order, and the inverse matrix will be scrambled
by columns.
if (irow != icol) {
for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
}
indxr[i]=irow; We are now ready to divide the pivot row by the
pivot element, located at irow and icol.indxc[i]=icol;
if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for (l=1;l<=n;l++) a[icol][l] *= pivinv;
for (l=1;l<=m;l++) b[icol][l] *= pivinv;
40
Chapter 2. Solution of Linear Algebraic Equations
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
visit website or call 1-800-872-7423 (North America only),or send email to (outside North America).



1 if i = j and i =2,4
1 if i =2,j=4
1 if i =4,j=2
0 otherwise
(2.1.5)
effects the interchange of rows 2 and 4. Gauss-Jordan elimination by row operations alone
(including the possibility of partial pivoting) consists of a series of such left-multiplications,
yielding successively
A · x = b
(···R
3
·R
2
·R
1
·A)·x =···R
3
·R
2
·R
1
·b
(1)·x =···R
3
·R
2
·R
1

· C
−1
1
· x = b
(A · C
1
· C
2
· C
3
···)···C
−1
3
·C
−1
2
·C
−1
1
·x =b
(1)···C
−1
3
·C
−1
2
·C
−1
1
·x =b