46 Chapter 2.Solution of Linear Algebraic Equations method.Only partial pivoting(interchange of rows)can be implemented efficiently. However this is enough to make the method stable.This means,incidentally,that we don't actually decompose the matrix A into LU form,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 (its final application)is exactly the same as equation(2.3.13)except for the division in the latter equation;in both cases the upper limit of the sum isk =j-1(=i-1).This means that we don't have to commit ourselves as to whether the diagonal element Bi is the one that happens 8 to fall on the diagonal in the first instance,or whether one of the (undivided)ai's below it in the column,i=j+1,...,N,is to be"promoted"to become the diagonal B.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 B (pivot element),then do all the divisions by that element en masse.This is Crout's method with partial pivoting.Our implementation has one additional wrinkle:It MERICAI initially finds the largest element in each row,and subsequently (when it is looking for the maximal pivot element)scales the comparison as if we had initially scaled all the equations to make their maximum coefficient equal to unity;this is the implicit pivoting mentioned in $2.1. Press. THE ART #include <math.h> #include "nrutil.h" 9 #define TINY 1.0e-20 A small number. Program void ludcmp(float *a,int n,int *indx,float *d) Given a matrix a[1..n][1..n],this routine replaces it by the LU decomposition of a rowwise permutation of itself.a and n are input.a is 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 to dir 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 lubksb to solve linear equations or invert a matrix. OF SCIENTIFIC COMPUTING(ISBN 1988-19920 int i,imax,,k; float big,dum,sum,tempi float *vv; vv stores the implicit scaling of each row vv=vector(1,n); *d=1.0; No row interchanges yet. Numerical Recipes 10-621 43106 for(i=1;1<=n;i++)[ Loop over rows to get the implicit scaling informa- big=0.0; tion. for (j=1;j<=n;j++) (outside 膜 if ((temp=fabs(a[i][j]))>big)big-temp; oftware. if (big ==0.0)nrerror("Singular matrix in routine ludcmp"); No nonzero largest element. v[1]=1.0/b1g; Save the scaling. for(j=1;j<=n;j++){ This is the loop over columns of Crout's method. for (i=1;i<j;i++) This is equation (2.3.12)except for i=j. sum=a[i][j]; for (k=1;k<i;k++)sum -a[i][k]*a[k][j]; a[i][j]=sum; b1g=0.0 Initialize for the search for largest pivot element. for (i=j;i<-n;i++){ This is i=j of equation(2.3.12)andi=j+1...N sum=a[i][j]; of equation (2.3.13). for (k=1;k<j;k++)46 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). method. Only partial pivoting (interchange of rows) can be implemented efficiently. However this is enough to make the method stable. This means, incidentally, that we don’t actually decompose the matrix A into LU form, 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 (its final application) is exactly the same as equation (2.3.13) except for the division in the latter equation; in both cases the upper limit of the sum is k = j − 1 (= i − 1). This means that we don’t have to commit ourselves as to whether the diagonal element βjj is the one that happens to fall on the diagonal in the first instance, or whether one of the (undivided) α ij ’s below it in the column, i = j + 1,...,N, is to be “promoted” to become the diagonal β. 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 method with partial pivoting. Our implementation has one additional wrinkle: It initially finds the largest element in each row, and subsequently (when it is looking for the maximal pivot element) scales the comparison as if we had initially scaled all the equations to make their maximum coefficient equal to unity; this is the implicit pivoting mentioned in §2.1. #include <math.h> #include "nrutil.h" #define TINY 1.0e-20 A small number. void ludcmp(float **a, int n, int *indx, float *d) Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is 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; 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 lubksb to solve linear equations or invert a matrix. { int i,imax,j,k; float big,dum,sum,temp; float *vv; vv stores the implicit scaling of each row. vv=vector(1,n); *d=1.0; No row interchanges yet. for (i=1;i<=n;i++) { Loop over rows to get the implicit scaling informa￾big=0.0; tion. for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); No nonzero largest element. vv[i]=1.0/big; Save the scaling. } for (j=1;j<=n;j++) { This is the loop over columns of Crout’s method. for (i=1;i<j;i++) { This is equation (2.3.12) except for i = j. sum=a[i][j]; for (k=1;k<i;k++) sum -= a[i][k]*a[k][j]; a[i][j]=sum; } big=0.0; Initialize for the search for largest pivot element. for (i=j;i<=n;i++) { This is i = j of equation (2.3.12) and i = j + 1 ...N sum=a[i][j]; of equation (2.3.13). for (k=1;k<j;k++)