2.3 LU Decomposition and Its Applications 47 sum -a[i][k]*a[k][j]; a[i][i]=sum; if (dum=vv[i]*fabs(sum))>=big){ Is the figure of merit for the pivot better than the best so far? big=dum; imax=i; if (i !imax){ Do we need to interchange rows? for(k=1;k<=n;k++){ Yes,do so... dum=a [imax][k]; a[imax][k]=a[j][k]; http://www.nr. a[j][k]=dum; *d=-(*d) .and change the parity of d. 83g vv[imax]=vv[j]; Also interchange the scale factor. granted for 2 indx[i]=imax; 1.800 if (a[j][j]==0.0)a[j][j]=TINY; If the pivot element is zero the matrix is singular (at least to the precision of the algorithm).For some applications on singular matrices,it is desirable to substitute from NUMERICAL RECIPESI 19881992 9 TINY for zero if (j !n){ Now,finally,divide by the pivot element. dum=1.0/(a[j][j]); (Nort server for(1=j+1;1<=n;1+)a[i][j]*=dum; THE Go back for the next column in the reduction. America computer, free_vector(vv,1,n); ART 9 Programs Here is the routine for forward substitution and backsubstitution,implementing equations (2.3.6)and (2.3.7). void lubksb(float **a,int n,int *indx,float b[]) to dir Solves the set of n linear equations A.X=B.Here a[1..n][1..n]is input,not as the matrix A but rather as its LU decomposition,determined by the routine ludcmp.indx[1..n]is input as the permutation vector returned by ludcmp.b[1..n]is input as the right-hand side vector SCIENTIFIC COMPUTING(ISBN B,and returns with the solution vector X.a,n,and indx are not modified by this routine 1881892 and can be left in place for successive calls with different right-hand sides b.This routine takes into account the possibility that b will begin with many zero elements,so it is efficient for use in matrix inversion. 10-621 int i,ii=0,ip,j; Numerical Recipes -43108 float sum; for (i=1;i<=n;i++){ When ii is set to a positive value,it will become the ip=indx[i]; index of the first nonvanishing element of b.We now (outside sum=b[ip]; do the forward substitution,equation (2.3.6).The Software. b[ip]=b[i]; only new wrinkle is to unscramble the permutation 1f(11) as we go. for (j=ii;j<=i-1;j++)sum-=a[i][j]*b[j]; else if (sum)ii=i; A nonzero element was encountered,so from now on we b[i]-sum; will have to do the sums in the loop above. for(i=n;1>=1;1-)[ Now we do the backsubstitution,equation (2.3.7). sum=b[i]; for (j=i+1;j<=n;j++)sum -a[i][j]*b[j]; b[i]=sum/a[i][i]; Store a component of the solution vector X. All done!2.3 LU Decomposition and Its Applications 47 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). sum -= a[i][k]*a[k][j]; a[i][j]=sum; if ( (dum=vv[i]*fabs(sum)) >= big) { Is the figure of merit for the pivot better than the best so far? big=dum; imax=i; } } if (j != imax) { Do we need to interchange rows? for (k=1;k<=n;k++) { Yes, do so... dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; } *d = -(*d); ...and change the parity of d. vv[imax]=vv[j]; Also interchange the scale factor. } indx[j]=imax; if (a[j][j] == 0.0) a[j][j]=TINY; If the pivot element is zero the matrix is singular (at least to the precision of the algorithm). For some applications on singular matrices, it is desirable to substitute TINY for zero. if (j != n) { Now, finally, divide by the pivot element. dum=1.0/(a[j][j]); for (i=j+1;i<=n;i++) a[i][j] *= dum; } } Go back for the next column in the reduction. free_vector(vv,1,n); } Here is the routine for forward substitution and backsubstitution, implementing equations (2.3.6) and (2.3.7). void lubksb(float **a, int n, int *indx, float b[]) Solves the set of n linear equations A·X = B. Here a[1..n][1..n] is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx[1..n] is 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, and indx are 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 b will begin with many zero elements, so it is efficient for use in matrix inversion. { int i,ii=0,ip,j; float sum; for (i=1;i<=n;i++) { When ii is set to a positive value, it will become the index of the first nonvanishing element of b. We now do the forward substitution, equation (2.3.6). The only new wrinkle is to unscramble the permutation as we go. ip=indx[i]; sum=b[ip]; b[ip]=b[i]; if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; else if (sum) ii=i; A nonzero element was encountered, so from now on we b[i]=sum; will have to do the sums in the loop above. } for (i=n;i>=1;i--) { Now we do the backsubstitution, equation (2.3.7). sum=b[i]; for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; b[i]=sum/a[i][i]; Store a component of the solution vector X. } All done! }