2.3 LU Decomposition and Its Applications 43 Isaacson,E.,and Keller,H.B.1966,Analysis of Numerical Methods (New York:Wiley),82.1. Johnson,L.W.,and Riess,R.D.1982,Numerica/Analysis,2nd ed.(Reading,MA:Addison- Wesley),$2.2.1. Westlake,J.R.1968,A Handbook of Numerical Matrix Inversion and Solution of Linear Equations (New York:Wiley). 2.3 LU Decomposition and Its Applications Suppose we are able to write the matrix A as a product of two matrices, 8= L.U=A (2.3.1) where L is lower triangular (has elements only on the diagonal and below)and U is upper triangular (has elements only on the diagonal and above).For the case of a 4 x 4 matrix A,for example,equation(2.3.1)would look like this: 011 0 0 0 311 12 613 B14 [a11 a12 a13 a14 server 2 a21 022 0 0 022 23 24 a21 a22 a23 a24 (Nort a31 a32 a33 0 0 0 334 a31 a32 a33 a34 041 a42 a43 044 0 0 344 a41 a42 a43 a44 America University Press. THE (2.3.2) ART We can use a decomposition such as(2.3.1)to solve the linear set Progra A·x=(LU)x=L·(Ux)=b (2.3.3) OF SCIENTIFIC( by first solving for the vector y such that 61 L·y=b (2.3.4) and then solving U·x=y (2.3.5) COMPUTING (ISBN 188810920 What is the advantage of breaking up one linear set into two successive ones? The advantage is that the solution of a triangular set of equations is quite trivial,as we have already seen in 82.2 (equation 2.2.4).Thus,equation(2.3.4)can be solved by forward substitution as follows, uurggoglrion Numerical Recipes 10-621 43106 b1 1= 011 (outside i- (2.3.6) North Software. = bi i=2,3,,N j=1 while(2.3.5)can then be solved by backsubstitution exactly as in equations(2.2.2)- (2.2.4), yN TN= BNN (2.3.7) i= i=N-1,N-2,,1 Bu
2.3 LU Decomposition and Its Applications 43 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §2.1. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §2.2.1. Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations (New York: Wiley). 2.3 LU Decomposition and Its Applications Suppose we are able to write the matrix A as a product of two matrices, L · U = A (2.3.1) where L is lower triangular (has elements only on the diagonal and below) and U is upper triangular (has elements only on the diagonal and above). For the case of a 4 × 4 matrix A, for example, equation (2.3.1) would look like this: α11 000 α21 α22 0 0 α31 α32 α33 0 α41 α42 α43 α44 · β11 β12 β13 β14 0 β22 β23 β24 0 0 β33 β34 000 β44 = a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44 (2.3.2) We can use a decomposition such as (2.3.1) to solve the linear set A · x = (L · U) · x = L · (U · x) = b (2.3.3) by first solving for the vector y such that L · y = b (2.3.4) and then solving U · x = y (2.3.5) What is the advantage of breaking up one linear set into two successive ones? The advantage is that the solution of a triangular set of equations is quite trivial, as we have already seen in §2.2 (equation 2.2.4). Thus, equation (2.3.4) can be solved by forward substitution as follows, y1 = b1 α11 yi = 1 αii bi − i−1 j=1 αijyj i = 2, 3,...,N (2.3.6) while (2.3.5) can then be solved by backsubstitution exactly as in equations (2.2.2)– (2.2.4), xN = yN βNN xi = 1 βii yi − N j=i+1 βijxj i = N − 1, N − 2,..., 1 (2.3.7)
44 Chapter 2.Solution of Linear Algebraic Equations Equations(2.3.6)and(2.3.7)total (for each right-hand side b)N2 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 N3 to N3,while (2.3.7)is unchanged at N3 Notice that,once we have the LU decomposition of A,we can solve with as 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,ith component of equation(2.3.1)or(2.3.2).That component always is a sum beginning with ICAL 01B1j十···=ai对 RECIPES 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. 9 ij: a1B1j+a2P2j+…+afi=a (2.3.10) 色 9 Equations(2.3.8)-(2.3.10)total N2 equations for the N2+N unknown a's and B's(the diagonal being represented twice).Since the number of unknowns is greater than the number ofequations,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 ai≡1i=1,.,N (2.3.11) A surprising procedure,now,is Crout's algorithm,which quite trivially solves the set of N2+N equations(2.3.8)-(2.3.11)for all the a's and B's by just arranging 兰丝母 10621 the equations in a certain order!That order is as follows: Numerica .Set ai =1,i =1,...,N (equation 2.3.11). 43126 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 Bij,namely i1 North =a1 (2.3.12) k= (Wheni=1 in 2.3.12 the summation term is taken to mean zero.)Second. fori=j+1,j+2,...,N use (2.3.10)to solve for aii,namely 2.3.13) Be sure to do both procedures before going on to the next j
44 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 2N3 to 1 6N3, while (2.3.7) is unchanged at 1 2N3. Notice that, once we have the LU decomposition of A, we can solve with as 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 N2 + 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 βij = aij − i−1 k=1 αikβkj . (2.3.12) (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 αij = 1 βjj aij − j−1 k=1 αikβkj . (2.3.13) Be sure to do both procedures before going on to the next j
2.3 LU Decomposition and Its Applications 45 a etc. http://www.nr.com or call 1-800-872- (including this one) internet ( ① diagonal elements subdiagonal elements -7423(North America to any server computer,is users to make one paper etc. only),or send to dir Copyright (C) from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING(ISBN 1988-1992 by Cambridge University Press.Programs Figure 2.3.1.Crout's algorithm for LU decomposition of a matrix.Elements of the original matrix are modified in the order indicated by lower case letters:a,b,c.etc.Shaded boxes show the previously modified elements that are used in modifying two typical elements,each indicated by an"x". 1788-1982 If you work through a few iterations of the above procedure,you will see that the a's and B's that occur on the right-hand side of equations(2.3.12)and(2.3.13) .Further reproduction, Numerical Recipes 10-621 are already determined by the time they are needed.You will also see that every a is used only once and never again.This means that the corresponding a ij or Bij can be stored in the location that the a used to occupy:the decomposition is"in place." 43195 [The diagonal unity elements a(equation 2.3.11)are not stored at all.]In brief, (outside Crout's method fills in the combined matrix of a's and B's. North Software. 11 012613 3141 ing of 021 322323 524 (2.3.14) Q31 032 333 34 visit website machine L041 Q42 043 344」 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
2.3 LU Decomposition and Its Applications 45 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). c g i b d f h j diagonal elements subdiagonal elements etc. etc. x x a e Figure 2.3.1. Crout’s algorithm for LU decomposition of a matrix. Elements of the original matrix are modified in the order indicated by lower case letters: a, b, c, etc. Shaded boxes show the previously modified 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 a ij 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 fills in the combined matrix of α’s and β’s, β11 β12 β13 β14 α21 β22 β23 β24 α31 α32 β33 β34 α41 α42 α43 β44 (2.3.14) 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 efficiently. 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 equation (2.3.12)in the case of i=j (its final 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 isk =j-1(=i-1).This means that we don't have to commit ourselves as to whether the diagonal element Bi is the one that happens 8 to fall on the diagonal in the first instance,or whether one of the (undivided)ai's below it in the column,i=j+1,...,N,is to be"promoted"to become the diagonal B.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 B (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 MERICAI initially finds 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 coefficient equal to unity;this is the implicit pivoting mentioned in $2.1. Press. THE ART #include #include "nrutil.h" 9 #define TINY 1.0e-20 A small number. Program 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 effected by the partial to dir 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. OF SCIENTIFIC COMPUTING(ISBN 1988-19920 int i,imax,,k; float big,dum,sum,tempi float *vv; vv stores the implicit scaling of each row vv=vector(1,n); *d=1.0; No row interchanges yet. Numerical Recipes 10-621 43106 for(i=1;1big)big-temp; oftware. if (big ==0.0)nrerror("Singular matrix in routine ludcmp"); No nonzero largest element. v[1]=1.0/b1g; Save the scaling. for(j=1;j<=n;j++){ This is the loop over columns of Crout's method. for (i=1;i<j;i++) This is equation (2.3.12)except for i=j. sum=a[i][j]; for (k=1;k<i;k++)sum -a[i][k]*a[k][j]; a[i][j]=sum; b1g=0.0 Initialize for the search for largest pivot element. for (i=j;i<-n;i++){ This is i=j of equation(2.3.12)andi=j+1...N sum=a[i][j]; of equation (2.3.13). for (k=1;k<j;k++)
46 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). method. Only partial pivoting (interchange of rows) can be implemented efficiently. 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 equation (2.3.12) in the case of i = j (its final 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 first 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 finds 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 coefficient 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 effected 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 big) big=temp; if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); No nonzero largest element. vv[i]=1.0/big; Save the scaling. } for (j=1;j<=n;j++) { This is the loop over columns of Crout’s method. for (i=1;i<j;i++) { This is equation (2.3.12) except for i = j. sum=a[i][j]; for (k=1;k<i;k++) sum -= a[i][k]*a[k][j]; a[i][j]=sum; } big=0.0; Initialize for the search for largest pivot element. for (i=j;i<=n;i++) { This is i = j of equation (2.3.12) and i = j + 1 ...N sum=a[i][j]; of equation (2.3.13). for (k=1;k<j;k++)
2.3 LU Decomposition and Its Applications 47 sum -a[i][k]*a[k][j]; a[i][i]=sum; if (dum=vv[i]*fabs(sum))>=big){ Is the figure of merit for the pivot better than the best so far? big=dum; imax=i; if (i !imax){ Do we need to interchange rows? for(k=1;k=1;1-)[ Now we do the backsubstitution,equation (2.3.7). sum=b[i]; for (j=i+1;j<=n;j++)sum -a[i][j]*b[j]; b[i]=sum/a[i][i]; Store a component of the solution vector X. All done!
2.3 LU Decomposition and Its Applications 47 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). sum -= a[i][k]*a[k][j]; a[i][j]=sum; if ( (dum=vv[i]*fabs(sum)) >= big) { Is the figure of merit for the pivot better than the best so far? big=dum; imax=i; } } if (j != imax) { Do we need to interchange rows? for (k=1;k=1;i--) { Now we do the backsubstitution, equation (2.3.7). sum=b[i]; for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; b[i]=sum/a[i][i]; Store a component of the solution vector X. } All done! }
48 Chapter 2.Solution of Linear Algebraic Equations The LU decomposition in ludcmp requires about N3 executions of the inner 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 82.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 as discussed following equation 2.3.7 above)is (++)N3N3,the same as gaussj To summarize,this is the preferred way to solve the linear set of equations A·x=b nted for 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 R 令 Press. different right-hand side b,you repeat only 9 lubksb(a,n,indx,b); SCIENTIFIC not,of course,with the original matrix A.but with a and indx as were already 6 set by ludcmp. Inverse of a Matrix 1920 COMPUTING (ISBN Using the above L0 decomposition and backsubstitution routines,it is com- Numerical Recipes 021 pletely straightforward to find the inverse of a matrix column by column. 43108 #define N… float **a,**y,d,*col; int i,j,*indx; (outside Software. ludcmp(a,N,indx,&d); Decompose the matrix just once. for(j=1;j<=N;j++){ Find inverse by columns. for(1=1;i<=N;i++) co1[1]=0.0: Amer ying of co1[j]=1.0: lubksb(a,N,indx,col); for(i=1;i<=N;i+)y[i][j]=co1[1]; The matrix y will now contain the inverse of the original matrix a,which will have been destroyed.Alternatively,there is nothing wrong with using a Gauss-Jordan routine like gaussj (82.1)to invert a matrix in place,again destroying the original Both methods have practically the same operations count
48 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The LU decomposition in ludcmp requires about 1 3N3 executions of the inner 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 as discussed following equation 2.3.7 above) is ( 1 3 + 1 6 + 1 2 )N3 = N3, the same 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 completely straightforward to find 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<=N;j++) { Find inverse by columns. for(i=1;i<=N;i++) col[i]=0.0; col[j]=1.0; lubksb(a,N,indx,col); for(i=1;i<=N;i++) y[i][j]=col[i]; } The matrix y will now contain the inverse of the original matrix a, which will have been destroyed. Alternatively, there is nothing wrong with using a Gauss-Jordan routine like gaussj (§2.1) to invert a matrix in place, again destroying the original. Both methods have practically the same operations count
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 The determinant of an Lo decomposed matrix is just the product of the diagonal elements, det (2.3.15) 1=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 RECIPES I 2 product by the corresponding sign.(You now finally know the purpose of setting d in the routine ludcmp.) Calculation of a determinant thus requires one call to ludcmp,with no subse- UnN电.U quent backsubstitutions by lubksb. #define N... Programs float米*a,d; int j,*indx; SCIENTIFIC ludcmp(a,N,indx,&d); This returns d as士l. for(j=1;j<=;j+)d*a[j][j]; The variable d now contains the determinant of the original matrix a,which will have been destroyed. COMPUTING (ISBN 19891892 For a matrix of any substantial size,it is quite likely that the determinant will overflow or underflow your computer's floating-point dynamic range.In this case you can modify the loop of the above fragment and (e.g.divide by powers of ten, to keep track of the scale separately,or(e.g.)accumulate the sum of logarithms of Fuurggoglrion Numerical Recipes 10-621 43106 the absolute values of the factors and the sign separately. (outside Complex Systems of Equations North If your matrix A is real,but the right-hand side vector is complex,say b+id,then (i) LU decompose A in the usual way,(i)backsubstitute b to get the real part of the solution vector,and(iii)backsubstitute d to get the imaginary part of the solution vector. If the matrix itself is complex,so that you want to solve the system (A+C)·(x+y)=(b+id) (2.3.16) then there are two possible ways to proceed.The best way is to rewrite ludcmp and lubksb as complex routines.Complex modulus substitutes for absolute value in the construction of the scaling vector vv and in the search for the largest pivot elements.Everything else goes through in the obvious way,with complex arithmetic used as needed.(See $81.2 and 5.4 for discussion of complex arithmetic in C.)
2.3 LU Decomposition and Its Applications 49 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 The determinant of an LU decomposed matrix is just the product of the diagonal elements, det = N j=1 βjj (2.3.15) 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 finally know the purpose of setting d in the routine ludcmp.) Calculation of a determinant thus requires one call to ludcmp, with no subsequent 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<=N;j++) d *= a[j][j]; The variable d now contains the determinant of the original matrix a, which will have been destroyed. For a matrix of any substantial size, it is quite likely that the determinant will overflow or underflow your computer’s floating-point dynamic range. In this case you can modify the loop of the above fragment and (e.g.) divide by powers of ten, to keep track of the scale separately, or (e.g.) accumulate the sum of logarithms of the absolute values of the factors and the sign separately. Complex Systems of Equations If your matrix A is real, but the right-hand side vector is complex, say b + id, then (i) LU decompose A in the usual way, (ii) backsubstitute b to get the real part of the solution vector, and (iii) backsubstitute d to get the imaginary part of the solution vector. If the matrix itself is complex, so that you want to solve the system (A + iC) · (x + iy)=(b + id) (2.3.16) then there are two possible ways to proceed. The best way is to rewrite ludcmp and lubksb as complex routines. Complex modulus substitutes for absolute value in the construction of the scaling vector vv and in the search for the largest pivot elements. Everything else goes through in the obvious way, with complex arithmetic used as needed. (See §§1.2 and 5.4 for discussion of complex arithmetic in C.)
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 x 2N set of real equations, (è)(（)=() (2.3.18) and then solved with ludcmp and lubksb in their present forms.This scheme is a factor of 2 inefficient in storage,since A and C are stored twice.It is also a factor of 2 inefficient in time,since the complex multiplies in a complexified version of the routines would each 81 use 4 real multiplies,while the solution of a 2N x 2N problem involves 8 times the work of an N x N one.If you can tolerate these factor-of-two inefficiencies,then equation (2.3.18) is an easy way to proceed. ICAL 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. RECIPES 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),83.3,and p.50. Forsythe.G.E.,and Moler,C.B.1967.Computer Solution of Linear Algebraic Systems(Engle- 代》 Press. 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). 84.2. SCIENTIFIC Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),89.11. 6 Horn,R.A.,and Johnson,C.R.1985,Matrix Analysis(Cambridge:Cambridge University Press). 2.4 Tridiagonal and Band Diagonal Systems of Equations Numerica 10621 uctio 43106 The special case of a system of linear equations that is tridiagonal,that is,has (outside Recipes 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). North Software. 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 x N matrix,but only for the nonzero components,stored as three vectors.The set ofequations to be solved is [b1c10- u1 a2b2c2··· u2 T2 (2.4.1) ·aN-1bN-1 CN-1 UN-1 TN-1 . uN TN
50 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin 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) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 C · x + A · y = d (2.3.17) which can be written as a 2N × 2N set of real equations, A −C C A · x y = b d (2.3.18) and then solved with ludcmp and lubksb in their present forms. This scheme is a factor of 2 inefficient in storage, since A and C are stored twice. It is also a factor of 2 inefficient in time, since the complex multiplies in a complexified 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 inefficiencies, 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 (Englewood 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 backsubstitution 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 b1 c1 0 ··· a2 b2 c2 ··· ··· ··· aN−1 bN−1 cN−1 ··· 0 aN bN · u1 u2 ··· uN−1 uN = r1 r2 ··· rN−1 rN (2.4.1)
