Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 28
Текст из файла (страница 28)
Then we get(4.25)andHence, if we let L = (lij) be the unit lower-triangular matrix given in termsof the final content of the work array W by(4.26)then we can write these equations (for i = 1, . . . , n) in matrix form simplyas162MATRICES AND SYSTEMS OF LINEAR EQUATIONSandThis demonstrates the triangular factorization, in case no interchangesoccurred. If, on the other hand, interchanges did occur, then the finalcontent of W would have been unchanged had we carried out theseinterchanges at the outset and then applied Algorithm 4.2 without anyinterchanges.
This is so because all operations in the algorithm involve thesubtraction of a certain multiple of one row from certain other rows inorder to produce a zero in those other rows, and, for this, it does notmatter in which order we have written down the rows. The only thing thatmatters is that, once a row has been used as pivotal row, it is not modifiedany further, and, for this, we must keep apart from the others those rowsnot yet used as pivotal rows.Consequently, if interchanges do occur during execution of Algorithm4.2, then the matrices L and U obtained by the algorithm satisfyfor some appropriate permutation matrix P, i.e., then(4.27)and also(4.28)In terms of the vector p used in Algorithm 4.2 to record the interchanges made, the pkth equation is used as pivot equation during the kthstep. Hence P-1 should carry row p k to row k, all k.
This means thatPik = iPk, all k [see Theorem 4.5(iii)], if one really wanted to know. All thatmatters to us, though, is that, in terms of the output p from Algorithm 4.2,As a first application of the factorization point of view, we now look atthe possibility of splitting the process of solving Ax = b into two phases,the factorization phase in which the triangular factors L and U (and apossibly different order p of the rows) are derived, and the solving phaseduring which one first solves the triangular system(4.29)for y and then solves the triangular systemfor x, by back-substitution. Note that the right-hand side b enters only thesecond phase.
Hence, if the system is also to be solved for some otherright-hand sides, only the second phase needs to be repeated.According to (4.24), one solves (4.29) in Algorithm 4.2 by the steps4.4THE TRIANGULAR FACTORIZATION163In effect, this is like the back-substitution Algorithm 4.1 for solvingUx = y for x, except that the equations are gone through from first to last,since L is lower-triangular.We record the entire solving phase in the following:Algorithm 4.4: Forward- and back-substitution Given the final contentsof the first n columns of the working array W and the n-vector p ofAlgorithm 4.2 (applied to the system Ax = b); also, given the rightside b.The vector x = (xi ) now contains the solution of Ax = b.Note that, once again, both sums are sometimes empty.The practical significance of the preceding discussion becomes clearwhen we count (floating-point) operations in Algorithms 4.2 and 4.4. ByExercise 4.2-2, it takes n divisions, n(n - 1)/2 multiplications, and n(n 1)/2 additions to carry out the second loop in Algorithm 4.4.
The first looptakes the same number of operations, except that no divisions are required.Hence Algorithm 4.4 takesBy Exercise 4.2-4, this is exactly the number of operations required tomultiply an n × n matrix with an n-vector.By contrast,are necessary to calculate the first n columns of the final contents of theworking matrix W by Algorithm 4.2 (see Exercise 4.2-5). Hence the bulk ofthe work in solving Ax = b by elimination is needed to obtain the finalcontent of the working matrix W, namely,additions and the samenumber of multiplications/divisions, for large n.
The subsequent forwardand back-substitution takes an order of magnitude less operations, namely,additions and the same number of multiplications, per right side.Hence we can solve Ax = b for many different right sides (once we knowthe final content of W) in the time it takes to calculate the final contentof W.In this accounting of the work, we have followed tradition andcounted only floating-point operations. In particular, we have ignored164MATRICES AND SYSTEMS OF LINEAR EQUATIONSindex calculations, the cost of managing DO loops and other bookkeepingcosts, since these latter calculations used to be much faster than floatingpoint operations. This is not the case anymore on today’s computers, andthis way of accounting the work done may give an inaccurate picture (seeExercise 4.2-7). On the other hand, just how the work (as measured bycomputing time required) depends on the bookkeeping aspect of a program varies strongly from computer to computer and is therefore hard todiscuss in the generality of this textbook.A FORTRAN subroutine, called SUBST, which incorporates thesubstitution Algorithm 4.4, follows.SUBROUTINE SUBST ( W, IPIVOT, B, N, X )INTEGER IPIVOT(N)I,IP,JREAL B(N) ,W(N,N) ,X(N), SUMc****** I N P U T ******C W, IPIVOT, N ARE AS ON OUTPUT FROM F A C T 0 R , APPLIED TO THEMATRIX A OF ORDER N .CC B IS AN N-VECTOR, GIVING THE RIGHT SIDE OF THE SYSTEM TO BE SOLVED.C****** O U T P U T ******C X IS THE N-VECTOR SATISFYING A*X = B .C****** M E T H O D ******C ALGORITHM 4.4 IS USED, I.E., THE FACTORIZATION OF A CONTAINED INC W AND IPIVOT (AS GENERATED IN FACTOR ) IS USED TO SOLVE A*X = BC FOR X BY SOLVING TWO TRIANGULAR SYSTEMS.CIF [N .LE.
1) THENX(1) = B(1)/W(1,1)RETURNEND IFIP = IPIVOT(1)X(1) = B(IP)DO 15 I=2,NSUM = 0.DO 14 J=I,I-1SUM = W(I,J)*X(J) + SUM14IP = IPIVOT(I)X(I) = B(IP) - SUM15CX(N) = X(N)/W(N,N)DO 20 I=N-1,1,-lSUM = 0.DO 19 J=I+1,NSUM = W(I,J)*X(J) + SUM19X(I) = (X(I) - SUM)/W(I,I)20RETURNENDNext, we give a FORTRAN subroutine called FACTOR, which usesthe elimination Algorithm 4.2, with the pivoting strategy dictated by scaledpartial pivoting, to calculate a triangular factorization (if possible) for agiven N × N matrix A, storing the factorization in an N × N matrix W,and storing the pivoting strategy in an N-vector IPIVOT, ready for use inthe subroutine SUBST given earlier.
The user must provide an additionalN-vector D as a working space needed to store the “size” of the rows of A.If there is no further need for the matrix A and storage is scarce, then Aitself can be used for W in the argument list of the CALL statement (this isillegal in some FORTRAN dialects). The factorization will then replacethe original matrix in the array A.4.4THE TRIANGULAR FACTORIZATION165SUBROUTINE FACTOR ( W, N, D, IPIVOT, IFLAG )I,ISTAR,J,KINTEGER IFLAG,IPIVOT(N),REAL D(N) ,W(N,N), AWIKOD,COLMAX,RATIO,ROWMAX,TEMPC****** I N P U T ******C W ARRAY OF SIZE (N,N) CONTAINING THE MATRIX A OF ORDER N TO BEFACTORED.CC N THE ORDER OF THE MATRIXC****** W O R KA R E A ******C D A REAL VECTOR OF LENGTH N, TO HOLD ROW SIZESC****** O U T P U T ******C W ARRAY OF SIZE (N,N) CONTAINING THE LU FACTORIZATION OF P*A FORCSOME PERMUTATION MATRIX P SPECIFIED BY IPIVOT .C IPIVOT INTEGER VECTOR OF LENGTH N INDICATING THAT ROW IPIVOT(K)CWAS USED TO ELIMINATE X(K) , K=l,...,N .C IFLAGAN INTEGER,C= 1, IF AN EVEN NUMBER OF INTERCHANGES WAS CARRIED OUT,= -1, IF AN ODD NUMBER OF INTERCHANGES WAS CARRIED OUT,CC= 0, IF THE UPPER TRIANGULAR FACTOR HAS ONE OR MORE ZERO DIACGONAL ENTRIES.CTHUS, DETERMINANT(A) = IFLAG*W(1,1)*...*W(N,N) .IF IFLAG .NE.
0, THEN THE LINEAR SYSTEM A*X = B CAN BE SOLVED FORCCX BY ACALL SUBST (W, IPIVOT, B, N, X )CC****** M E T H O D ******C THE PROGRAM FOLLOWS ALGORITHM 4.2, USING SCALED PARTIAL PIVOTING.CIFLAG = 1INITIALIZE IPIVOT, DCDO 9 I=1,NIPIVOT(I) = IROWMAX = 0.DO 5 J=1,NROWMAX = AMAX1(ROWMAX,ABS(W(I,J)))5IF (ROWMAX .EQ. 0.) THENIFLAG = 0ROWMAX = 1.END IFD(I) = ROWMAX9IF (N .LE. 1)RETURNCFACTORIZATIONDO 20 K=1,N-1DETERMINE PIVOT ROW, THE ROW ISTAR .CCOLMAX = ABS(W(K,K) )/D(K)ISTAR = KDO 13 I=K+1,NAWIKOD = ABS(W(I,K) )/D(I)IF (AWIKOD .GT.
COLMAX) THENCOLMAX = AWIKODISTAR = IEND IF13 CONTINUEIF (COLMAX .EQ. 0.) THENIFLAG = 0ELSEIF (ISTAR .GT. K) THENMAKE K THE PIVOT ROW BY INTERCHANGING IT WITHCTHE CHOSEN ROW ISTAR .CIFLAG = -IFLAGI = IPIVOT(ISTAR)IPIVOT(ISTAR) = IPIVOT(K)IPIVOT(K) = ITEMP = D(ISTAR)D(ISTAR) = D(K)D(K) = TEMPDO 15 J=l,NTEMP = W(ISTAR,J)W(ISTAR,J) = W(K,J)15W(K,J) = TEMPEND IFCELIMINATE X(K) FROM ROWS K+1,...,N .16DO 19 I=K+1,NW(I,K) = W(I,K)/W(K,K)166MATRICES AND SYSTEMS OF LINEAR EQUATIONSRATIO = W(I,K)DO 19 J=K+1,NW(I,J) = W(I,J) - RATIO*W(K,J)CONTINUE19END IF20 CONTINUEIFLAG = 0IF (W(N,N) .EQ.
0.)RETURNENDThe preceding discussion points toward an efficient way to calculatethe inverse for a given invertible matrix A of order n. As was pointed out inSec. 4.1, for j = 1, . . . , n, the jth column A-1 ij of the inverse matrix A-1is the solution of the linear systemHence, to calculate A-1, one calls on FACTOR once, then solves each ofthe n systems Ax = ij, j = 1, .
. . , n, by Algorithm 4.4, that is, usingSUBST. Therefore, once the elimination is carried out, it takes only n · n2multiplications, and about the same number of additions, to find A-1 .Having given this simple prescription for calculating the inverse of amatrix, we hasten to point out that there is usually no good reason for evercalculating the inverse. It does at times happen in certain problems that theentries of A-1 have some special physical significance.
In the statisticaltreatment of the fitting of a function to observed data by the method ofleast squares, for example, the entries of a certain A-1 give informationabout the kinds and magnitudes of errors in the data. But whenever A-1 isneeded merely to calculate a vector A-1b ( a s i n s o l v i n g Ax = b) or amatrix product A-1B , A - 1 s h o u l d never be calculated explicitly. Rather,the substitution Algorithm 4.4 should be used to form these products.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
