44 Chapter 2. Solution of Linear Algebraic Equations
Equations (2.3.6) and (2.3.7) total (for each right-hand side b) N 2 executions
of an inner loop containing one multiply and one add. If we have N right-hand
sides which are the unit column vectors (which is the case when we are inverting a
matrix), then taking into account the leading zeros reduces the total execution count
of (2.3.6) from 1 N 3 to 1 N 3 , while (2.3.7) is unchanged at 1 N 3 .
2 6 2
Notice that, once we have the LU decomposition of A, we can solve with as
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
many right-hand sides as we then care to, one at a time. This is a distinct advantage
over the methods of §2.1 and §2.2.
Performing the LU Decomposition
How then can we solve for L and U, given A? First, we write out the
i, jth component of equation (2.3.1) or (2.3.2). That component always is a sum
beginning with
αi1 β1j + · · · = aij
The number of terms in the sum depends, however, on whether i or j is the smaller
number. We have, in fact, the three cases,
ij: αi1 β1j + αi2 β2j + · · · + αij βjj = aij (2.3.10)
Equations (2.3.8)–(2.3.10) total N 2 equations for the N 2 + N unknown α’s and
β’s (the diagonal being represented twice). Since the number of unknowns is greater
than the number of equations, we are invited to specify N of the unknowns arbitrarily
and then try to solve for the others. In fact, as we shall see, it is always possible to take
αii ≡ 1 i = 1, . . . , N (2.3.11)
A surprising procedure, now, is Crout’s algorithm, which quite trivially solves
the set of N 2 + N equations (2.3.8)–(2.3.11) for all the α’s and β’s by just arranging
the equations in a certain order! That order is as follows:
• Set αii = 1, i = 1, . . . , N (equation 2.3.11).
• For each j = 1, 2, 3, . . . , N do these two procedures: First, for i =
1, 2, . . ., j, use (2.3.8), (2.3.9), and (2.3.11) to solve for βij , namely
i−1
βij = aij − αik βkj . (2.3.12)
k=1
(When i = 1 in 2.3.12 the summation term is taken to mean zero.) Second,
for i = j + 1, j + 2, . . . , N use (2.3.10) to solve for αij , namely
j−1
1
αij = aij − αik βkj . (2.3.13)
βjj
k=1
Be sure to do both procedures before going on to the next j.

2.3 LU Decomposition and Its Applications 45
a
c
e
g
i
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
etc.
x
b
d
f di
ag
h on
j al
su el
bd em
ia en
go ts
x na
le
le
m
etc. en
ts
Figure 2.3.1. Crout’s algorithm for LU decomposition of a matrix. Elements of the original matrix are
modiﬁed in the order indicated by lower case letters: a, b, c, etc. Shaded boxes show the previously
modiﬁed elements that are used in modifying two typical elements, each indicated by an “x”.
If you work through a few iterations of the above procedure, you will see that
the α’s and β’s that occur on the right-hand side of equations (2.3.12) and (2.3.13)
are already determined by the time they are needed. You will also see that every aij
is used only once and never again. This means that the corresponding αij or βij can
be stored in the location that the a used to occupy: the decomposition is “in place.”
[The diagonal unity elements αii (equation 2.3.11) are not stored at all.] In brief,
Crout’s method ﬁlls in the combined matrix of α’s and β’s,
 
β11 β12 β13 β14
 α21 β22 β23 β24 
  (2.3.14)
α31 α32 β33 β34
α41 α42 α43 β44
by columns from left to right, and within each column from top to bottom (see
Figure 2.3.1).
What about pivoting? Pivoting (i.e., selection of a salubrious pivot element
for the division in equation 2.3.13) is absolutely essential for the stability of Crout’s

46 Chapter 2. Solution of Linear Algebraic Equations
method. Only partial pivoting (interchange of rows) can be implemented efﬁciently.
However this is enough to make the method stable. This means, incidentally, that
we don’t actually decompose the matrix A into LU form, but rather we decompose
a rowwise permutation of A. (If we keep track of what that permutation is, this
decomposition is just as useful as the original one would have been.)
Pivoting is slightly subtle in Crout’s algorithm. The key point to notice is that
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
equation (2.3.12) in the case of i = j (its ﬁnal application) is exactly the same as
equation (2.3.13) except for the division in the latter equation; in both cases the
upper limit of the sum is k = j − 1 (= i − 1). This means that we don’t have to
commit ourselves as to whether the diagonal element βjj is the one that happens
to fall on the diagonal in the ﬁrst instance, or whether one of the (undivided) αij ’s
below it in the column, i = j + 1, . . . , N , is to be “promoted” to become the diagonal
β. This can be decided after all the candidates in the column are in hand. As you
should be able to guess by now, we will choose the largest one as the diagonal β
(pivot element), then do all the divisions by that element en masse. This is Crout’s
method with partial pivoting. Our implementation has one additional wrinkle: It
initially ﬁnds the largest element in each row, and subsequently (when it is looking
for the maximal pivot element) scales the comparison as if we had initially scaled all
the equations to make their maximum coefﬁcient equal to unity; this is the implicit
pivoting mentioned in §2.1.
#include
#include "nrutil.h"
#define TINY 1.0e-20; A small number.
void ludcmp(float **a, int n, int *indx, float *d)
Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise
permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above;
indx[1..n] is an output vector that records the row permutation eﬀected by the partial
pivoting; d is output as ±1 depending on whether the number of row interchanges was even
or odd, respectively. This routine is used in combination with lubksb to solve linear equations
or invert a matrix.
{
int i,imax,j,k;
float big,dum,sum,temp;
float *vv; vv stores the implicit scaling of each row.
vv=vector(1,n);
*d=1.0; No row interchanges yet.
for (i=1;i

2.3 LU Decomposition and Its Applications 47
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big) {
Is the ﬁgure of merit for the pivot better than the best so far?
big=dum;
imax=i;
}
}
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
if (j != imax) { Do we need to interchange rows?
for (k=1;k

48 Chapter 2. Solution of Linear Algebraic Equations
The LU decomposition in ludcmp requires about 1 N 3 executions of the inner
3
loops (each with one multiply and one add). This is thus the operation count
for solving one (or a few) right-hand sides, and is a factor of 3 better than the
Gauss-Jordan routine gaussj which was given in §2.1, and a factor of 1.5 better
than a Gauss-Jordan routine (not given) that does not compute the inverse matrix.
For inverting a matrix, the total count (including the forward and backsubstitution
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
as discussed following equation 2.3.7 above) is ( 1 + 1 + 1 )N 3 = N 3 , the same
3 6 2
as gaussj.
To summarize, this is the preferred way to solve the linear set of equations
A · x = b:
float **a,*b,d;
int n,*indx;
...
ludcmp(a,n,indx,&d);
lubksb(a,n,indx,b);
The answer x will be given back in b. Your original matrix A will have
been destroyed.
If you subsequently want to solve a set of equations with the same A but a
different right-hand side b, you repeat only
lubksb(a,n,indx,b);
not, of course, with the original matrix A, but with a and indx as were already
set by ludcmp.
Inverse of a Matrix
Using the above LU decomposition and backsubstitution routines, it is com-
pletely straightforward to ﬁnd the inverse of a matrix column by column.
#define N ...
float **a,**y,d,*col;
int i,j,*indx;
...
ludcmp(a,N,indx,&d); Decompose the matrix just once.
for(j=1;j

2.3 LU Decomposition and Its Applications 49
Incidentally, if you ever have the need to compute A−1 · B from matrices A
and B, you should LU decompose A and then backsubstitute with the columns of
B instead of with the unit vectors that would give A’s inverse. This saves a whole
matrix multiplication, and is also more accurate.
Determinant of a Matrix
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
The determinant of an LU decomposed matrix is just the product of the
diagonal elements,
N
det = βjj (2.3.15)
j=1
We don’t, recall, compute the decomposition of the original matrix, but rather a
decomposition of a rowwise permutation of it. Luckily, we have kept track of
whether the number of row interchanges was even or odd, so we just preface the
product by the corresponding sign. (You now ﬁnally know the purpose of setting
d in the routine ludcmp.)
Calculation of a determinant thus requires one call to ludcmp, with no subse-
quent backsubstitutions by lubksb.
#define N ...
float **a,d;
int j,*indx;
...
ludcmp(a,N,indx,&d); This returns d as ±1.
for(j=1;j

50 Chapter 2. Solution of Linear Algebraic Equations
A quick-and-dirty way to solve complex systems is to take the real and imaginary
parts of (2.3.16), giving
A·x−C·y=b
(2.3.17)
C·x+A·y=d
which can be written as a 2N × 2N set of real equations,
visit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).
readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMs
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machine-
Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.
Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)
A −C x b
· = (2.3.18)
C A y d
and then solved with ludcmp and lubksb in their present forms. This scheme is a factor of
2 inefﬁcient in storage, since A and C are stored twice. It is also a factor of 2 inefﬁcient in
time, since the complex multiplies in a complexiﬁed version of the routines would each use
4 real multiplies, while the solution of a 2N × 2N problem involves 8 times the work of
an N × N one. If you can tolerate these factor-of-two inefﬁciencies, then equation (2.3.18)
is an easy way to proceed.
CITED REFERENCES AND FURTHER READING:
Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins
University Press), Chapter 4.
Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.).
Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical
Computations (Englewood Cliffs, NJ: Prentice-Hall), §3.3, and p. 50.
Forsythe, G.E., and Moler, C.B. 1967, Computer Solution of Linear Algebraic Systems (Engle-
wood Cliffs, NJ: Prentice-Hall), Chapters 9, 16, and 18.
Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations
(New York: Wiley).
Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),
§4.2.
Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:
McGraw-Hill), §9.11.
Horn, R.A., and Johnson, C.R. 1985, Matrix Analysis (Cambridge: Cambridge University Press).
2.4 Tridiagonal and Band Diagonal Systems
of Equations
The special case of a system of linear equations that is tridiagonal, that is, has
nonzero elements only on the diagonal plus or minus one column, is one that occurs
frequently. Also common are systems that are band diagonal, with nonzero elements
only along a few diagonal lines adjacent to the main diagonal (above and below).
For tridiagonal sets, the procedures of LU decomposition, forward- and back-
substitution each take only O(N ) operations, and the whole solution can be encoded
very concisely. The resulting routine tridag is one that we will use in later chapters.
Naturally, one does not reserve storage for the full N × N matrix, but only for
the nonzero components, stored as three vectors. The set of equations to be solved is
     
b 1 c1 0 · · · u1 r1
 a2 b 2 c2 · · ·   u2   r2 
     
 ···  ·  · · ·  =  · · ·  (2.4.1)
     
· · · aN−1 bN−1 cN−1 uN−1 rN−1
··· 0 aN bN uN rN