f2-1

2.1Gauss-JordanElimination 27Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press). Gill,P.E., Murray, W.,andWright, M.H. 1991, NumericalLinearAlgebraandOptimization , vol. 1 (Redwood City, CA: Addison-Wesley). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 4. Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.). Coleman,T.F.,andVanLoan,C.1988, HandbookforMatrixComputations (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. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (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. 2.1 Gauss-Jordan Elimination For inverting a matrix, Gauss-Jordan elimination is about as efficient as any other method. For solving sets of linear equations, Gauss-Jordan elimination produces boththe solution of the equations for one or more right-handside vectors b, and also the matrixinverse A−1. However,its principalweaknesses are (i) that it requires all the right-hand sides to be stored and manipulated at the same time, and (ii) that when the inverse matrix is notdesired, Gauss-Jordan is three times slower thanthebestalternativetechniqueforsolvingasinglelinearset( §2.3). Themethod’s principal strength is that it is as stable as any other direct method, perhaps even a bit more stable when full pivoting is used (see below). If you come along later with an additional right-hand side vector, you can multiplyitbytheinversematrix,ofcourse. Thisdoesgiveananswer,butonethatis quite susceptible to roundofferror,not nearlyas goodas if the new vector had beenincluded with the set of right-hand side vectors in the first instance. Forthesereasons,Gauss-Jordaneliminationshouldusuallynotbeyourmethod of first choice, either for solving linear equations or for matrix inversion. The decompositionmethodsin §2.3arebetter. WhydowegiveyouGauss-Jordanatall? Because it is straightforward, understandable, solid as a rock, and an exceptionallygood“psychological”backupforthosetimesthatsomethingisgoingwrongandyou think itmightbe your linear-equation solver. Some people believe that the backup is more than psychological, that Gauss- Jordan elimination is an “independent” numerical method. This turns out to be mostly myth. Except for the relatively minor differences in pivoting, described below, the actual sequence of operations performed in Gauss-Jordan elimination is very closely related to that performedby the routines in the next two sections. 28 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Forclarity,andtoavoidwritingendlessellipses( ···)wewillwriteoutequations onlyforthecaseoffourequationsandfourunknowns,andwiththreedifferentright-hand side vectors that are known in advance. You can write bigger matrices and extend the equations to the case of N×Nmatrices, with Msets of right-hand side vectors, in completely analogous fashion. The routine implemented below is,of course, general. EliminationonColumn-AugmentedMatrices Consider the linear matrix equation  a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 a41 a42 a43 a44·  x11 x21 x31 x41/unionsq x12 x22 x32 x42/unionsq x13 x23 x33 x43/unionsq y11 y12 y13 y14 y21 y22 y23 y24 y31 y32 y33 y34 y41 y42 y43 y44  =  b11 b21 b31 b41/unionsq b12 b22 b32 b42/unionsq b13 b23 b33 b43/unionsq 1000 0100 0010 0001   (2.1.1 ) Here the raised dot ( ·) signifies matrix multiplication, while the operator /unionsqjust signifies column augmentation, that is, removing the abutting parentheses andmaking a wider matrix out of the operands of the /unionsqoperator. Itshouldnottakeyoulongtowriteoutequation(2.1.1)andtoseethatitsimply states that x ijis 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 bij,i=1 ,2,3,4; and that the matrix of unknown coefficients yijis the inverse matrix of aij. In other words, the matrix solution of [A]·[x1/unionsqx2/unionsqx3/unionsqY]=[b1/unionsqb2/unionsqb3/unionsq1]( 2.1.2 ) whereAandYare square matrices, the bi’s andxi’s are column vectors, and 1is the identity matrix, simultaneously solves the linear sets A·x1=b1A·x2=b2A·x3=b3 (2.1.3 ) and A·Y=1 (2.1.4 ) Now it is also elementary to verify the following facts about (2.1.1): •Interchanging any two rowsofAand the corresponding rowsof theb’s and of1, does not change (or scramble in any way) the solution x’s and Y. Rather, it just corresponds to writing the same set of linear equations in a different order. •Likewise, the solution set is unchanged and in no way scrambled if we replace any row in Aby a linear combination of itself and any other row, as longas we do the same linear combinationofthe rows of the b’s and1 (which then is no longer the identity matrix, of course). 2.1Gauss-JordanElimination 29Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).•Interchanging any two columnsofAgives the same solution set only if we simultaneously interchange corresponding rowsof thex’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 Ato the identity matrix. When this is accomplished, the right-hand side becomes the solution set, as one sees instantly from (2.1.2). 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 a11(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 ai1’s zero. The first column of Anow agrees with the identity matrix. We move to the second column and divide the second row bya 22, thensubtracttherightamountofthe secondrowfromrows1,3, and4,so as to make their entries in the second column zero. The second column is now reduced to the identity form. And so on for the third and fourth columns. As we do theseoperationsto A, we ofcourse also do the correspondingoperationsto the b’s and to 1(which by now no longer resembles the identity matrix in any way!). Obviously we will run into trouble if we ever encounter a zero element on the (then current) diagonal when we are going to divide by the diagonal element. (The element that we divide by, incidentally,is called the pivot element orpivot.) Not so obvious,buttrue,isthefactthatGauss-Jordaneliminationwithnopivoting(nouseof thefirst orthirdproceduresin theabovelist) is numericallyunstablein thepresence ofanyroundofferror,evenwhenazeropivotisnotencountered. Youmust neverdo Gauss-Jordanelimination(or Gaussian elimination,see below)without pivoting! Sowhatisthismagicpivoting? Nothingmorethaninterchangingrows( partial pivoting) or rows and columns ( full pivoting ), so as to put a particularly desirable element in the diagonalpositionfrom whichthe pivot is aboutto be selected. Since wedon’twanttomessupthepartoftheidentitymatrixthatwehavealreadybuiltup,we can choose among elements that are both (i) on rows below (or on) the one that is about to be normalized, and also (ii) on columns to the right (or on) the column we are about to eliminate. Partial pivoting is easier than full pivoting, because wedon’t have to keep track of the permutation of the solution vector. Partial pivoting makes available as pivots only the elements already in the correct column. It turns out that partial pivoting is “almost” as good as full pivoting, in a sense that can be made mathematically precise, but which need not concern us here (for discussion andreferences,see [1]). Toshowyoubothvariants,wedofullpivotingintheroutine in this section, partial pivoting in §2.3. We have to state how to recognize a particularly desirable pivot when we see one. The answer to this is not completely known theoretically. It is known, boththeoreticallyandinpractice,thatsimplypickingthelargest(inmagnitude)available elementasthepivotisaverygoodchoice. Acuriosityofthisprocedure,however,is thatthechoiceofpivotwilldependontheoriginalscalingoftheequations. Ifwetake thethirdlinearequationin ouroriginalset andmultiplyit bya factorofa million,it 30 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).is almostguaranteedthatit willcontributethefirstpivot;yettheunderlyingsolution oftheequationsisnotchangedbythismultiplication! Onethereforesometimesseesroutines which choose as pivot that element which wouldhave been largest if the original equationshad all been scaled to have their largest coefficient normalizedto unity. Thisiscalled implicitpivoting . Thereissomeextrabookkeepingtokeeptrack of the scale factors by whichthe rows wouldhave beenmultiplied. (Theroutinesin §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 reflectionyouwill see that at everystage ofthe algorithm, eitheran element of Ais predictablyaoneorzero(ifitisalreadyinapartofthematrixthathasbeenreduced toidentityform) or elsethe exactlycorrespondingelementofthematrixthatstarted as1ispredictablyaoneorzero(ifitsmatein Ahasnotbeenreducedtotheidentity form). Thereforethematrix 1doesnothaveto existas separatestorage: Thematrix inverse of Ais gradually built up in Aas the original Ais destroyed. Likewise, the solution vectors xcan gradually replace the right-hand side vectors band share the same storage, since after each column in Ais reduced, the corresponding row entry in the b’s is never again used. Here is the routine for Gauss-Jordan elimination with full pivoting: SUBROUTINE gaussj(a,n,np,b,m,mp) INTEGER m,mp,n,np,NMAX REAL a(np,np),b(np,mp)PARAMETER (NMAX=50) Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) above. a(1:n,1:n) is an input matrix stored in an array of physical dimensions npbynp.b(1:n,1:m) is an in- put matrix containing the mright-hand side vectors, stored in an array of physical dimensions npbymp. On output, a(1:n,1:n) is replaced by its matrix inverse, and b(1:n,1:m) is replaced by the corresponding set of solution vectors. Parameter: NMAX is the largest anticipated value of n. INTEGER i,icol,irow,j,k,l,ll,indxc(NMAX),indxr(NMAX), * ipiv(NMAX) The integer arrays ipiv ,indxr ,a n d indxc are used for bookkeeping on the pivoting. REAL big,dum,pivinv do11j=1,n ipiv(j)=0 enddo 11 do22i=1,n This is the main loop over the columns to be re- duced. big=0. do13j=1,n This is the outer loop of the search for a pivot ele- ment. if(ipiv(j).ne.1)then do12k=1,n if (ipiv(k).eq.0) then if (abs(a(j,k)).ge.big)then big=abs(a(j,k)) irow=jicol=k endif endif enddo 12 endif enddo 13 ipiv(icol)=ipiv(icol)+1 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: 2.1Gauss-JordanElimination 31Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).indxc(i) ,t h ec o l u m no ft h e 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) /negationslash= 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.ne.icol) then do14l=1,n dum=a(irow,l)a(irow,l)=a(icol,l) a(icol,l)=dum enddo 14 do15l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum enddo 15 endifindxr(i)=irow W ea r en o wr e a d yt od i v i d et h ep i v o tr o wb yt h ep i v o t element, located at irow andicol . indxc(i)=icol if (a(icol,icol).eq.0.) pause ’singular matrix in gaussj’ pivinv=1./a(icol,icol) a(icol,icol)=1.do 16l=1,n a(icol,l)=a(icol,l)*pivinv enddo 16 do17l=1,m b(icol,l)=b(icol,l)*pivinv enddo 17 do21ll=1,n Next, we reduce the rows... if(ll.ne.icol)then ...except for the pivot one, of course. dum=a(ll,icol) a(ll,icol)=0.do 18l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum enddo 18 do19l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum enddo 19 endif enddo 21 enddo 22 This is the end of the main loop over columns of the reduction. do24l=n,1,-1 It only remains to unscramble the solution in view of the column interchanges. We do this by in-terchanging pairs of columns in the reverse orderthat the permutation was built up.if(indxr(l).ne.indxc(l))then do 23k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum enddo 23 endif enddo 24 return And we are done. END Rowversus ColumnEliminationStrategies The above discussion can be amplified by a modest amount of formalism. Row operations on a matrix Acorrespond to pre- (that is, left-) multiplication by some simple 32 Chapter2. SolutionofLinearAlgebraicEquationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).matrixR. For example, the matrix Rwith components Rij=  1ifi=jand i/negationslash=2 ,4 1ifi=2,j=4 1ifi=4,j=2 0otherwise(2.1.5 ) effects the interchange of rows 2and 4. Gauss-Jordan elimination by row operations alone (including the possibility of partialpivoting) consists of a series of such left-multiplications, yielding successively A·x=b (···R3·R2·R1·A)·x= ···R3·R2·R1·b (1)·x= ···R3·R2·R1·b x= ···R3·R2·R1·b(2.1.6 ) The key point is that since the R’s build from right to left, the right-hand side is simply transformed at each stage from one vector to another. Column operations, on the other hand, correspond to post-, or right-, multiplications by simple matrices, call them C. The matrix in equation (2.1.5), if right-multiplied onto a matrixA,willinterchange A’ssecond and fourth columns. Eliminationby column operations involves (conceptually) inserting a column operator, and also its inverse, between the matrix Aand the unknown vector x: A·x=b A·C1·C−1 1 ·x=b A·C1·C2·C−1 2 ·C−1 1 ·x=b (A·C1·C2·C3···)···C−1 3 ·C−1 2 ·C−1 1 ·x=b (1)···C−1 3 ·C−1 2 ·C−1 1 ·x=b(2.1.7 ) which (peeling of the C−1’s one at a time) implies a solution x=C1·C2·C3···b (2.1.8 ) Notice the essential difference between equation (2.1.8) and equation (2.1.6). In the latter case, the C’s must be applied to bin thereverse order from that in which they become known. That is, they must all be stored along the way. This requirement greatly reducesthe usefulness of column operations, generally restricting them to simple permutations, for example in support of full pivoting. CITED REFERENCES AND FURTHER READING: Wilkinson,J.H.1965, TheAlgebraicEigenvalueProblem (NewYork:OxfordUniversityPress).[1] Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), Example 5.2, p. 282. Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York: McGraw-Hill), Program B-2, p. 298. Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.3–1. 2.2GaussianEliminationwithBacksubstitution 33Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).2.2 GaussianEliminationwithBacksubstitution The usefulness of Gaussian elimination with backsubstitution is primarily pedagogical. It stands between full elimination schemes such as Gauss-Jordan, and triangular decomposition schemes such as will be discussed in the next section. Gaussian elimination reduces a matrix not all the way to the identity matrix, but onlyhalfway,toamatrixwhosecomponentsonthediagonalandabove(say)remain nontrivial. Let us now see what advantages accrue. Suppose that in doing Gauss-Jordan elimination, as described in §2.1, we at eachstagesubtractawayrowsonly belowthethen-currentpivotelement. When a22 is the pivot element, forexample,we dividethe secondrow by its value (as before), but now use the pivot row to zero only a32and a42, not a12(see equation 2.1.1). Suppose,also,thatwedoonlypartialpivoting,neverinterchangingcolumns,sothat the order of the unknowns never needs to be modified. Then, when we have done this for all the pivots, we will be left with a reduced equation that looks like this (in the case of a single right-handside vector):  a/prime 11 a/prime12 a/prime13 a/prime14 0 a/prime 22 a/prime23 a/prime24 00 a/prime 33 a/prime34 000 a/prime 44 · x1 x2 x3 x4 = b/prime 1 b/prime2 b/prime 3 b/prime 4  (2.2.1 ) Here the primes signify that the a’s and b’s do not have their original numerical values, but have been modified by all the row operations in the elimination to this point. The procedure up to this point is termed Gaussian elimination . Backsubstitution But how do we solve for the x’s? The last x(x4in this example) is already isolated, namely x4=b/prime 4/a/prime44(2.2.2 ) With the last xknown we can move to the penultimate x, x3=1 a/prime 33[b/prime 3−x4a/prime 34]( 2.2.3 ) and then proceed with the xbefore that one. The typical step is xi=1 a/prime ii b/prime i−N/summationdisplay j=i+1a/prime ijxj  (2.2.4 ) The procedure defined by equation (2.2.4) is called backsubstitution . The com- bination of Gaussian elimination and backsubstitution yields a solution to the set of equations.