Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 21
Текст из файла (страница 21)
Let’s call it τj , j = 1, . . . , n. This array will keep a record of the columninterchanges that we do as we do them, so that in the end we will be able to identify theoutput. Initially, we put τj = j for j = 1, . . . , n. If at a certain moment we are aboutto interchange, say, the pth column and the qth column, then we will also interchange theentries τp and τq . At all times then, τj will hold the number of the column where the currentjth column really belongs.3.3 Building blocks for the linear equation solver93It must be noted that there is a fundamental difference between the interchange of rowsand the interchange of columns.
An interchange of two rows corresponds simply to listingthe equations that we are trying to solve in a slightly different sequence, but has no effect onthe solutions. On the other hand, an interchange of two columns amounts to renumberingtwo of the unknowns. Hence we must keep track of the column interchanges while we aredoing them, so we’ll be able to tell which unknown is which at output time, but we don’tneed to record row interchanges.At the end of the calculation then, the output arrays will have to be shuffled. The readermight want to think about how do carry out that rearrangement, and we will return to itin section 3.6 under the heading of “to unscramble the eggs”.The next item to consider is that we would like our program to be able to solve notjust one system Ax = b, but several systems of simultaneous equations, each of the formAx = b, where the left-hand sides are all the same, but the right-hand sides are different.The data for our program will therefore be an m by n + p matrix whose first n columns willcontain the coefficient matrix A and whose last p columns will be the p different right-handside vectors b.Why are we allowing several different right sides? Some of the main customers for ourprogram will be matrices A whose inverses we want to calculate.
To find, say, the firstcolumn of the inverse of A we want to solve Ax = b, where b is the first column of theidentity matrix. For the second column of the inverse, b would be the second column ofthe identity matrix, and so on. Hence, to find A−1 , if A is an n × n matrix, we must solven systems of simultaneous equations each having the same left-hand side A, but with ndifferent right-hand side vectors.It is convenient to solve all n of these systems at once because the reduction that weapply to A itself to bring it into reduced echelon form is useful in solving all n of thesesystems, and we avoid having to repeat that part of the job n times. Thus, for matrixinversion, and for other purposes too, it is very handy to have the capability of solvingseveral systems with a common left-hand side at once.The next point concerns the linear array τ that we are going to use to keep track of thecolumn interchanges.
Instead of storing it in its own private array, it’s easier to adjoin itto the matrix A that we’re working on, as an extra row, for then when we interchange twocolumns we will automatically interchange the corresponding elements of τ , and therebyavoid separate programming.This means that the full matrix that we will be working with in our program will be(m + 1) × (n + p) if we are solving p systems of m equations in n unknowns with p righthand sides. In the program itself, let’s call this matrix C. So C will be thought of as beingpartitioned into blocks of sizes as shown below:C=A: m×nτ : 1×nRHS : m × p .(3.3.1)0: 1×pNow a good way to begin the writing of a program such as the general-purpose matrix94Numerical linear algebraanalysis program that we now have in mind is to consider the different procedures, ormodules, into which it may be broken up.
We suggest that the individual blocks that weare about to discuss should be written as separate subroutines, each with its own clearlydefined input and output, each with its own documentation, and each with its own localvariable names. They should then be tested one at a time, by giving them small, suitabletest problems. If this is done, then the main routine won’t be much more than a string ofcalls to the various blocks.1.
Procedure searchmat(C,r,s,i1,j1,i2,j2)This routine is given an r × s array C, and two positions in the matrix, say [i1 , j1 ] and[i2 , j2 ]. It then carries out a search of the rectangular submatrix of C whose northwestcorner is at [i1 , j1 ] and whose southeast corner is at [i2 , j2 ], inclusive, in order to find anelement of largest absolute value that lives in that rectangle. The subroutine returns thiselement of largest magnitude, as big, and the row and column in which it lives, as iloc,jloc.Subroutine searchmat will be called in at least two different places in the main routine.First, it will do the search for the next pivot element in the southeast rectangle. Second, itcan be used to determine if the equations are consistent by searching the right-hand sidesof equations r + 1, .
. . , m (r is the pseudorank) to see if they are all zero (i.e., below ourtolerance level).2. Procedure switchrow(C,r,s,i,j,k,l)The program is given an r × s matrix C, and four integers i, j, k and l. The subroutineinterchanges rows i and j of C, between columns k and l inclusive, and returns a variablecalled sign with a value of −1, unless i = j, in which case it does nothing to C and returnsa +1 in sign.3. Procedure switchcol(C,r,s,i,j,k,l)This subroutine is like the previous one, except it interchanges columns i and j of C,between rows k and l inclusive. It also returns a variable called sign with a value of −1,unless i = j, in which case it does nothing to C and returns a +1 in sign.The subroutines switchrow and switchcol are used during the forward solution in theobvious way, and again after the back solution has been done, to unscramble the output(see procedure 5 below).4.
Procedure pivot(C,r,s,i,k,u)Given the r × s matrix C, and three integers i, k and u, the subroutine assumes thatCii = 1. It then stores Cki in the local variable tm sets Cki to zero, and reduces row k ofC, in columns u to s, by doing the operation Ckq := Ckq − t ∗ Ciq for q = u, .
. . , s.The use of the parameter u in this subroutine allows the flexibility for economical operation in both the forward and back solution. In the forward solution, we take u = i + 1and it reduces the whole row k. In the back solution we use u = n + 1, because the rest ofrow k will have already been reduced.5. Procedure scale(C,r,s,i,u)3.3 Building blocks for the linear equation solver95Given an r × s matrix C and integers i and u, the routine stores Cii in a variable calledpiv.
It then does Cij := Cij /piv for j = u, . . . , s and returns the value of piv.6. Procedure scramb(C,r,s,n)This procedure permutes the first n rows of the r × s matrix C according to the permutation that occupies the positions Cr1 , Cr2 , . . . , Crn on input.The use of this subroutine is explained in detail in section 3.6 (q.v.). Its purpose is torearrange the rows of the output matrix that holds a basis for the kernel, and also the rows ofthe output matrix that holds particular solutions of the give system(s).
After rearrangementthe rows will correspond to the original numbering of the unknowns, thereby compensatingfor the renumbering that was induced by column interchanges during the forward solution.This subroutine poses some interesting questions if we require that it should not use anyadditional array space beyond the input matrix itself.7. Procedure ident(C,r,s,i,j,n,q)This procedure inserts q times the n × n identity matrix into the n × n submatrix whoseNorthwest corner is at position [i, j] of the r × s matrix C.Now let’s look at the assembly of these building blocks into a complete matrix analysisprocedure called matalg(C,r,s,m,n,p,opt,eps). Input items to it are:• An r ×s matrix C (as well as the values of r and s), whose Northwest m×n submatrixcontains the matrix of coefficients of the system(s) of equations that are about to besolved. The values of m and n must also be provided to the procedure.
It is assumedthat r = 1 + max(m, n). Unless the inverse of the coefficient matrix is wanted, theNortheast m × p submatrix of C holds p different right-hand side vectors for whichwe want solutions.• The numbers r, s, m, n and p.• A parameter option that is equal to 1 if we want an inverse, equal to 2 if we wantto see the determinant of the coefficient matrix (if square) as well as a basis for thekernel (if it is nontrivial) and a set of p particular solution vectors.• A real parameter eps that is used to bound roundoff error.Output items from the procedure matalg are:• The pseudorank r• The determinant det if m = n• An n × r matrix basis , whose columns are a basis for the kernel of the coefficientmatrix.• An n × p matrix partic, whose columns are particular solution vectors for the givensystems.96Numerical linear algebraIn case opt = 1 is chosen, the procedure will fill the last m columns and rows of C withan m × m identity matrix, set p = n = m, and proceed as before, leaving the inverse matrixin the same place, on output.Let’s remark on how the determinant is calculated.
The reduction of the input matrixto echelon form in the forward solution phase entails the use of three kinds of operations.First we divide a row by a pivot element. Second, we multiply a row by a number and addit to another row. Third, we exchange a pair of rows or columns.The first operation divides the determinant by that same pivot element. The secondhas no effect on the determinant. The third changes the sign of the determinant, at anyrate if the rows or columns are distinct. At the end of the forward solution the matrix isupper triangular, with 1’s on the diagonal, hence its determinant is clearly 1.What must have been the value of the determinant of the input matrix? Clearly itmust have been equal to the product of all the pivot elements that were used during thereduction, together with a plus or minus sign from the row or column interchanges.Hence, to compute the determinant, we begin by setting det to 1. Then, each time a newpivot element is selected, we multiply det by it.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
