Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 15
Текст из файла (страница 15)
In otherwords, the matrix solution of[A] · [x1 x2 x3 Y] = [b1 b2 b3 1](2.1.2)where A and Y are square matrices, the bi ’s and xi ’s are column vectors, and 1 isthe identity matrix, simultaneously solves the linear setsA · x1 = b 1A · x2 = b2A · x3 = b3(2.1.3)andA·Y = 1(2.1.4)Now it is also elementary to verify the following facts about (2.1.1):• Interchanging any two rows of A and the corresponding rows of the b’sand of 1, does not change (or scramble in any way) the solution x’s andY. Rather, it just corresponds to writing the same set of linear equationsin a different order.• Likewise, the solution set is unchanged and in no way scrambled if wereplace any row in A by a linear combination of itself and any other row,as long as we do the same linear combination of the rows of the b’s and 1(which then is no longer the identity matrix, of course).• Interchanging any two columns of A gives the same solution set onlyif we simultaneously interchange corresponding rows of the x’s and ofY.
In other words, this interchange scrambles the order of the rows inthe solution. If we do this, we will need to unscramble the solution byrestoring the rows to their original order.Gauss-Jordan elimination uses one or more of the above operations to reducethe matrix A to the identity matrix. When this is accomplished, the right-hand sidebecomes the solution set, as one sees instantly from (2.1.2).38Chapter 2.Solution of Linear Algebraic EquationsPivotingIn “Gauss-Jordan elimination with no pivoting,” only the second operation inthe above list is used. The first row is divided by the element a11 (this being atrivial linear combination of the first row with any other row — zero coefficient forthe other row).
Then the right amount of the first row is subtracted from each otherrow to make all the remaining ai1 ’s zero. The first column of A now agrees withthe identity matrix. We move to the second column and divide the second row bya22 , then subtract the right amount of the second row from rows 1, 3, and 4, so as tomake their entries in the second column zero. The second column is now reducedto the identity form. And so on for the third and fourth columns. As we do theseoperations to A, we of course also do the corresponding operations to the b’s and to1 (which by now no longer resembles the identity matrix in any way!).Obviously we will run into trouble if we ever encounter a zero element on the(then current) diagonal when we are going to divide by the diagonal element.
(Theelement that we divide by, incidentally, is called the pivot element or pivot.) Not soobvious, but true, is the fact that Gauss-Jordan elimination with no pivoting (no use ofthe first or third procedures in the above list) is numerically unstable in the presenceof any roundoff error, even when a zero pivot is not encountered. You must never doGauss-Jordan elimination (or Gaussian elimination, see below) without pivoting!So what is this magic pivoting? Nothing more than interchanging rows (partialpivoting) or rows and columns (full pivoting), so as to put a particularly desirableelement in the diagonal position from which the pivot is about to be selected.
Sincewe don’t want to mess up the part of the identity matrix that we have already built up,we can choose among elements that are both (i) on rows below (or on) the one thatis about to be normalized, and also (ii) on columns to the right (or on) the columnwe are about to eliminate. Partial pivoting is easier than full pivoting, because wedon’t have to keep track of the permutation of the solution vector. Partial pivotingmakes available as pivots only the elements already in the correct column.
It turnsout that partial pivoting is “almost” as good as full pivoting, in a sense that can bemade mathematically precise, but which need not concern us here (for discussionand references, see [1]). To show you both variants, we do full pivoting in the routinein this section, partial pivoting in §2.3.We have to state how to recognize a particularly desirable pivot when we seeone. The answer to this is not completely known theoretically. It is known, boththeoretically and in practice, that simply picking the largest (in magnitude) availableelement as the pivot is a very good choice. A curiosity of this procedure, however, isthat the choice of pivot will depend on the original scaling of the equations. If we takethe third linear equation in our original set and multiply it by a factor of a million, itis almost guaranteed that it will contribute the first pivot; yet the underlying solutionof the equations is not changed by this multiplication! One therefore sometimes seesroutines which choose as pivot that element which would have been largest if theoriginal equations had all been scaled to have their largest coefficient normalized tounity.
This is called implicit pivoting. There is some extra bookkeeping to keep trackof the scale factors by which the rows would have been multiplied. (The routines in§2.3 include implicit pivoting, but the routine in this section does not.)Finally, let us consider the storage requirements of the method. With a littlereflection you will see that at every stage of the algorithm, either an element of A is2.1 Gauss-Jordan Elimination39predictably a one or zero (if it is already in a part of the matrix that has been reducedto identity form) or else the exactly corresponding element of the matrix that startedas 1 is predictably a one or zero (if its mate in A has not been reduced to the identityform). Therefore the matrix 1 does not have to exist as separate storage: The matrixinverse of A is gradually built up in A as the original A is destroyed.
Likewise,the solution vectors x can gradually replace the right-hand side vectors b and sharethe same storage, since after each column in A is reduced, the corresponding rowentry in the b’s is never again used.Here is the routine for Gauss-Jordan elimination with full pivoting:#include <math.h>#include "nrutil.h"#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}void gaussj(float **a, int n, float **b, int m)Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) above. a[1..n][1..n]is the input matrix.
b[1..n][1..m] is input containing the m right-hand side vectors. Onoutput, a is replaced by its matrix inverse, and b is replaced by the corresponding set of solutionvectors.{int *indxc,*indxr,*ipiv;int i,icol,irow,j,k,l,ll;float big,dum,pivinv,temp;indxc=ivector(1,n);The integer arrays ipiv, indxr, and indxc areindxr=ivector(1,n);used for bookkeeping on the pivoting.ipiv=ivector(1,n);for (j=1;j<=n;j++) ipiv[j]=0;for (i=1;i<=n;i++) {This is the main loop over the columns to bebig=0.0;reduced.for (j=1;j<=n;j++)This is the outer loop of the search for a pivotif (ipiv[j] != 1)element.for (k=1;k<=n;k++) {if (ipiv[k] == 0) {if (fabs(a[j][k]) >= big) {big=fabs(a[j][k]);irow=j;icol=k;}} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");}++(ipiv[icol]);We now have the pivot element, so we interchange rows, if needed, to put the pivotelement on the diagonal.
The columns are not physically interchanged, only relabeled:indxc[i], the column of the ith pivot element, is the ith column that is reduced, whileindxr[i] is the row in which that pivot element was originally located. If indxr[i] =indxc[i] there is an implied column interchange. With this form of bookkeeping, thesolution b’s will end up in the correct order, and the inverse matrix will be scrambledby columns.if (irow != icol) {for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])}indxr[i]=irow;We are now ready to divide the pivot row by theindxc[i]=icol;pivot element, located at irow and icol.if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");pivinv=1.0/a[icol][icol];a[icol][icol]=1.0;for (l=1;l<=n;l++) a[icol][l] *= pivinv;for (l=1;l<=m;l++) b[icol][l] *= pivinv;40Chapter 2.Solution of Linear Algebraic Equationsfor (ll=1;ll<=n;ll++)Next, we reduce the rows...if (ll != icol) {...except for the pivot one, of course.dum=a[ll][icol];a[ll][icol]=0.0;for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;}}This is the end of the main loop over columns of the reduction.
It only remains to unscramble the solution in view of the column interchanges. We do this by interchanging pairs ofcolumns in the reverse order that the permutation was built up.for (l=n;l>=1;l--) {if (indxr[l] != indxc[l])for (k=1;k<=n;k++)SWAP(a[k][indxr[l]],a[k][indxc[l]]);}And we are done.free_ivector(ipiv,1,n);free_ivector(indxr,1,n);free_ivector(indxc,1,n);}Row versus Column Elimination StrategiesThe above discussion can be amplified by a modest amount of formalism. Rowoperations on a matrix A correspond to pre- (that is, left-) multiplication by some simplematrix R. For example, the matrix R with components 1 if i = j and i = 2, 4(2.1.5)Rij = 1 if i = 2, j = 4 1 if i = 4, j = 20 otherwiseeffects the interchange of rows 2 and 4.
Gauss-Jordan elimination by row operations alone(including the possibility of partial pivoting) consists of a series of such left-multiplications,yielding successivelyA·x =b(· · · R3 · R2 · R1 · A) · x = · · · R3 · R2 · R1 · b(1) · x = · · · R3 · R2 · R1 · b(2.1.6)x = · · · R3 · R2 · R1 · bThe key point is that since the R’s build from right to left, the right-hand side is simplytransformed at each stage from one vector to another.Column operations, on the other hand, correspond to post-, or right-, multiplicationsby simple matrices, call them C. The matrix in equation (2.1.5), if right-multiplied onto amatrix A, will interchange A’s second and fourth columns. Elimination by column operationsinvolves (conceptually) inserting a column operator, and also its inverse, between the matrixA and the unknown vector x:A·x= bA · C1 · C−11 ·x = b−1A · C1 · C2 · C−12 · C1 · x = b(A · C1 · C2 ·C3 · · ·) · · · C−13(1) · · · C−13·C−12·C−12·C−11·x= b·C−11·x= b(2.1.7)2.2 Gaussian Elimination with Backsubstitution41which (peeling of the C−1 ’s one at a time) implies a solutionx = C 1 · C2 · C3 · · · b(2.1.8)Notice the essential difference between equation (2.1.8) and equation (2.1.6).
In thelatter case, the C’s must be applied to b in the reverse order from that in which they becomeknown. That is, they must all be stored along the way. This requirement greatly reducesthe usefulness of column operations, generally restricting them to simple permutations, forexample in support of full pivoting.CITED REFERENCES AND FURTHER READING:Wilkinson, J.H. 1965, The Algebraic Eigenvalue Problem (New York: Oxford University Press). [1]Carnahan, B., Luther, H.A., and Wilkes, J.O.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
