Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 23
Текст из файла (страница 23)
The matrix C (which is called byname in the procedure, which means that the input matrix is altered by the procedure) willcontain a basis for the kernel of the coefficient matrix in columns psrank + 1 to n, and pparticular solution vectors, one for each input right-hand side, in columns n + 1 to n + p.Two new sub-procedures are called by this procedure, namely scaler and pivotr.These are called immediately before the action of scale or pivot, respectively, and theirmission is to update the R matrix in accordance with equations (3.4.2) or (3.4.4) to takeaccount of the impending scaling or pivoting operation .Exercises 3.41.
Break off from the complete algorithm above, the forward solution process. State itformally as algorithm forwd, list its global variables, and describe precisely its effecton them. Do the same for the backwards solution.2. When the program runs, it gives the solutions and their roundoff error estimates.Work out an elegant way to print the answers and the error estimates. For instance,there’s no point in giving 12 digits of roundoff error estimate. That’s too much. Justprint the number of digits of the answers that can be trusted. How would you dothat? Write subroutine prnt that will carry it out.102Numerical linear algebra3. Suppose you want to re-run a problem, with a different set of random numbers inthe roundoff matrix initialization.
How would you do that? Run one problem threeor four times to see how sensitive the roundoff estimates are to the choice of randomvalues that start them off.4. Show your program to a person who is knowledgeable in programming, but who isnot one of your classmates.
Ask that person to use your program to solve some set ofthree simultaneous equations in five unknowns.Do not answer any questions verbally about how the program works or how to use it.Refer all such questions to your written documentation that accompanies the program.If the other person is able to run your program and understand the answers, awardyourself a gold medal in documentation. Otherwise, improve your documentation andlet the person try again. When successful, try again on someone else.5. Select two vectors f , g of length 10 by choosing their elements at random. Form the10 × 10 matrix of rank 1 whose elements are fi gj . Do this three times and add theresulting matrices to get a single 10 × 10 matrix of rank three.Run your program on the coefficient matrix you just constructed, in order to see ifthe program is smart enough to recognize a matrix of rank three when it sees one, byhalting the forward solution with pseudorank = 3.Repeat the above experiment 50 times, and tabulate the frequencies with which yourprogram “thought” that the 10 × 10 matrix had various ranks.3.5Operation countWith any numerical algorithm it is important to know how much work is involved in carryingit out.
In this section we are going to estimate the labor involved in solving linear systemsby the method of the previous sections.Let’s recognize two kinds of labor: arithmetic operations, and other operations, both asapplied to elements of the matrix. The arithmetic operations are +, −, ×, ÷, all lumpedtogether, and by other operations we mean comparisons of size, movement of data, andother operations performed directly on the matrix elements. Of course there are many“other operations,” not involving the matrix elements directly, such as augmenting counters,testing for completion of loops, etc., that go on during the reduction of the matrix, but thetwo categories above represent a good measure of the work done.
We’re not going to includethe management of the roundoff error matrix R in our estimates, because its effect wouldbe simply to double the labor involved. Hence, remember to double all of the estimates ofthe labor that we are about to derive if you’re using the R matrix.We consider a generic stage in the forward solution where we have been given m equations in n unknowns with p right-hand sides, and during the forward solution phase we havejust arrived at the [i, i] element.3.5 Operation count103The next thing to do is to search the Southeast rectangle for the largest element, Therectangle contains about (m − i) ∗ (m − i) elements.
Hence the search requires that manycomparisons.Then we exchange two rows (n+p−i operations), exchange two columns (m operations)and divide a row by the pivot element (n + p − i arithmetic operations).Next, for each of the m−i−1 rows below the ith, and for each of the n+p−i elements ofone of those rows, we do two arithmetic operations when we do the elementary row operationthat produces a zero in the ith column.
This requires, therefore, 2(n + p − i)(m − i − 1)arithmetic operations.For the forward phase of the solution, therefore, we countAf =rX{2(n + p − i)(m − i − 1) + (n + p − i)}(3.5.1)i=1arithmetic operations altogether, where r is the pseudorank of the matrix, because theforward solution halts after the rth row with only zeros below.The non-arithmetic operations in the forward phase amount toNf =rX{(m − i)(n − i) + (n + p − i) + m}.(3.5.2)i=1Let’s leave these sums for a while, and go to the backwards phase of the solution.
Wedo the columns in reverse order, from column r back to 1, and when we have arrived at ageneric column j, we want to create zeroes in all of the positions above the 1 in the [j, j]position.To do this we perform the elementary row operation−−→ := −−→ − Aij ∗ −−→row(i)row(i)row(j)(3.5.3)to each of the j − 1 rows above row j. Let i be the number of one of these rows. Then,−→ are acted upon by the elementary row operation above?exactly how many elements of −row(i)−→ that lie in columns 1 through j − 1 are unaffected, becauseCertainly the elements in −row(i)−−→only zero elements are in row(j)below them, thanks to the forward reduction process.−→ that lie in columns j + 1 through r are unaffected,Furthermore, the elements in −row(i)for a different reason. Indeed, any such entry is already zero, because it lies above an entryof 1 in some diagonal position that has already had its turn in the back solution (rememberthat we’re doing the columns in the sequence r, r − 1, .
. . , 1). Not only is such an entry zero,−→but it remains zero, because the entry of −row(j)below it is also zero, having previouslybeen deleted by the action of a diagonal element below it.−−→Hence in row(i),the elements that are affected by the elementary row operation (3.5.3)are those that lie in columns j, n − r, . . . , n, n + 1, .
. . , n + p (be sure to write [or modify!]the program so that the row reduction (3.5.3) acts only on those columns!). We have now104Numerical linear algebra−→shown that exactly N + p + r − 1 entries of each row above −row(j)are affected (note thatthe number is independent of j and i), soAb =rX(n + p − r + 1)(j − 1)(3.5.4)j=1arithmetic operations are done during the back solution, and no other operations.It remains only to do the various sums, and for this purpose we recall thatNXi=1NXi=N (N + 1)2(3.5.5)N (N + 1)(2N + 1)i =6i=12Then it is straightforward to find the total number of arithmetic operations from Af +AbasArith(m, n, p, r) =r3r2− (2m + n + p − 5) + ((n + p)(2m − 5/2) − m + 1/3)r62(3.5.6)and the total of the non-arithmetic operations from Nf asNonArith(m, n, p, r) =r3r2r− (m + n) + (6mn + 3m + 3n + 6p − 2) .326(3.5.7)Let’s look at a few important special cases.
First, suppose we are solving one systemof n equations in n unknowns that has a unique solution. Then we have m = n = r andp = 1. We find that2Arith(n, n, 1, n) = n3 + O(n2 )(3.5.8)3where O(n2 ) refers to some function of n that is bounded by a constant times n2 as n growslarge. Similarly, for the non-arithmetic operations on matrix elements we find 13 n3 + O(n2 )in this case.It follows that a system of n equations can be solved for about one third of the price,in terms of arithmetic operations, of one matrix multiplication, at least if matrices aremultiplied in the usual way (did you know that there is a faster way to multiply twomatrices? We will see one later on).Now what is the price of a matrix inversion by this method? Then we are solving nsystems of n equations in n unknowns, all with the same left-hand side.
Hence we haver = m = n = p, and we find thatArith(n, n, n, n) =13 3n + O(n2 ).6(3.5.9)Hence we can invert a matrix by this method for about the same price as solving 3.25systems of equations! At first glance, it may seem as if the cost should be n times as great3.6 To unscramble the eggs105because we are solving n systems.
The great economy results, of course, from the commonleft-hand sides.The cost of the non-arithmetic operations remains at 13 n3 + O(n2 ).If we want only the determinant of a square matrix A, or want only the rank of A thenwe need to do only the forward solution, and we can save the cost of the back solution. Weleave it to the reader to work out the cost of a determinant, or of finding the rank.3.6To unscramble the eggsNow we have reached the last of the issues that needs to be discussed in order to plan acomplete linear equation solving routine, and it concerns the rearrangement of the outputso that it ends up in the right order.During the operation of the forward solution algorithm we found it necessary to interchange rows and columns so that the largest element of the Southeast rectangle was broughtinto the pivot position.
As we mentioned previously, we don’t need to keep a record of therow interchanges, because they correspond simply to solving the equations in a differentsequence. We must remember the column interchanges that occur along the way though,because each time we do one of them we are, in effect, renumbering the unknowns.To remember the column interchanges we glue onto our array C an additional row, justfor bookkeeping purposes. Its elements are called τj , j = 1, . .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
