54 Chapter 2.Solution of Linear Algebraic Equations Some pivoting is possible within the storage limitations of bandec,and the above routine does take advantage of the opportunity.In general,when TINY is returned as a diagonal element of 0,then the original matrix (perhaps as modified by roundoff error) is in fact singular.In this regard,bandec is somewhat more robust than tridag above, which can fail algorithmically even for nonsingular matrices;bandec is thus also useful (with m=m2=1)for some ill-behaved tridiagonal systems. Once the matrix A has been decomposed,any number ofright-hand sides can be solved in turn by repeated calls to banbks,the backsubstitution routine whose analog in 82.3 is lubksb. #define SWAP(a,b){dum=(a)(a)=(b);(b)=dum; void banbks(float *a,unsigned long n,int mi,int m2,float **al, 81 nsigned1 Long indx☐，f1oatb[]) Given the arrays a,al,and indx as returned from bandec,and given a right-hand side vector b[1..n],solves the band diagonal linear equations A.x=b.The solution vector x overwrites b[1..n].The other input arrays are not modified,and can be left in place for successive calls with different right-hand sides. unsigned long i,k,1; 3 int mm; float dum; RECIPES mm=m1+m2+1; 1m1; for (k=1;k<=n;k++){ Forward substitution,unscrambling the permuted rows isindx[k]: as we go. Press. if (i !k)SWAP(b[k],b[i]) 1f(1<n)1+; for(1=k+1;i<=1;i++)b[i]-=a1k][i-k]*b[k]; 1=1; for(1n;1>=1;i--)[ Backsubstitution. SCIENTIFIC dum=b[i]; for(k=2;k<=1;k++)dum-=a[i][k]*b[k+i-1]; b[i]=dum/a[i][1]; to dir 1f(1<mm)1++: 2 192 COMPUTING (ISBN The routines bandec and banbks are based on the Handbook routines bandet/and bansoll in [1]. Fuunrggoirioh Numerica 10.621 43106 CITED REFERENCES AND FURTHER READING: Recipes Keller,H.B.1968,Numerical Methods for Two-Point Boundary-Value Problems (Waltham,MA: Blaisdell),p.74. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall), North Software. Example 5.4.3,p.166. Ralston,A.,and Rabinowitz,P.1978,A First Course in Numerical Analysis,2nd ed.(New York: McGraw-Hill),89.11. Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.ll of Handbook for Automatic Com- putation (New York:Springer-Verlag),Chapter 1/6.[1] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),84.3.54 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Some pivoting is possible within the storage limitations of bandec, and the above routine does take advantage of the opportunity. In general, when TINY is returned as a diagonal element of U, then the original matrix (perhaps as modified by roundoff error) is in fact singular. In this regard, bandec is somewhat more robust than tridag above, which can fail algorithmically even for nonsingular matrices; bandec is thus also useful (with m1 = m2 = 1) for some ill-behaved tridiagonal systems. Once the matrix A has been decomposed, any number of right-hand sides can be solved in turn by repeated calls to banbks, the backsubstitution routine whose analog in §2.3 is lubksb. #define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;} void banbks(float **a, unsigned long n, int m1, int m2, float **al, unsigned long indx[], float b[]) Given the arrays a, al, and indx as returned from bandec, and given a right-hand side vector b[1..n], solves the band diagonal linear equations A · x = b. The solution vector x overwrites b[1..n]. The other input arrays are not modified, and can be left in place for successive calls with different right-hand sides. { unsigned long i,k,l; int mm; float dum; mm=m1+m2+1; l=m1; for (k=1;k<=n;k++) { Forward substitution, unscrambling the permuted rows i=indx[k]; as we go. if (i != k) SWAP(b[k],b[i]) if (l < n) l++; for (i=k+1;i<=l;i++) b[i] -= al[k][i-k]*b[k]; } l=1; for (i=n;i>=1;i--) { Backsubstitution. dum=b[i]; for (k=2;k<=l;k++) dum -= a[i][k]*b[k+i-1]; b[i]=dum/a[i][1]; if (l < mm) l++; } } The routines bandec and banbks are based on the Handbook routines bandet1 and bansol1 in [1]. CITED REFERENCES AND FURTHER READING: Keller, H.B. 1968, Numerical Methods for Two-Point Boundary-Value Problems (Waltham, MA: Blaisdell), p. 74. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Example 5.4.3, p. 166. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York: McGraw-Hill), §9.11. Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Com￾putation (New York: Springer-Verlag), Chapter I/6. [1] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §4.3