2.3 LU Decomposition and Its Applications 49 Incidentally,if you ever have the need to compute A-1.B from matrices A and B,you should LU decompose A and then backsubstitute with the columns of B instead of with the unit vectors that would give A's inverse.This saves a whole matrix multiplication,and is also more accurate. Determinant of a Matrix The determinant of an Lo decomposed matrix is just the product of the diagonal elements, det (2.3.15) 1=1 We don't,recall,compute the decomposition of the original matrix,but rather a decomposition of a rowwise permutation of it.Luckily,we have kept track of whether the number of row interchanges was even or odd,so we just preface the RECIPES I 2 product by the corresponding sign.(You now finally know the purpose of setting d in the routine ludcmp.) Calculation of a determinant thus requires one call to ludcmp,with no subse- UnN电.U quent backsubstitutions by lubksb. #define N... Programs float米*a,d; int j,*indx; SCIENTIFIC ludcmp(a,N,indx,&d); This returns d as士l. for(j=1;j<=;j+)d*a[j][j]; The variable d now contains the determinant of the original matrix a,which will have been destroyed. COMPUTING (ISBN 19891892 For a matrix of any substantial size,it is quite likely that the determinant will overflow or underflow your computer's floating-point dynamic range.In this case you can modify the loop of the above fragment and (e.g.divide by powers of ten, to keep track of the scale separately,or(e.g.)accumulate the sum of logarithms of Fuurggoglrion Numerical Recipes 10-621 43106 the absolute values of the factors and the sign separately. (outside Complex Systems of Equations North If your matrix A is real,but the right-hand side vector is complex,say b+id,then (i) LU decompose A in the usual way,(i)backsubstitute b to get the real part of the solution vector,and(iii)backsubstitute d to get the imaginary part of the solution vector. If the matrix itself is complex,so that you want to solve the system (A+C)·(x+y)=(b+id) (2.3.16) then there are two possible ways to proceed.The best way is to rewrite ludcmp and lubksb as complex routines.Complex modulus substitutes for absolute value in the construction of the scaling vector vv and in the search for the largest pivot elements.Everything else goes through in the obvious way,with complex arithmetic used as needed.(See $81.2 and 5.4 for discussion of complex arithmetic in C.)2.3 LU Decomposition and Its Applications 49 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). Incidentally, if you ever have the need to compute A −1 · B from matrices A and B, you should LU decompose A and then backsubstitute with the columns of B instead of with the unit vectors that would give A’s inverse. This saves a whole matrix multiplication, and is also more accurate. Determinant of a Matrix The determinant of an LU decomposed matrix is just the product of the diagonal elements, det = N j=1 βjj (2.3.15) We don’t, recall, compute the decomposition of the original matrix, but rather a decomposition of a rowwise permutation of it. Luckily, we have kept track of whether the number of row interchanges was even or odd, so we just preface the product by the corresponding sign. (You now finally know the purpose of setting d in the routine ludcmp.) Calculation of a determinant thus requires one call to ludcmp, with no subse￾quent backsubstitutions by lubksb. #define N ... float **a,d; int j,*indx; ... ludcmp(a,N,indx,&d); This returns d as ±1. for(j=1;j<=N;j++) d *= a[j][j]; The variable d now contains the determinant of the original matrix a, which will have been destroyed. For a matrix of any substantial size, it is quite likely that the determinant will overflow or underflow your computer’s floating-point dynamic range. In this case you can modify the loop of the above fragment and (e.g.) divide by powers of ten, to keep track of the scale separately, or (e.g.) accumulate the sum of logarithms of the absolute values of the factors and the sign separately. Complex Systems of Equations If your matrix A is real, but the right-hand side vector is complex, say b + id, then (i) LU decompose A in the usual way, (ii) backsubstitute b to get the real part of the solution vector, and (iii) backsubstitute d to get the imaginary part of the solution vector. If the matrix itself is complex, so that you want to solve the system (A + iC) · (x + iy)=(b + id) (2.3.16) then there are two possible ways to proceed. The best way is to rewrite ludcmp and lubksb as complex routines. Complex modulus substitutes for absolute value in the construction of the scaling vector vv and in the search for the largest pivot elements. Everything else goes through in the obvious way, with complex arithmetic used as needed. (See §§1.2 and 5.4 for discussion of complex arithmetic in C.)