c2-1 (779459), страница 2
Текст из файла (страница 2)
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] 6=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;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).predictably 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.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;}}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 6= 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·C−12·C−11·x= b−1−1(1) · · · C−13 · C2 · C1 · x = b(2.1.7)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).}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);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)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. 1969, Applied Numerical Methods (New York:Wiley), Example 5.2, p. 282.Bevington, P.R. 1969, Data Reduction and Error Analysis for the Physical Sciences (New York:McGraw-Hill), Program B-2, p. 298.Westlake, J.R. 1968, A Handbook of Numerical Matrix Inversion and Solution of Linear Equations(New York: Wiley).Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis, 2nd ed. (New York:McGraw-Hill), §9.3–1.2.2 Gaussian Elimination with BacksubstitutionThe usefulness of Gaussian elimination with backsubstitution is primarilypedagogical. It stands between full elimination schemes such as Gauss-Jordan, andtriangular decomposition schemes such as will be discussed in the next section.Gaussian elimination reduces a matrix not all the way to the identity matrix, butonly halfway, to a matrix whose components on the diagonal and above (say) remainnontrivial.
Let us now see what advantages accrue.Suppose that in doing Gauss-Jordan elimination, as described in §2.1, we ateach stage subtract away rows only below the then-current pivot element. When a22is the pivot element, for example, we divide the second row by its value (as before),but now use the pivot row to zero only a32 and a42 , not a12 (see equation 2.1.1).Suppose, also, that we do only partial pivoting, never interchanging columns, so thatthe order of the unknowns never needs to be modified.Then, when we have done this for all the pivots, we will be left with a reducedequation that looks like this (in the case of a single right-hand side vector): 0 0x1b1a11 a012 a013 a014 0 a022 a023 a024 x2 b02 (2.2.1)· = 0 00 a033 a034x3b3000 a044x4b04Here the primes signify that the a’s and b’s do not have their original numericalvalues, but have been modified by all the row operations in the elimination to thispoint.
The procedure up to this point is termed Gaussian elimination.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).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..
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
