2.4 Tridiagonal and Band Diagonal Systems of Equations 53 It is not possible to store the LU decomposition of a band diagonal matrix A quite as compactly as the compact form of A itself.The decomposition (essentially by Crout's method,see 82.3)produces additional nonzero"fill-ins."One straightforward storage scheme is to return the upper triangular factor (in the same space that A previously occupied,and to return the lower triangular factor(L)in a separate compact matrix of size N x m.The diagonal elements of U (whose product,times d=+1,gives the determinant)are returned in the first column of A's storage space. The following routine,bandec,is the band-diagonal analog of ludcmp in 82.3: #include <math.h> #define SWAP(a,b){dum=(a);(a)=(b);(b)=dum;} #define TINY 1.0e-20 83 void bandec(float **a,unsigned long n,int m1,int m2,float **al unsigned long indx,float *d) 鱼 granted for 19881992 Given an n x n band diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal rows, 1800 豆 compactly stored in the array a[1..n][1..m1+m2+1]as described in the comment for routine banmul,this routine constructs an LU decomposition of a rowwise permutation of A.The upper triangular matrix replaces a,while the lower triangular matrix is returned in al[1..n][1..m1] indx[1..n]is an output vector which records the row permutation effected by the partial from NUMERICAL RECIPES IN pivoting:d is output as +1 depending on whether the number of row interchanges was even or odd,respectively.This routine is used in combination with banbks to solve band-diagonal sets of equations. unsigned long i,j,k,1; server computer, (North America to make one paper THE int mm; float dum; ART 是 mm=m1+m2+1; st st Programs 1=m1; for(i=1;1<=m1;i++)[ Rearrange the storage a bit for(j=m1+2-i;j<=mm;j++)a[i][j-1]=a[i][j]; 1--; email for(j=mm-1;j<=mm;j++)a[i][j]=0.0; Copyright (C) *d=1.0: 1=m1; for (k=1;k<=nik++){ For each row... to directcustsen OF SCIENTIFIC COMPUTING(ISBN dum=a[k][1]; 1=k; 1f(1<n)1++; for（j=k+1;j<=1;j+){ Find the pivot element if (fabs(a[i][1])>fabs(dum)){ dum=a[j][1]; @cambridge.org 1988-1992 by Numerical Recipes 10-:6211 43106 1=j; indx[k]=i; (outside if (dum =0.0)a[k][1]=TINY; Software. Matrix is algorithmically singular,but proceed anyway with TINY pivot (desirable in some applications). 1f(1I=k)[ Interchange rows Ame *d=-(*d); for (j=1;j<=mm;j++)SWAP(a[k][j],a[i][i]) for(1=k+1;1<=1;i+) Do the elimination dum=a[i][1]/a[k][1]; al[k][i-k]=dum; for (j=2;j<=mm;j++)a[i][j-1]=a[i][j]-dum*a[k][j]; a[i][mm]=0.0;2.4 Tridiagonal and Band Diagonal Systems of Equations 53 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). It is not possible to store the LU decomposition of a band diagonal matrix A quite as compactly as the compact form of A itself. The decomposition (essentially by Crout’s method, see §2.3) produces additional nonzero “fill-ins.” One straightforward storage scheme is to return the upper triangular factor (U) in the same space that A previously occupied, and to return the lower triangular factor (L) in a separate compact matrix of size N × m1. The diagonal elements of U (whose product, times d = ±1, gives the determinant) are returned in the first column of A’s storage space. The following routine, bandec, is the band-diagonal analog of ludcmp in §2.3: #include <math.h> #define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;} #define TINY 1.0e-20 void bandec(float **a, unsigned long n, int m1, int m2, float **al, unsigned long indx[], float *d) Given an n × n band diagonal matrix A with m1 subdiagonal rows and m2 superdiagonal rows, compactly stored in the array a[1..n][1..m1+m2+1] as described in the comment for routine banmul, this routine constructs an LU decomposition of a rowwise permutation of A. The upper triangular matrix replaces a, while the lower triangular matrix is returned in al[1..n][1..m1]. indx[1..n] is an output vector which records the row permutation effected by the partial pivoting; d is output as ±1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with banbks to solve band-diagonal sets of equations. { unsigned long i,j,k,l; int mm; float dum; mm=m1+m2+1; l=m1; for (i=1;i<=m1;i++) { Rearrange the storage a bit. for (j=m1+2-i;j<=mm;j++) a[i][j-l]=a[i][j]; l--; for (j=mm-l;j<=mm;j++) a[i][j]=0.0; } *d=1.0; l=m1; for (k=1;k<=n;k++) { For each row... dum=a[k][1]; i=k; if (l < n) l++; for (j=k+1;j<=l;j++) { Find the pivot element. if (fabs(a[j][1]) > fabs(dum)) { dum=a[j][1]; i=j; } } indx[k]=i; if (dum == 0.0) a[k][1]=TINY; Matrix is algorithmically singular, but proceed anyway with TINY pivot (desirable in some applications). if (i != k) { Interchange rows. *d = -(*d); for (j=1;j<=mm;j++) SWAP(a[k][j],a[i][j]) } for (i=k+1;i<=l;i++) { Do the elimination. dum=a[i][1]/a[k][1]; al[k][i-k]=dum; for (j=2;j<=mm;j++) a[i][j-1]=a[i][j]-dum*a[k][j]; a[i][mm]=0.0; } } }