f2-7

2.7SparseLinearSystems 63Sample 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).REAL absa,absb absa=abs(a) absb=abs(b) if(absa.gt.absb)then pythag=absa*sqrt(1.+(absb/absa)**2) else if(absb.eq.0.)then pythag=0. else pythag=absb*sqrt(1.+(absa/absb)**2) endif endifreturn END (Double precision versions of svdcmp,svbksb, and pythag, named dsvdcmp, dsvbksb, and dpythag, are used by the routine ratlsqin§5.13. You can easily maketheconversions,orelsegettheconvertedroutinesfromthe NumericalRecipes diskette.) CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §8.3 and Chapter 12. Lawson, C.L., and Hanson, R. 1974, Solving Least Squares Problems (Englewood Cliffs, NJ: Prentice-Hall), Chapter 18. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), Chapter 9. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(NewYork: Springer-Verlag), Chapter I.10 by G.H. Golub andC. Reinsch. [2] Dongarra, J.J., et al. 1979, LINPACK User’s Guide (Philadelphia: S.I.A.M.), Chapter 11. [3] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §6.7. [4] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §5.2.6. [5] 2.7 Sparse Linear Systems A system of linear equations is called sparseif only a relatively small number of its matrix elements aijare nonzero. It is wasteful to use general methods of linear algebra on such problems, because most of the O(N3)arithmetic operations devotedtosolvingthesetofequationsorinvertingthematrixinvolvezerooperands. Furthermore, you might wish to work problems so large as to tax your available memory space, and it is wasteful to reserve storage for unfruitful zero elements.Note that there are two distinct (and not always compatible) goals for any sparse matrix method: saving time and/or saving space. We have already considered one archetypal sparse form in §2.4, the band diagonal matrix. In the tridiagonal case, e.g., we saw that it was possible to save 64 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).both time (order Ninstead of N3) and space (order Ninstead of N2). The method of solution was not different in principle from the general method of LU decomposition;itwasjustappliedcleverly,andwithdueattentiontothebookkeeping ofzeroelements. Manypracticalschemesfordealingwithsparseproblemshavethis samecharacter. Theyarefundamentallydecompositionschemes,orelseeliminationschemesakintoGauss-Jordan,butcarefullyoptimizedsoastominimizethenumber of so-called fill-ins, initially zero elements which must become nonzero during the solution process, and for which storage must be reserved. Direct methods for solving sparse equations, then, depend crucially on the precise pattern of sparsity of the matrix. Patterns that occur frequently, or that areuseful as way-stations in the reduction of more general forms, already have special names and special methods of solution. We do not have space here for any detailed review ofthese. Referenceslisted at the end of this sectionwill furnishyou with an“in” to the specialized literature, and the following list of buzz words (and Figure 2.7.1) will at least let you hold your own at cocktail parties: •tridiagonal •band diagonal (or banded) with bandwidth M •band triangular •block diagonal •block tridiagonal •block triangular •cyclic banded •singly (or doubly) bordered block diagonal •singly (or doubly) bordered block triangular •singly (or doubly) bordered band diagonal •singly (or doubly) bordered band triangular •other (!) You should also be aware of some of the special sparse forms that occur in the solutionofpartialdifferentialequationsintwoormoredimensions. SeeChapter19. If your particular pattern of sparsity is not a simple one, then you may wish to tryananalyze/factorize/operate package,whichautomatestheprocedureoffiguring out how fill-ins are to be minimized. The analyzestage is done once only for each pattern of sparsity. The factorize stage is done once for each particular matrix that fits the pattern. The operatestage is performed once for each right-hand side to be used with the particular matrix. Consult [2,3]for references on this. The NAG library[4]has an analyze/factorize/operate capability. A substantial collection of routines for sparse matrix calculation is also available from IMSL [5]as theYale Sparse Matrix Package [6]. You should be aware that the special order of interchanges and eliminations, prescribed by a sparse matrix method so as to minimize fill-ins and arithmeticoperations, generally acts to decrease the method’s numerical stability as compared to, e.g., regular LUdecomposition with pivoting. Scaling your problem so as to make its nonzero matrix elements have comparable magnitudes (if you can do it)will sometimes ameliorate this problem. Intheremainderofthissection,wepresentsomeconceptswhichareapplicable to some general classes of sparse matrices, and which do not necessarily depend on details of the pattern of sparsity. 2.7SparseLinearSystems 65Sample 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).(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k)zeroszeros zeros Figure2.7.1. Somestandard formsforsparsematrices. (a)Banddiagonal; (b)block triangular; (c)block tridiagonal; (d) singly bordered block diagonal; (e) doubly bordered block diagonal; (f) singly bordered block triangular; (g) bordered band-triangular; (h) and (i) singly and doubly bordered band diagonal; (j)and (k) other! (after Tewarson) [1]. Sherman-MorrisonFormula Supposethatyouhavealreadyobtained,byherculeaneffort,theinversematrix A−1of a square matrix A. Now you want to make a “small”change in A, for example change one element aij, or a few elements, or one row, or one column. Is there any way of calculating the correspondingchange in A−1without repeating 66 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).your difficult labors? Yes, if your change is of the form A→ (A+u⊗v)( 2.7.1 ) for some vectors uandv.I fuis a unit vector ei, then (2.7.1) adds the components ofvtothe ithrow. (Recallthat u⊗vis amatrixwhose i, jthelementis theproduct of the ith componentof uandthe jth componentof v.) Ifvis a unit vector ej, then (2.7.1)addsthecomponentsof utothe jthcolumn. Ifboth uandvareproportional to unitvectors eiandejrespectively,thena term is addedonlyto the element aij. TheSherman-Morrison formulagivestheinverse (A+u⊗v)−1,andisderived briefly as follows: (A+u⊗v)−1=(1+A−1·u⊗v)−1·A−1 =(1−A−1·u⊗v+A−1·u⊗v·A−1·u⊗v−...)·A−1 =A−1−A−1·u⊗v·A−1(1−λ+λ2−...) =A−1−(A−1·u)⊗(v·A−1) 1+λ (2.7.2 ) where λ≡v·A−1·u (2.7.3 ) The second line of (2.7.2) is a formal power series expansion. In the third line, the associativity of outer and inner products is used to factor out the scalars λ. The use of (2.7.2) is this: Given A−1and the vectors uandv, we need only perform two matrix multiplications and a vector dot product, z≡A−1·uw ≡(A−1)T·v λ=v·z (2.7.4 ) to get the desired change in the inverse A−1→A−1−z⊗w 1+λ(2.7.5 ) The whole procedure requires only 3N2multiplies and a like number of adds (an even smaller number if uorvis a unit vector). The Sherman-Morrison formula can be directly applied to a class of sparse problems. If you already have a fast way of calculating the inverse of A(e.g., a tridiagonal matrix, or some other standard sparse form), then (2.7.4) –(2.7.5) allow you to build up to your related but more complicated form, adding for example a row or column at a time. Notice that you can apply the Sherman-Morrisonformulamore than once successively, using at each stage the most recent update of A −1 (equation 2.7.5). Of course, if you have to modify everyrow, then you are back to anN3method. The constant in front of the N3is only a few times worse than the better direct methods, but you have deprived yourself of the stabilizing advantages of pivoting —so be careful. For some other sparse problems, the Sherman-Morrison formula cannot be directly applied for the simple reason that storage of the whole inverse matrix A−1 2.7SparseLinearSystems 67Sample 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 not feasible. If you want to add only a single correction of the form u⊗v, and solve the linear system (A+u⊗v)·x=b (2.7.6 ) then you proceed as follows. Using the fast method that is presumed available for the matrix A, solve the two auxiliary problems A·y=bA ·z=u (2.7.7 ) for the vectors yandz. In terms of these, x=y−/bracketleftbiggv·y 1+(v·z)/bracketrightbigg z (2.7.8 ) as we see by multiplying (2.7.2) on the right by b. Cyclic TridiagonalSystems So-called cyclic tridiagonal systems occur quite frequently, and are a good exampleof howtouse the Sherman-Morrisonformulain themannerjust described. The equations have the form  b 1c1 0··· β a2b2c2··· ··· ··· aN−1bN−1cN−1 α ··· 0 aN bN · x 1 x2 ··· xN−1 xN = r 1 r2 ··· rN−1 rN (2.7.9 ) This is a tridiagonal system, except for the matrix elements αandβin the corners. Forms like this are typically generated by finite-differencing differential equations with periodic boundary conditions ( §19.4). We use the Sherman-Morrisonformula, treating the system as tridiagonal plus a correction. In the notation of equation (2.7.6),de fine vectors uandvto be u= γ 0... 0 α v= 1 0... 0 β/γ (2.7.10 ) Here γis arbitraryfor the moment. Then the matrix Ais the tridiagonal part of the matrix in (2.7.9), with two terms modi fied: b /prime 1=b1−γ, b/prime N=bN−αβ/γ (2.7.11 ) We now solve equations (2.7.7) with the standard tridiagonal algorithm, and then get the solution from equation (2.7.8). Theroutine cyclicbelowimplementsthis algorithm. We choosethe arbitrary parameter γ=−b1toavoid loss ofprecisionbysubtractionin the first ofequations (2.7.11). In the unlikely event that this causes loss of precision in the second of these equations, you can make a different choice. 68 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).SUBROUTINE cyclic(a,b,c,alpha,beta,r,x,n) INTEGER n,NMAX REAL alpha,beta,a(n),b(n),c(n),r(n),x(n) PARAMETER (NMAX=500) C USES tridag Solvesforavector x(1:n)the “cyclic”setoflinearequations givenbyequation (2.7.9). a,b,c,and rareinputvectors,while alphaandbetaarethecornerentriesinthematrix. The input is not modified. INTEGER i REAL fact,gamma,bb(NMAX),u(NMAX),z(NMAX) if(n.le.2)pause ’n too small in cyclic’if(n.gt.NMAX)pause ’NMAX too small in cyclic’gamma=-b(1) Avoidsubtractionerrorinforming bb(1). bb(1)=b(1)-gamma Setupthediagonalofthemodifiedtridiagonalsystem. bb(n)=b(n)-alpha*beta/gammado 11i=2,n-1 bb(i)=b(i) enddo 11 call tridag(a,bb,c,r,x,n) SolveA·x=r. u(1)=gamma Set up the vector u. u(n)=alpha do12i=2,n-1 u(i)=0. enddo 12 call tridag(a,bb,c,u,z,n) SolveA·z=u. fact=(x(1)+beta*x(n)/gamma)/(1.+z(1)+beta*z(n)/gamma) Formv·x/(1 +v·z). do13i=1,n Nowgetthesolution vector x. x(i)=x(i)-fact*z(i) enddo 13 return END WoodburyFormula If you want to add more than a single correction term, then you cannot use (2.7.8) repeatedly, since without storing a new A−1you will not be able to solve the auxiliary problems(2.7.7)ef ficientlyafterthe firststep. Instead,youneedthe Woodburyformula ,which is the block-matrix version of the Sherman-Morrison formula, (A+U·VT)−1 =A−1−/bracketleftBig A−1·U·(1+VT·A−1·U)−1·VT·A−1/bracketrightBig (2.7.12 ) Here Ais, as usual, an N×Nmatrix, while UandVareN×Pmatrices with P<N and usually P/lessmuchN. The inner piece of the correction term may become clearer if written as the tableau,  U · 1+V T·A−1·U −1 · VT (2.7.13 ) where you cansee thatthematrix whose inverseisneeded isonly P×Pratherthan N×N. 2.7SparseLinearSystems 69Sample 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).TherelationbetweentheWoodburyformulaandsuccessiveapplicationsoftheSherman- Morrisonformulaisnowclari fiedbynotingthat,if Uisthematrixformedbycolumnsoutofthe Pvectorsu1,...,uP,andVisthematrixformedbycolumnsoutofthe Pvectorsv1,...,vP, U≡ u 1 ··· u P V≡ v 1 ··· v P (2.7.14 ) then two ways of expressing the same correction to Aare /parenleftBigg A+P/summationdisplay k=1uk⊗vk/parenrightBigg =(A+U·VT)( 2.7.15 ) (Note that the subscripts on uandvdonotdenote components, but rather distinguish the different column vectors.) Equation (2.7.15) reveals that, if you have A−1in storage, then you can either make the Pcorrections in one fell swoop by using (2.7.12), inverting a P×Pmatrix, or else make them by applying (2.7.5) Psuccessive times. If you don ’t have storage for A−1, then you mustuse (2.7.12) in the following way: To solve the linear equation /parenleftBigg A+P/summationdisplay k=1uk⊗vk/parenrightBigg ·x=b (2.7.16 ) first solve the Pauxiliary problems A·z1=u1 A·z2=u2 ··· A·zP=uP(2.7.17 ) and construct the matrix Zby columns from the z’s obtained, Z≡ z 1 ··· z P (2.7.18 ) Next, do the P×Pmatrix inversion H≡(1+VT·Z)−1(2.7.19 ) Finally, solve the one further auxiliary problem A·y=b (2.7.20 ) In terms of these quantities, the solution is given by x=y−Z·/bracketleftBig H·(VT·y)/bracketrightBig (2.7.21 ) 70 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).Inversionby Partitioning Once in a while, you will encounter a matrix (not even necessarily sparse) that can be inverted ef ficiently by partitioning. Suppose that the N×Nmatrix Ais partitioned into A=/bracketleftbigg PQ RS/bracketrightbigg (2.7.22 ) wherePandSaresquarematricesofsize p×pands×srespectively( p+s=N). The matrices QandRare not necessarily square, and have sizes p×sands×p, respectively. If the inverse of Ais partitioned in the same manner, A−1=/bracketleftBigg/tildewideP/tildewideQ /tildewideR/tildewideS/bracketrightBigg (2.7.23 ) then/tildewideP,/tildewideQ,/tildewideR,/tildewideS, which have the same sizes as P,Q,R,S, respectively, can be found by either the formulas /tildewideP=(P−Q·S−1·R)−1 /tildewideQ=−(P−Q·S−1·R)−1·(Q·S−1) /tildewideR=−(S−1·R)·(P−Q·S−1·R)−1 /tildewideS=S−1+(S−1·R)·(P−Q·S−1·R)−1·(Q·S−1)(2.7.24 ) or else by the equivalent formulas /tildewideP=P−1+(P−1·Q)·(S−R·P−1·Q)−1·(R·P−1) /tildewideQ=−(P−1·Q)·(S−R·P−1·Q)−1 /tildewideR=−(S−R·P−1·Q)−1·(R·P−1) /tildewideS=(S−R·P−1·Q)−1(2.7.25 ) The parentheses in equations (2.7.24) and (2.7.25) highlight repeated factors that you may wish to compute only once. (Of course, by associativity, you can instead do the matrix multiplications in any order you like.) The choice between usingequation (2.7.24) and (2.7.25) depends on whether you want /tildewidePor/tildewideSto have the simplerformula;oronwhethertherepeatedexpression (S−R·P −1·Q)−1iseasier to calculate than the expression (P−Q·S−1·R)−1; or on the relative sizes of P andS; or on whether P−1orS−1is already known. Another sometimes useful formula is for the determinant of the partitioned matrix, detA=d e tPdet(S−R·P−1·Q)=d e tSdet(P−Q·S−1·R)(2.7.26 ) 2.7SparseLinearSystems 71Sample 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).Indexed Storage of Sparse Matrices Wehavealreadyseen( §2.4)thattri-orband-diagonalmatricescanbestoredinacompact formatthatallocatesstorageonlytoelementswhichcanbenonzero,plusperhapsafewwastedlocations tomakethe bookkeeping easier. Whatabout moregeneral sparsematrices? Whenasparse matrixoflogical size N×Ncontains only a fewtimes Nnonzero elements (atypical case), it is surely inef ficient—and often physically impossible —to allocate storage for all N 2elements. Even if one did allocate such storage, it would be inef ficient or prohibitive in machine time to loop over all of it in search of nonzero elements. Obviouslysomekindofindexedstorageschemeisrequired,onethatstoresonlynonzero matrix elements, along with suf ficient auxiliary information to determine where an element logically belongs and how the various elements can be looped over in common matrixoperations. Unfortunately,thereisnoonestandardschemeingeneraluse. Knuth [7]describes one method. The Yale Sparse Matrix Package [6]and ITPACK [8]describe several other methods. For most applications, we favor the storage scheme used by PCGPACK [9], which isalmostthesameasthatdescribedbyBentley [10],andalsosimilartooneoftheYaleSparse Matrix Package methods. The advantage of this scheme, which can be called row-indexed sparse sto ragemode,isthatitrequiresstorageofonlyabout twotimesthenumberofnonzero matrix elements. (Other methods can require as much as three or five times.) For simplicity, we will treat only the case of square matrices, which occurs most frequently in practice. To represent a matrix Aof logical size N×N, the row-indexed scheme sets up two one-dimensional arrays,call them saandija. Thefirstof these stores matrixelement values insingleordoubleprecisionasdesired;thesecondstoresintegervalues. Thestoragerulesare: •ThefirstNlocationsof sastoreA’sdiagonalmatrixelements,inorder. (Notethat diagonal elements are stored even if they are zero; this is at most a slight storageinefficiency, since diagonal elements are nonzero in most realistic applications.) •Each of the firstNlocations of ijastores the index of the array sathat contains thefirstoff-diagonal element of the corresponding row of the matrix. (If there are no off-diagonal elements for that row, it is one greater than the index in saof the most recently stored element of a previous row.) •Location 1 of ijais always equal to N+2. (It can be read to determine N.) •Location N+1ofijais one greater than the index in saof the last off-diagonal element of the last row. (It can be read to determine the number of nonzeroelements in the matrix, or the logical length of the arrays saandija.) Location N+1ofsais not used and can be set arbitrarily. •Entries in saat locations ≥N+2containA’s off-diagonal values, ordered by rows and, within each row, ordered by columns. •Entriesin ijaatlocations ≥N+2containthecolumnnumberofthecorresponding element in sa. While these rules seem arbitrary at first sight, they result in a rather elegant storage scheme. As an example, consider the matrix  3.0.1.0.0. 0.4.0.0.0. 0.7.5.9.0. 0.0.0.0.2. 0.0.0.6.5. (2.7.27 ) In row-indexed compact storage, matrix (2.7.27) is represented by the two arrays of length 11, as follows index k 1 2 3 4 5 6 7 8 910 11 ija(k) 7 8 810 11 12 3 2 4 5 4 sa(k) 3.4.5.0.5. x 1.7.9.2.6.(2.7.28 ) Here xis an arbitrary value. Notice that, according to the storage rules, the value of N (namely 5) is ija(1)-2 , and the length of each array is ija(ija(1)-1)-1 , namely 11. 72 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).The diagonal element in row iissa(i), and the off-diagonal elements in that row are in sa(k)where kloops from ija(i)toija(i+1)-1 , if the upper limit is greater or equal to the lower one (as in FORTRAN do loops). Hereisaroutine, sprsin,thatconvertsamatrixfromfullstoragemodeintorow-indexed sparse storage mode, throwing away any elements that are less than a speci fied threshold. Of course, the principal use of sparse storage mode is for matrices whose full storage modewon’tfitinto your machine at all; then you have to generate them directly into sparse format. Nevertheless sprsinis useful as a precise algorithmic de finition of the storage scheme, for subscale testing oflarge problems, and forthe case whereexecution time,ratherthan storage,furnishes the impetus to sparse storage. SUBROUTINE sprsin(a,n,np,thresh,nmax,sa,ija) INTEGER n,nmax,np,ija(nmax)REAL thresh,a(np,np),sa(nmax) Convertsasquarematrix a(1:n,1:n)withphysicaldimension npintorow-indexedsparse storage mode. Only elements of awith magnitude ≥threshare retained. Output is in two linear arrays with physical dimension nmax(an input parameter): sa(1:)contains array values, indexed by ija(1:). The logical sizes of saandijaon output are both ija(ija(1)-1)-1 (see text). INTEGER i,j,kdo 11j=1,n Storediagonal elements. sa(j)=a(j,j) enddo 11 ija(1)=n+2 Indexto1strowoff-diagonalelement,ifany. k=n+1 do13i=1,n Loop over rows. do12j=1,n Loop over columns. if(abs(a(i,j)).ge.thresh)then if(i.ne.j)then Storeoff-diagonalelementsandtheircolumns. k=k+1 if(k.gt.nmax)pause ’nmax too small in sprsin’sa(k)=a(i,j) ija(k)=j endif endif enddo 12 ija(i+1)=k+1 Aseachrowiscompleted,storeindextonext. enddo 13 return END The single most important use of a matrix in row-indexed sparse storage mode is to multiply a vector to its right. In fact, the storage mode is optimized for just this purpose.The following routine is thus very simple. SUBROUTINE sprsax(sa,ija,x,b,n) INTEGER n,ija(*) REAL b(n),sa(*),x(n) Multiplyamatrixinrow-indexsparsestoragearrays saandijabyavector x(1:n),giving a vector b(1:n). INTEGER i,k if (ija(1).ne.n+2) pause ’mismatched vector and matrix in sprsax’ do12i=1,n b(i)=sa(i)*x(i) Startwithdiagonalterm. do11k=ija(i),ija(i+1)-1 Loopoveroff-diagonalterms. b(i)=b(i)+sa(k)*x(ija(k)) enddo 11 enddo 12 returnEND 2.7SparseLinearSystems 73Sample 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).Itisalsosimpletomultiplythe transpose ofamatrixbyavectortoitsright. (Wewilluse this operation laterin this section.) Note that the transpose matrixis not actually constructed. SUBROUTINE sprstx(sa,ija,x,b,n) INTEGER n,ija(*) REAL b(n),sa(*),x(n) Multiply the transpose of a matrix in row-index sparse storage arrays saandijaby a vector x(1:n), giving a vector b(1:n). INTEGER i,j,k if (ija(1).ne.n+2) pause ’mismatched vector and matrix in sprstx’ do11i=1,n Start withdiagonal terms. b(i)=sa(i)*x(i) enddo 11 do13i=1,n Loop overoff-diagonal terms. do12k=ija(i),ija(i+1)-1 j=ija(k) b(j)=b(j)+sa(k)*x(i) enddo 12 enddo 13 returnEND (Double precision versions of sprsaxandsprstx, named dsprsax anddsprstx, are used by the routine atimeslater in this section. You can easily make the conversion, or else get the converted routines from the Numerical Recipes diskettes.) In fact, because the choice of row-indexed storage treats rows and columns quite differently, it is quite an involved operation to construct the transpose of a matrix, given the matrix itself in row-indexed sparse storage mode. When the operation cannot be avoided, it is done as follows: An index of all off-diagonal elements by their columns is constructed(see§8.4). The elements are then written to the output array in column order. As each element is written, its row is determined and stored. Finally, the elements in each columnare sorted by row. SUBROUTINE sprstp(sa,ija,sb,ijb) INTEGER ija(*),ijb(*)REAL sa(*),sb(*) C USES iindexx Versionof indexxwithall REALvariableschangedto INTEGER. Constructthetransposeofasparsesquarematrix,fromrow-indexsparsestoragearrays sa andijainto arrays sbandijb. INTEGER j,jl,jm,jp,ju,k,m,n2,noff,inc,iv REAL v n2=ija(1) Linear sizeofmatrixplus2. do11j=1,n2-2 Diagonal elements. sb(j)=sa(j) enddo 11 call iindexx(ija(n2-1)-ija(1),ija(n2),ijb(n2)) Index all off-diagonal elements by their columns. jp=0 do13k=ija(1),ija(n2-1)-1 Loopoveroutput off-diagonalelements. m=ijb(k)+n2-1 Useindextabletostoreby(former)columns. sb(k)=sa(m) do12j=jp+1,ija(m) Fillintheindextoanyomittedrows. ijb(j)=k enddo 12 jp=ija(m) Usebisectiontofindwhichrowelement misinandputthat into ijb(k). jl=1 ju=n2-1 5 if (ju-jl.gt.1) then jm=(ju+jl)/2 if(ija(jm).gt.m)then ju=jm else 74 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).jl=jm endif goto 5 endifijb(k)=jl enddo 13 do14j=jp+1,n2-1 ijb(j)=ija(n2-1) enddo 14 MakeafinalpasstosorteachrowbyShellsortalgorithm. do16j=1,n2-2 jl=ijb(j+1)-ijb(j)noff=ijb(j)-1inc=1 1 inc=3*inc+1 if(inc.le.jl)goto 1 2 continue inc=inc/3 do 15k=noff+inc+1,noff+jl iv=ijb(k)v=sb(k) m=k 3 if(ijb(m-inc).gt.iv)then ijb(m)=ijb(m-inc) sb(m)=sb(m-inc) m=m-incif(m-noff.le.inc)goto 4 goto 3 endif 4 ijb(m)=iv sb(m)=v enddo 15 if(inc.gt.1)goto 2 enddo 16 return END Theaboveroutineembedsinternallyasortingalgorithmfrom §8.1,butcallstheexternal routine iindexx to construct the initialcolumn index. Thisroutine isidentical to indexx,as listed in §8.4, except that the latter ’st w o REALdeclarations should be changed to integer. (TheNumerical Recipes diskettes include both indexxandiindexx.) In fact, you can often use indexxwithoutmaking these changes, since many computers have the property that numerical values will sort correctly independently of whether they are interpreted asfloating or integer values. Asfinal examples of the manipulation of sparse matrices, we give two routines for the multiplicationoftwosparsematrices. Theseareusefulfortechniquestobedescribedin §13.10. In general, the product of two sparse matrices is not itself sparse. One therefore wants tolimitthesizeoftheproduct matrixinoneoftwoways: eithercompute only thoseelementsoftheproduct thatarespeci fiedinadvance byaknown pattern ofsparsity,orelsecompute all nonzero elements, but store only those whose magnitude exceeds some threshold value. Theformer technique, when it can be used, is quite ef ficient. The pattern of sparsity is speci fied by furnishing an index array in row-index sparse storage format (e.g., ija). The program then constructs a corresponding value array(e.g., sa). The lattertechnique runs the danger of excessive compute times and unknown output sizes, so it must be used cautiously. With row-index storage, it is much more natural to multiply a matrix (on the left) by thetranspose of a matrix (on the right), so that one is crunching rows on rows, rather than rows on columns. Our routines therefore calculate A·B T, rather than A·B. This means that you have to run your right-hand matrix through the transpose routine sprstpbefore sending it to the matrix multiply routine. Thetwoimplementingroutines, sprspmfor“patternmultiply ”andsprstmfor“threshold multiply”are quite similar in structure. Both are complicated by the logic of the various 2.7SparseLinearSystems 75Sample 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).combinationsofdiagonaloroff-diagonalelementsforthetwoinputstreamsandoutputstream. SUBROUTINE sprspm(sa,ija,sb,ijb,sc,ijc) INTEGER ija(*),ijb(*),ijc(*) REAL sa(*),sb(*),sc(*) Matrixmultiply A·BTwhereAandBaretwosparsematricesinrow-indexstoragemode, andBTisthetransposeof B. Here, saandijastorethematrix A;sbandijbstorethe matrixB. Thisroutinecomputesonlythosecomponentsofthematrixproductthatare pre- specifiedbytheinputindexarray ijc,whichisnotmodified. Onoutput,thearrays scand ijcgivetheproduct matrixinrow-indexstoragemode. Forsparsematrixmultiplication, this routine willoften bepreceded byacall to sprstp, soas toconstruct thetranspose of a known matrix into sb,ijb. INTEGER i,ijma,ijmb,j,m,ma,mb,mbb,mn REAL sum if (ija(1).ne.ijb(1).or.ija(1).ne.ijc(1)) * pause ’sprspm sizes do not match’ do13i=1,ijc(1)-2 Loop over rows. j=i Setupsothatfirstpassthroughloopdoesthediag- onal component. m=i mn=ijc(i) sum=sa(i)*sb(i) 1 continue Mainloopovereachcomponenttobeoutput. mb=ijb(j)do 11ma=ija(i),ija(i+1)-1 Loopthroughelementsin A’srow. Convolutedlogic, following,accountsforthevariouscombinations ofdiagonalandoff-diagonalelements.ijma=ija(ma) if(ijma.eq.j)then sum=sum+sa(ma)*sb(j) else 2 if(mb.lt.ijb(j+1))then ijmb=ijb(mb)if(ijmb.eq.i)then sum=sum+sa(i)*sb(mb) mb=mb+1goto 2 else if(ijmb.lt.ijma)then mb=mb+1 goto 2 else if(ijmb.eq.ijma)then sum=sum+sa(ma)*sb(mb) mb=mb+1goto 2 endif endif endif enddo 11 do12mbb=mb,ijb(j+1)-1 Exhausttheremainderof B’srow. if(ijb(mbb).eq.i)then sum=sum+sa(i)*sb(mbb) endif enddo 12 sc(m)=sum sum=0.e0 Resetindicesfornextpassthroughloop. if(mn.ge.ijc(i+1))goto 3 m=mn mn=mn+1j=ijc(m) goto 1 3 continue enddo 13 return END 76 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).SUBROUTINE sprstm(sa,ija,sb,ijb,thresh,nmax,sc,ijc) INTEGER nmax,ija(*),ijb(*),ijc(nmax) REAL thresh,sa(*),sb(*),sc(nmax) Matrixmultiply A·BTwhereAandBaretwosparsematricesinrow-indexstoragemode, andBTisthetransposeof B. Here, saandijastorethematrix A;sbandijbstorethe matrixB. Thisroutinecomputesallcomponentsofthematrixproduct(whichmaybenon- sparse!), but stores only those whose magnitude exceeds thresh. On output, the arrays scandijc(whosemaximumsizeisinputas nmax)givetheproductmatrixinrow-index storagemode. Forsparsematrixmultiplication,thisroutinewilloftenbeprecededbyacall tosprstp,soastoconstructthetransposeofaknownmatrixinto sb,ijb. INTEGER i,ijma,ijmb,j,k,ma,mb,mbbREAL sumif (ija(1).ne.ijb(1)) pause ’sprstm sizes do not match’ k=ija(1) ijc(1)=kdo 14i=1,ija(1)-2 Loop overrows of A, do13j=1,ijb(1)-2 and rows of B. if(i.eq.j)then sum=sa(i)*sb(j) else sum=0.e0 endifmb=ijb(j) do 11ma=ija(i),ija(i+1)-1 Loopthroughelementsin A’srow. Convolutedlogic, following,accountsforthevariouscombinationsofdiagonalandoff-diagonalelements.ijma=ija(ma) if(ijma.eq.j)then sum=sum+sa(ma)*sb(j) else 2 if(mb.lt.ijb(j+1))then ijmb=ijb(mb) if(ijmb.eq.i)then sum=sum+sa(i)*sb(mb)mb=mb+1goto 2 else if(ijmb.lt.ijma)then mb=mb+1goto 2 else if(ijmb.eq.ijma)then sum=sum+sa(ma)*sb(mb) mb=mb+1goto 2 endif endif endif enddo 11 do12mbb=mb,ijb(j+1)-1 Exhausttheremainderof B’srow. if(ijb(mbb).eq.i)then sum=sum+sa(i)*sb(mbb) endif enddo 12 if(i.eq.j)then Whereto puttheanswer... sc(i)=sum else if(abs(sum).gt.thresh)then if(k.gt.nmax)pause ’sprstm: nmax to small’sc(k)=sum ijc(k)=j k=k+1 endif enddo 13 ijc(i+1)=k enddo 14 returnEND 2.7SparseLinearSystems 77Sample 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).Conjugate GradientMethod for a Sparse System So-called conjugate gradient methods provide a quite general means for solving the N×Nlinear system A·x=b (2.7.29 ) The attractiveness of these methods for large sparse systems is that they reference Aonly through its multiplication of a vector, or the multiplication of its transpose and a vector. Aswe have seen, these operations can be very ef ficient for a properly stored sparse matrix. You, the“owner”of the matrix A, can be asked to provide subroutines that perform these sparse matrixmultiplicationsasef ficientlyaspossible. We,the “grandstrategists ”supplythegeneral routine, linbcgbelow,thatsolvesthesetoflinearequations,(2.7.29),usingyoursubroutines. Thesimplest, “ordinary”conjugategradientalgorithm [11-13]solves(2.7.29)onlyinthe casethatAissymmetricandpositivede finite. Itisbasedontheideaofminimizingthefunction f(x)=1 2x·A·x−b·x (2.7.30 ) This function is minimized when its gradient ∇f=A·x−b (2.7.31 ) is zero, which is equivalent to (2.7.29). The minimization is carried out by generating a succession of search directions pkand improved minimizers xk. At each stage a quantity αk is found that minimizes f(xk+αkpk), andxk+1is set equal to the new point xk+αkpk. Thepkandxkare built up in such a way that xk+1is also the minimizer of fover the whole vector space of directions already taken, {p1,p2,...,pk}. After Niterations you arrive at the minimizer over the entire vector space, i.e., the solution to (2.7.29). Later, in §10.6, we will generalize this “ordinary”conjugate gradient algorithm to the minimization of arbitrary nonlinear functions. Here, where our interest is in solving linear,but not necessarily positive de finite or symmetric, equations, a different generalization is important, the biconjugate gradient method . This method does not, in general, have a simple connection with function minimization. It constructs four sequences of vectors, r k,rk,pk, pk,k=1,2,.... You supply the initial vectors r1andr1, and setp1=r1,p1=r1. Then you carry out the following recurrence: αk=rk·rk pk·A·pk rk+1=rk−αkA·pk rk+1=rk−αkAT·pk βk=rk+1·rk+1 rk·rk pk+1=rk+1+βkpk pk+1=rk+1+βkpk(2.7.32 ) This sequence of vectors satis fies thebiorthogonality condition ri·rj=ri·rj=0,j < i (2.7.33 ) and thebiconjugacy condition pi·A·pj=pi·AT·pj=0,j < i (2.7.34 ) There is also a mutual orthogonality, ri·pj=ri·pj=0,j < i (2.7.35 ) The proof of these properties proceeds by straightforward induction [14]. As long as the recurrence does not break down earlier because one of the denominators is zero, it must 78 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).terminate after m≤Nstepswith rm+1=rm+1=0. Thisisbasically because afteratmost Nsteps you run out of new orthogonal directions to the vectors you ’ve already constructed. To use the algorithm to solve the system (2.7.29), make an initial guess x1for the solution. Choose r1to be the residual r1=b−A·x1 (2.7.36 ) and choose r1=r1. Then form the sequence of improved estimates xk+1=xk+αkpk (2.7.37 ) while carrying out the recurrence (2.7.32). Equation (2.7.37) guarantees that rk+1from the recurrence is in fact the residual b−A·xk+1corresponding to xk+1. Sincerm+1=0, xm+1is the solution to equation (2.7.29). While there is no guarantee that this whole procedure will not break down or become unstable for general A, in practice this is rare. More importantly, the exact termination in at most Niterations occurs only with exact arithmetic. Roundoff error means that you should regard the process as a genuinely iterative procedure, to be halted when some appropriateerror criterion is met. Theordinaryconjugate gradientalgorithmisthespecialcaseofthebiconjugate gradient algorithm when Ais symmetric, and we choose r1=r1. Thenrk=rkandpk=pkfor all k;youcanomitcomputing themandhalvetheworkofthealgorithm. Thisconjugategradient version has the interpretation of minimizing equation (2.7.30). If Ais positive de finite as well as symmetric, the algorithm cannot break down (in theory!). The routine linbcgbelow indeed reduces to the ordinary conjugate gradient method if you input a symmetric A,b u t it does all the redundant computations. Another variant of the general algorithm corresponds to a symmetric but non-positive definiteA, with the choice r1=A·r1instead of r1=r1. In this case rk=A·rkand pk=A·pkfor all k. This algorithm is thus equivalent to the ordinary conjugate gradient algorithm,butwithalldotproducts a·breplacedby a·A·b. Itiscalledthe minimumresidual algorithm, because it corresponds to successive minimizations of the function Φ(x)=1 2r·r=1 2|A·x−b|2(2.7.38 ) wherethesuccessiveiterates xkminimize Φoverthesamesetofsearchdirections pkgenerated in the conjugate gradient method. This algorithm has been generalized in various ways forunsymmetric matrices. The generalized minimum residual method (GMRES; see [9,15])i s probably the most robust of these methods. Note that equation (2.7.38) gives ∇Φ(x)=AT·(A·x−b)( 2.7.39 ) Forany nonsingular matrix A,AT·Aissymmetric andpositivede finite. Youmighttherefore be tempted to solve equation (2.7.29) by applying the ordinary conjugate gradient algorithmto the problem (AT·A)·x=AT·b (2.7.40 ) Don’t! The condition number of the matrix AT·Ais the square of the condition number of A(see§2.6 for de finition of condition number). A large condition number both increases the number of iterations required, and limits the accuracy to which a solution can be obtained. Itis almost always better to apply the biconjugate gradient method to the original matrix A. So far we have said nothing about the rateof convergence of these methods. The ordinary conjugate gradient method works well for matrices that are well-conditioned, i.e.,“close”to the identity matrix. This suggests applying these methods to the preconditioned form of equation (2.7.29), (/tildewideA −1·A)·x=/tildewideA−1·b (2.7.41 ) Theidea isthat you might already be able tosolve your linear system easilyfor some /tildewideAclose toA, in which case /tildewideA−1·A≈1, allowing the algorithm to converge in fewer steps. The 2.7SparseLinearSystems 79Sample 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).matrix /tildewideAis called a preconditioner [11], and the overall scheme given here is known as the preconditioned biconjugate gradient method orPBCG. Forefficientimplementation,thePBCGalgorithmintroducesanadditionalsetofvectors zkandzkdefined by /tildewideA·zk=rkand /tildewideAT·zk=rk (2.7.42 ) and modi fies the definitions of αk,βk,pk, andpkin equation (2.7.32): αk=rk·zk pk·A·pk βk=rk+1·zk+1 rk·zk pk+1=zk+1+βkpk pk+1=zk+1+βkpk(2.7.43 ) Forlinbcg, below, we willask you to supply routines that solve the auxiliary linear systems (2.7.42). If you have no idea what to use for the preconditioner /tildewideA, then use the diagonal part ofA, or even the identity matrix, in which case the burden of convergence will be entirely on the biconjugate gradient method itself. Theroutine linbcg,below,isbasedonaprogramoriginallywrittenbyAnneGreenbaum. (See[13]for a different, less sophisticated, implementation.) There are a few wrinkles you should know about. What constitutes “good”convergence is rather application dependent. The routine linbcgtherefore provides for four possibilities, selected by setting the flagitolon input. Ifitol=1, iteration stops when the quantity |A·x−b|/|b|is less than the input quantity tol.I f itol=2, the required criterion is |/tildewideA−1·(A·x−b)|/|/tildewideA−1·b|<tol (2.7.44 ) Ifitol=3, the routine uses its own estimate of the error in x, and requires its magnitude, dividedbythemagnitudeof x,tobelessthan tol. Thesetting itol=4isthesameas itol=3, except that the largest (in absolute value) component of the error and largest component of x are used instead of the vector magnitude (that is, the L∞norm instead of the L2norm). You may need to experiment to find which of these convergence criteriais best for your problem. On output, erris the tolerance actually achieved. If the returned count iterdoes not indicate that the maximum number of allowed iterations itmaxwas exceeded, then err should be less than tol. If you want to do further iterations, leave all returned quantities as they are and call the routine again. The routine loses its memory of the spanned conjugategradient subspace between calls, however, so you should not force it to return more oftenthan about every Niterations. Finally, note that linbcgis furnished in double precision, since it will be usually be used when Nis quite large. SUBROUTINE linbcg(n,b,x,itol,tol,itmax,iter,err) INTEGER iter,itmax,itol,n,NMAX DOUBLE PRECISION err,tol,b(*),x(*),EPS Double precision isagood ideainthis rou- tine. PARAMETER (NMAX=1024,EPS=1.d-14) C USES atimes,asolve,snrm SolvesA·x=bforx(1:n),given b(1:n),bytheiterativebiconjugategradientmethod. On input x(1:n)should beset to aninitial guess of the solution (or allzeros); itolis 1,2,3,or4,specifyingwhichconvergencetestisapplied(seetext); itmaxisthemaximum number of allowed iterations; and tolis the desired convergence tolerance. On output, x(1:n)isresettotheimprovedsolution, iteristhenumberofiterationsactuallytaken, anderristheestimatederror. Thematrix Aisreferencedonlythroughtheuser-supplied routines atimes,whichcomputestheproductofeither Aoritstransposeonavector;and asolve,whichsolves /tildewideA·x=bor/tildewideAT·x=bforsomepreconditioner matrix /tildewideA(possibly the trivial diagonal part of A). INTEGER j 80 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).DOUBLE PRECISION ak,akden,bk,bkden,bknum,bnrm,dxnrm, * xnrm,zm1nrm,znrm,p(NMAX),pp(NMAX),r(NMAX),rr(NMAX), * z(NMAX),zz(NMAX),snrm iter=0 Calculateinitialresidual. call atimes(n,x,r,0) Inputto atimesisx(1:n),outputis r(1:n); thefinal 0indicatesthatthematrix(not itstranspose)istobeused.do11j=1,n r(j)=b(j)-r(j)rr(j)=r(j) enddo 11 C call atimes(n,r,rr,0) Uncomment this line to get the “minimum residual”variantofthealgorithm. if(itol.eq.1) then bnrm=snrm(n,b,itol)call asolve(n,r,z,0) Inputto asolveisr(1:n),outputis z(1:n); the final 0indicates that the matrix /tildewideA (notitstranspose)istobeused.else if (itol.eq.2) then call asolve(n,b,z,0)bnrm=snrm(n,z,itol) call asolve(n,r,z,0) else if (itol.eq.3.or.itol.eq.4) then call asolve(n,b,z,0)bnrm=snrm(n,z,itol) call asolve(n,r,z,0) znrm=snrm(n,z,itol) else pause ’illegal itol in linbcg’ endif 100 if (iter.le.itmax) then Mainloop. iter=iter+1 call asolve(n,rr,zz,1) Final 1indicatesuseoftransposematrix /tildewideA T. bknum=0.d0 do12j=1,n Calculatecoefficient bkanddirectionvectors pand pp. bknum=bknum+z(j)*rr(j) enddo 12 if(iter.eq.1) then do13j=1,n p(j)=z(j)pp(j)=zz(j) enddo 13 else bk=bknum/bkdendo 14j=1,n p(j)=bk*p(j)+z(j) pp(j)=bk*pp(j)+zz(j) enddo 14 endifbkden=bknum Calculate coefficient ak, new iterate x,a n d newresiduals randrr. call atimes(n,p,z,0) akden=0.d0 do 15j=1,n akden=akden+z(j)*pp(j) enddo 15 ak=bknum/akdencall atimes(n,pp,zz,1)do 16j=1,n x(j)=x(j)+ak*p(j) r(j)=r(j)-ak*z(j) rr(j)=rr(j)-ak*zz(j) enddo 16 call asolve(n,r,z,0) Solve /tildewideA·z=randcheckstoppingcriterion. if(itol.eq.1)then err=snrm(n,r,itol)/bnrm else if(itol.eq.2)then err=snrm(n,z,itol)/bnrm else if(itol.eq.3.or.itol.eq.4)then zm1nrm=znrm 2.7SparseLinearSystems 81Sample 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).znrm=snrm(n,z,itol) if(abs(zm1nrm-znrm).gt.EPS*znrm) then dxnrm=abs(ak)*snrm(n,p,itol) err=znrm/abs(zm1nrm-znrm)*dxnrm else err=znrm/bnrm Errormaynotbeaccurate,soloopagain. goto 100 endifxnrm=snrm(n,x,itol) if(err.le.0.5d0*xnrm) then err=err/xnrm else err=znrm/bnrm Errormaynotbeaccurate,soloopagain. goto 100 endif endif write (*,*) ’ iter=’,iter,’ err=’,err if(err.gt.tol) goto 100endifreturn END The routine linbcguses this short utility for computing vector norms: FUNCTION snrm(n,sx,itol) INTEGER n,itol,i,isamax DOUBLE PRECISION sx(n),snrm Computeoneoftwonormsforavector sx(1:n),assignaledby itol.U s e db y linbcg. if (itol.le.3)then snrm=0.do 11i=1,n Vector magnitude norm. snrm=snrm+sx(i)**2 enddo 11 snrm=sqrt(snrm) else isamax=1 do12i=1,n Largest component norm. if(abs(sx(i)).gt.abs(sx(isamax))) isamax=i enddo 12 snrm=abs(sx(isamax)) endifreturnEND So that the speci fications for the routines atimesandasolveare clear, we list here simple versions that assume a matrix Astored somewhere in row-index sparse format. SUBROUTINE atimes(n,x,r,itrnsp) INTEGER n,itrnsp,ija,NMAX DOUBLE PRECISION x(n),r(n),sa PARAMETER (NMAX=1000)COMMON /mat/ sa(NMAX),ija(NMAX) Thematrixisstoredsomewhere. C USES dsprsax,dsprstx DOUBLE PRECISION versionsof sprsaxandsprstx. if (itrnsp.eq.0) then call dsprsax(sa,ija,x,r,n) else call dsprstx(sa,ija,x,r,n) endifreturnEND 82 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).SUBROUTINE asolve(n,b,x,itrnsp) INTEGER n,itrnsp,ija,NMAX,i DOUBLE PRECISION x(n),b(n),sa PARAMETER (NMAX=1000)COMMON /mat/ sa(NMAX),ija(NMAX) Thematrixisstoredsomewhere. do 11i=1,n x(i)=b(i)/sa(i) The matrix /tildewideAis the diagonal part of A,s t o r e di n the first nelements of sa.S i n c et h et r a n s p o s e matrixhasthesamediagonal,theflag itrnspis not used.enddo 11 return END CITED REFERENCES AND FURTHER READING: Tewarson, R.P. 1973, Sparse Matrices (New York: Academic Press). [1] Jacobs,D.A.H.(ed.)1977, TheStateoftheArtinNumericalAnalysis (London:AcademicPress), Chapter I.3 (by J.K. Reid). [2] George,A.,andLiu,J.W.H.1981, ComputerSolutionofLargeSparsePositiveDefiniteSystems (Englewood Cliffs, NJ: Prentice-Hall). [3] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.). [4] IMSL Math/Library UsersManual (IMSL Inc., 2500CityWest Boulevard, HoustonTX77042).[5] Eisenstat,S.C.,Gursky,M.C.,Schultz,M.H.,andSherman,A.H.1977, YaleSparseMatrixPack- age,TechnicalReports112and114(YaleUniversityDepartmentofComputerScience).[6] Knuth,D.E.1968, FundamentalAlgorithms ,vol.1ofTheArtofComputerProgramming (Reading, MA: Addison-Wesley), §2.2.6. [7] Kincaid,D.R., Respess, J.R., Young, D.M., andGrimes, R.G. 1982, ACMTransactions onMath- ematical Software , vol. 8, pp. 302–322. [8] PCGPAK User’s Guide (New Haven: Scientific Computing Associates, Inc.). [9] Bentley, J. 1986, Programming Pearls (Reading, MA: Addison-Wesley), §9. [10] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), Chapters 4 and 10, particularly §§10.2–10.3. [11] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 8. [12] Baker, L. 1991, More C Tools for Scientists and Engineers (New York: McGraw-Hill). [13] Fletcher,R.1976,in NumericalAnalysisDundee1975 ,Lecture NotesinMathematics, vol.506, A. Dold and B Eckmann, eds. (Berlin: Springer-Verlag), pp. 73–89. [14] Saad, Y., and Schulz, M. 1986, SIAM Journal on Scientific and Statistical Computing , vol. 7, pp. 856–869. [15] Bunch, J.R., and Rose, D.J. (eds.) 1976, Sparse Matrix Computations (New York: Academic Press). Duff, I.S., and Stewart, G.W. (eds.) 1979, Sparse Matrix Proceedings 1978 (Philadelphia: S.I.A.M.). 2.8 Vandermonde Matrices and Toeplitz Matrices In§2.4 the case of a tridiagonal matrix was treated specially, because that particular type of linear system admits a solution in only of order Noperations, rather than of order N3for the general linear problem. When such particular types