f2-3

34 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).TheadvantageofGaussianeliminationandbacksubstitutionoverGauss-Jordan elimination is simply that the former is faster in raw operations count: Theinnermost loops of Gauss-Jordan elimination, each containing one subtraction and one multiplication, are executed N 3andN2Mtimes (where there are Nequations andMunknowns). The corresponding loops in Gaussian elimination are executed only1 3N3times (only half the matrix is reduced, and the increasing numbers of predictable zeros reduce the count to one-third), and1 2N2Mtimes, respectively. Each backsubstitutionof a right-handside is1 2N2executionsof a similar loop (one multiplication plus one subtraction). For M/lessmuchN(only a few right-hand sides) Gaussian elimination thus has about a factor three advantage over Gauss-Jordan.(We couldreducethis advantagetoa factor1.5by notcomputingtheinversematrix as part of the Gauss-Jordan scheme.) For computing the inverse matrix (which we can view as the case of M =N right-hand sides, namely the Nunit vectors which are the columns of the identity matrix),Gaussianeliminationandbacksubstitutionatfirstglancerequire 1 3N3(matrix reduction) +1 2N3(right-hand side manipulations) +1 2N3(Nbacksubstitutions) =4 3N3loopexecutions,whichismorethanthe N3forGauss-Jordan. However,the unit vectors are quite special in containing all zeros except for one element. If this is taken intoaccount,the right-sidemanipulationscanbereducedto only1 6N3loop executions,and,for matrixinversion,thetwo methodshave identicalefficiencies. BothGaussianeliminationandGauss-Jordaneliminationsharethedisadvantage that all right-handsides must beknownin advance. The LUdecompositionmethod in the next section does not share that deficiency, and also has an equally small operations count, both for solution with any number of right-hand sides, and formatrix inversion. For this reason we will not implement the method of Gaussian elimination as a routine. CITED REFERENCES AND FURTHER READING: Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §9.3–1. 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: Addison- Wesley), §2.2.1. Westlake,J.R.1968, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). 2.3 LU Decomposition and Its Applications Suppose we are able to write the matrix Aas a product of two matrices, L·U=A (2.3.1 ) whereLislower triangular (has elements only on the diagonal and below) and U isupper triangular (has elements only on the diagonal and above). For the case of 2.3LUDecompositionandItsApplications 35Sample 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).a4×4matrixA, for example, equation (2.3.1) would look like this:  α11 000 α21α22 00 α31α32α33 0 α41α42α43α44 · β11β12β13β14 0β22β23β24 00 β33β34 000 β44 = a11a12a13a14 a21a22a23a24 a31a32a33a34 a41a42a43a44 (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 ysuch that L·y=b (2.3.4 ) andthensolving 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 (equation2.2.4). Thus, equation (2.3.4)can be solved byforward substitution as follows, y 1=b1 α11 yi=1 αii bi−i−1/summationdisplay j=1αijyj  i=2,3,...,N(2.3.6 ) while (2.3.5)canthenbe solvedby backsubstitution exactlyas inequations(2.2.2)– (2.2.4), xN=yN βNN xi=1 βii yi−N/summationdisplay j=i+1βijxj  i=N−1,N−2,..., 1(2.3.7 ) Equations (2.3.6) and (2.3.7) total (for each right-hand side b)N2executions of an inner loop containing one multiply and one add. If we have Nright-hand sides which are the unit column vectors (which is the case when we are inverting a matrix),thentakinginto accountthe leadingzerosreducesthe total executioncount of (2.3.6) from1 2N3to1 6N3, while (2.3.7) is unchanged at1 2N3. Notice that, once we have the LUdecomposition of A, we can solve with as manyright-handsides as we then care to, one at a time. This is a distinct advantage over the methods of §2.1 and §2.2. 36 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).Performingthe LU Decomposition How then can we solve for LandU,g i v e nA? 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 iorjis the smaller number. We have, in fact, the three cases, i<j : αi1β1j+αi2β2j+··· +αiiβij=aij (2.3.8 ) i=j: αi1β1j+αi2β2j+··· +αiiβjj=aij (2.3.9 ) i>j : αi1β1j+αi2β2j+··· +αijβjj=aij (2.3.10 ) Equations(2.3.8)–(2.3.10)total N2equationsforthe N2+Nunknown α’sand β’s(thediagonalbeingrepresentedtwice). Sincethenumberofunknownsisgreater thanthenumberofequations,weareinvitedtospecify Noftheunknownsarbitrarily andthentrytosolvefortheothers. Infact,asweshallsee,itisalwayspossibletotake αii≡1 i=1,...,N (2.3.11 ) A surprising procedure, now, is Crout’s algorithm , which quite trivially solves thesetof N2+Nequations(2.3.8)–(2.3.11)forallthe α’sand β’sbyjustarranging 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,...,Ndo 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/summationdisplay k=1αikβkj. (2.3.12 ) (When i=1in2.3.12thesummationtermistakentomeanzero.) Second, fori=j+1,j+2,...,Nuse (2.3.10)to solve for αij, namely αij=1 βjj/parenleftBigg aij−j−1/summationdisplay k=1αikβkj/parenrightBigg . (2.3.13 ) Be sure to do both procedures before going on to the next j. 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 alreadydeterminedbythe time theyareneeded. You will also see that every aij is usedonlyonceandneveragain. Thismeansthatthecorresponding αijorβijcan be stored in the location that the aused to occupy: the decompositionis “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). 2.3LUDecompositionandItsApplications 37Sample 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).c g i b d f h jdiagonal elements subdiagonal elements etc.etc. xxa e Figure 2.3.1. Crout ’s algorithm for LUdecomposition 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”. What aboutpivoting? Pivoting (i.e., selection of a salubriouspivot element for the division in equation 2.3.13) is absolutely essential for the stability of Crout ’s method. Onlypartialpivoting(interchangeofrows)can beimplementedef ficiently. However this is enough to make the method stable. This means, incidentally, that we don’t actually decompose the matrix AintoLUform, 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(itsfinal application) is exactly the same as equation (2.3.13) except for the division in the latter equation; in both cases theupper limit of the sum is k=j−1( = i−1). This means that we don ’th a v et o commit ourselves as to whether the diagonal element β jjis the one that happens to fall on the diagonal in the first instance, or whether one of the (undivided) αij’s belowitinthecolumn, i=j+1,...,N,istobe“promoted ”tobecomethediagonal β. 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 38 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).method with partial pivoting . Our implementation has one additional wrinkle: It initiallyfinds the largest element in each row, and subsequently (when it is looking forthemaximalpivotelement)scalesthecomparison asifwehadinitiallyscaledall the equations to make their maximum coef ficient equal to unity; this is the implicit pivotingmentioned in §2.1. SUBROUTINE ludcmp(a,n,np,indx,d) INTEGER n,np,indx(n),NMAX REAL d,a(np,np),TINY PARAMETER (NMAX=500,TINY=1.0e-20) Largest expected n, and a small number. Given a matrix a(1:n,1:n) , with physical dimension npbynp, this routine replaces it by theLUdecomposition of a rowwise permutation of itself. aandnare input. ais 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; dis output as ±1depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksbto solve linear equations or invert a matrix. INTEGER i,imax,j,kREAL aamax,dum,sum,vv(NMAX) vv stores the implicit scaling of each row. d=1. No row interchanges yet. do 12i=1,n Loop over rows to get the implicit scaling informa- tion. aamax=0. do11j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) enddo 11 if (aamax.eq.0.) pause ’singular matrix in ludcmp’ No nonzero largest element. vv(i)=1./aamax Save the scaling. enddo 12 do19j=1,n This is the loop over columns of Crout’s method. do14i=1,j-1 This is equation (2.3.12) except for i=j. sum=a(i,j) do13k=1,i-1 sum=sum-a(i,k)*a(k,j) enddo 13 a(i,j)=sum enddo 14 aamax=0. Initialize for the search for largest pivot element. do16i=j,n This is i=jof equation (2.3.12) and i=j+1...N of equation (2.3.13). sum=a(i,j) do15k=1,j-1 sum=sum-a(i,k)*a(k,j) enddo 15 a(i,j)=sum dum=vv(i)*abs(sum) Figure of merit for the pivot. if (dum.ge.aamax) then Is it better than the best so far? imax=i aamax=dum endif enddo 16 if (j.ne.imax)then Do we need to interchange rows? do17k=1,n Yes, do so... dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum enddo 17 d=-d ...and change the parity of d. vv(imax)=vv(j) Also interchange the scale factor. endifindx(j)=imax if(a(j,j).eq.0.)a(j,j)=TINY If the pivot element is zero the matrix is singular (at least to the precision of the al-gorithm). For some applications on singular matrices, it is desirable to substitute TINY for zero. 2.3LUDecompositionandItsApplications 39Sample 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).if(j.ne.n)then Now, finally, divide by the pivot element. dum=1./a(j,j) do18i=j+1,n a(i,j)=a(i,j)*dum enddo 18 endif enddo 19 Go back for the next column in the reduction. returnEND Hereistheroutineforforwardsubstitutionandbacksubstitution,implementing equations (2.3.6) and (2.3.7). SUBROUTINE lubksb(a,n,np,indx,b) INTEGER n,np,indx(n)REAL a(np,np),b(n) Solves the set of nlinear equations A·X=B.H e r e ais input, not as the matrix Abut rather as its LUdecomposition, determined by the routine ludcmp.indxis input as the permutation vector returned by ludcmp.b(1:n)is input as the right-hand side vector B, and returns with the solution vector X.a,n,np,a n d indxare not modified by this routine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that bwill begin with many zero elements, so it is efficient for use in matrix inversion. INTEGER i,ii,j,ll REAL sum ii=0 When iiis set to a positive value, it will become the in- dex of the first nonvanishing element of b.W en o wd o the forward substitution, equation (2.3.6). The only new wrinkle is to unscramble the permutation as we go.do12i=1,n ll=indx(i) sum=b(ll)b(ll)=b(i)if (ii.ne.0)then do 11j=ii,i-1 sum=sum-a(i,j)*b(j) enddo 11 else if (sum.ne.0.) then ii=i A nonzero element was encountered, so from now on we will have to do the sums in the loop above. endif b(i)=sum enddo 12 do14i=n,1,-1 Now we do the backsubstitution, equation (2.3.7). sum=b(i) do13j=i+1,n sum=sum-a(i,j)*b(j) enddo 13 b(i)=sum/a(i,i) Store a component of the solution vector X. enddo 14 return All done! END TheLUdecompositionin ludcmprequires about1 3N3executionsof 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 gaussjwhich 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 asgaussj. 40 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).To summarize, this is the preferred way to solve the linear set of equations A·x=b: call ludcmp(a,n,np,indx,d) call lubksb(a,n,np,indx,b) The answer xwill be returned in b. Your original matrix Awill have been destroyed. If you subsequently want to solve a set of equations with the same Abut a different right-hand side b, you repeat only call lubksb(a,n,np,indx,b) not, of course, with the original matrix A, but with aandindxas were already returned from ludcmp. Inverse of a Matrix Using the above LUdecomposition and backsubstitution routines, it is com- pletely straightforwardto find the inverse of a matrix column by column. INTEGER np,indx(np) REAL a(np,np),y(np,np)... do 12i=1,n Set up identity matrix. do11j=1,n y(i,j)=0. enddo 11 y(i,i)=1. enddo 12 call ludcmp(a,n,np,indx,d) Decompose the matrix just once. do13j=1,n Find inverse by columns. call lubksb(a,n,np,indx,y(1,j)) Note that FORTRAN stores two-dimensional matrices by column, so y(1,j)is the address of the jth column of y. enddo 13 The matrix ywill now contain the inverse of the original matrix a, which will have been destroyed. Alternatively, there is nothing wrong with using a Gauss-Jordanroutine like gaussj(§2.1) to inverta matrix in place, again destroyingthe original. Both methods have practically the same operations count. Incidentally, if you ever have the need to compute A −1·Bfrom matrices A andB, you should LUdecompose Aand then backsubstitute with the columns of Binstead of with the unit vectors that would give A’s inverse. This saves a whole matrix multiplication, and is also more accurate. 2.3LU DecompositionandIts Applications 41Sample 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).Determinant of a Matrix The determinant of an LUdecomposed matrix is just the product of the diagonal elements, det =N/productdisplay 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 ofwhether the number of row interchanges was even or odd, so we just preface the productby the correspondingsign. (You now finally know the purpose of returning din the routine ludcmp.) Calculation of a determinant thus requires one call to ludcmp, withnosubse- quent backsubstitutions by lubksb. INTEGER np,indx(np) REAL a(np,np)... call ludcmp(a,n,np,indx,d) This returns das±1. do 11j=1,n d=d*a(j,j) enddo 11 The variable dnow 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 under flow your computer ’sfloating-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. ComplexSystems of Equations If your matrix Ais real, but the right-hand side vector is complex, say b+id, then (i) LUdecompose Ain the usual way, (ii) backsubstitute bto get the real part of the solution vector, and (iii) backsubstitute dto 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 ludcmpandlubksb as complex routines. Complex modulus substitutes for absolute value in the construction ofthe scaling vector vvand in the search for the largest pivot elements. Everything else goes through in the obvious way, with complex arithmetic used as needed. 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×2Nset ofrealequations, /parenleftbigg A−C CA/parenrightbigg ·/parenleftbigg x y/parenrightbigg =/parenleftbigg b d/parenrightbigg (2.3.18 ) 42 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).and then solved with ludcmpandlubksbin their present forms. This scheme is a factor of 2 inefficient in storage, since AandCare stored twice. It is also a factor of 2 inef ficient in time, since the complex multiplies in a complexi fied version of the routines would each use 4 real multiplies, while the solution of a 2N×2Nproblem involves 8 times the work of anN×None. If you can tolerate these factor-of-two inef ficiencies, then equation (2.3.18) is an easy way to proceed. CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins 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, AHandbookofNumericalMatrixInversionandSolutionofLinearEquations (New York: Wiley). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork: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.,andJohnson,C.R.1985, MatrixAnalysis (Cambridge:CambridgeUniversityPress). 2.4 Tridiagonal and Band Diagonal Systems of Equations The special case of a system of linear equations that is tridiagonal , that is, has nonzeroelementsonlyonthediagonalplusorminusonecolumn,is onethat occursfrequently. Alsocommonaresystemsthatare banddiagonal ,withnonzeroelements onlyalonga fewdiagonallines adjacent to the maindiagonal(aboveandbelow). For tridiagonal sets, the procedures of LUdecomposition, forward- and back- substitutioneachtakeonly O(N)operations,andthewholesolutioncanbeencoded veryconcisely. Theresultingroutine tridagisonethatwewilluseinlaterchapters. Naturally,one does not reservestorage forthe full N×Nmatrix, but only for thenonzerocomponents,storedasthreevectors. Thesetofequationstobesolvedis  b 1c1 0··· a2b2c2······ ··· a N−1bN−1cN−1 ··· 0 aN bN · u 1 u2 ··· uN−1 uN = r 1 r2 ··· rN−1 rN (2.4.1 ) Noticethat a 1andcNareundefinedandarenotreferencedbytheroutinethatfollows.