48 Chapter 2.Solution of Linear Algebraic Equations The LU decomposition in ludcmp requires about N3 executions of the inner loops (each with one multiply and one add).This is thus the operation count for solving one (or a few)right-hand sides,and is a factor of 3 better than the Gauss-Jordan routine gaussj which was given in 82.1,and a factor of 1.5 better than a Gauss-Jordan routine (not given)that does not compute the inverse matrix For inverting a matrix,the total count (including the forward and backsubstitution as discussed following equation 2.3.7 above)is (++)N3N3,the same as gaussj To summarize,this is the preferred way to solve the linear set of equations A·x=b nted for float **a,*b,d; int n,*indx; ludcmp(a,n,indx,&d); lubksb(a,n,indx,b); The answer x will be given back in b.Your original matrix A will have been destroyed. If you subsequently want to solve a set of equations with the same A but a R 令 Press. different right-hand side b,you repeat only 9 lubksb(a,n,indx,b); SCIENTIFIC not,of course,with the original matrix A.but with a and indx as were already 6 set by ludcmp. Inverse of a Matrix 1920 COMPUTING (ISBN Using the above L0 decomposition and backsubstitution routines,it is com- Numerical Recipes 021 pletely straightforward to find the inverse of a matrix column by column. 43108 #define N… float **a,**y,d,*col; int i,j,*indx; (outside Software. ludcmp(a,N,indx,&d); Decompose the matrix just once. for(j=1;j<=N;j++){ Find inverse by columns. for(1=1;i<=N;i++) co1[1]=0.0: Amer ying of co1[j]=1.0: lubksb(a,N,indx,col); for(i=1;i<=N;i+)y[i][j]=co1[1]; The matrix y will now contain the inverse of the original matrix a,which will have been destroyed.Alternatively,there is nothing wrong with using a Gauss-Jordan routine like gaussj (82.1)to invert a matrix in place,again destroying the original Both methods have practically the same operations count.48 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). The LU decomposition in ludcmp requires about 1 3N3 executions of the inner loops (each with one multiply and one add). This is thus the operation count for solving one (or a few) right-hand sides, and is a factor of 3 better than the Gauss-Jordan routine gaussj which was given in §2.1, and a factor of 1.5 better than a Gauss-Jordan routine (not given) that does not compute the inverse matrix. For inverting a matrix, the total count (including the forward and backsubstitution as discussed following equation 2.3.7 above) is ( 1 3 + 1 6 + 1 2 )N3 = N3, the same as gaussj. To summarize, this is the preferred way to solve the linear set of equations A · x = b: float **a,*b,d; int n,*indx; ... ludcmp(a,n,indx,&d); lubksb(a,n,indx,b); The answer x will be given back in b. Your original matrix A will have been destroyed. If you subsequently want to solve a set of equations with the same A but a different right-hand side b, you repeat only lubksb(a,n,indx,b); not, of course, with the original matrix A, but with a and indx as were already set by ludcmp. Inverse of a Matrix Using the above LU decomposition and backsubstitution routines, it is com￾pletely straightforward to find the inverse of a matrix column by column. #define N ... float **a,**y,d,*col; int i,j,*indx; ... ludcmp(a,N,indx,&d); Decompose the matrix just once. for(j=1;j<=N;j++) { Find inverse by columns. for(i=1;i<=N;i++) col[i]=0.0; col[j]=1.0; lubksb(a,N,indx,col); for(i=1;i<=N;i++) y[i][j]=col[i]; } The matrix y will now contain the inverse of the original matrix a, which will have been destroyed. Alternatively, there is nothing wrong with using a Gauss-Jordan routine like gaussj (§2.1) to invert a matrix in place, again destroying the original. Both methods have practically the same operations count