Nash - Compact Numerical Methods for Computers (523163), страница 19
Текст из файла (страница 19)
These algorithms are complicated by ALGOL'S lack of double-length arithmetic, necessitating calls to machinecode procedures. (This comment applies to ALGOL-60, not ALGOL-68.)By and large I have avoided scaling within my programs because of the greatdifficulty of making any reliable general recommendations. Indeed, given any twonon-singular diagonal matrices D and E, the system of equationsDAEE - 1x = D b(6.33)has the same solution x as the equationsAx = b.(2.2)In scaling the equations by row multiplication we are adjusting D, which adjuststhe pivot selection. It is often recommended that the scaling factors in D bechosen to equilibrate the matrix A, that is, so thatmax |( DA)ij| = 1jfor i = 1, 2, . .
. , n(6.34)where for the moment E is assumed to be a unit matrix. This is simply a dose ofcommon sense which attempts to avoid arithmetic involving numbers widelydifferent in magnitude. However, as Dahlquist and Björck (1974, pp 181–3)point out, the scaling E-1 of the solution x can frustrate our efforts to stabilise thecomputation. Furthermore, optimal scaling depends on knowledge of the matrixA-1, which is not known. They therefore suggest E be chosen to reflect ‘theimportance of the unknowns’. This statement is suitably amorphous to coverwhatever situations arise, so I shall venture the opinion that the magnitudes of thesolution elementsy = E- 1x(6.35)should be roughly equivalent. That is to say, the variables in the problem at handshould be measured in units which give the expected solution elementsLinear equations—a direct approach81approximately the same size.
Is this worth the bother? I can only add that I rarelyscale sets of equations unless there is some very obvious and natural way to do it.Similar comments apply to iterative improvement of a computed solution(Dahlquist and Björck 1974, pp 183-5, Bowdler et al 1966). Given a computedsolution(6.36)if(6.37)then a triangular decomposition of A permits solution ofAc = LRc = r(6.38)for c and computation of a new solution(6.39)The process can be repeated until the computed solution has converged, but invirtually all cases the improvement of the solution occurs in the first application of(6.36) to (6.39).
A similar algorithm can be used to improve least-squaressolutions (Dahlquist and Björck 1974, pp 204–5). Unfortunately these improvement procedures are dependent on accurate computation of residuals, anddouble-precision arithmetic must be called into play. As the systems for which thisbook is designed often lack this feature, one may fall into the habit of not usingiterative improvement of linear-equation or least-squares solutions.Even when computing in an environment where double-length arithmetic isavailable, I generally do not bother to employ it.
Personally, very little of mywork has concerned numbers which have arisen out of precise measurement. Infact, my clients are often only sure of the first one or two digits of their data, sothat it is unnecessary to provide an extremely accurate solution, though it isimportant in many cases to identify near-linear dependencies (hence singularities)by means of a technique such as the singular-value decomposition.So far in this section I have mentioned only techniques of which I do notgenerally make use. To finish, then, consider the one excursion from algorithms 5and 6 which has diverted me from time to time.
This is the purely organisationalquestion of how the information in the working array A should be organised andaccessed. In the algorithms as presented, I have chosen to perform interchangesexplicitly and store the coefficient matrix and right-hand sides together in a singletwo-dimensional array. The choice of a single working array with solutionsoverwriting the right-hand sides b I feel to be the sensible one for small-computerimplementations. The choice of method for accessing the elements of this array isless simple. Besides the direct, two-dimensional method which has been used, it ispossible to perform pivot interchanges implicitly if the pivot positions are saved,for instance in an integer vector q so that the ith pivot is stored in A[q[i,i].
Thusif the algorithm is started so thatq[i] = ifor i = 1, 2, . . . , n(6.40)Compact numerical methods for computers82then Gauss elimination and back-substitution can be carried out exactly as inalgorithms 5 and 6 if every array reference is made with A[ , ] replaced byA[q[ , ]. However, a simplification occurs in the interchange step 3, which canbe replaced by a simple interchange of the row indices.
That is, at step j, if thepivot is in row q [k] q[j], or k j, then the indices are simply interchanged ratherthan the entire rows. However, all array access operations are complicated.Some overall increases in efficiency may be obtained if we take over thecompiler or interpreter function in accessing two-dimensional arrays.
That is, westore the working array A which is m = (n + p) by n in a single vector a of mnelements. We can do this columnwise, so thatA[ i ,j] = a[n * (j – 1) + i](6.41)A[i, j] = a[m * (i – 1) + j].(6.42)or row-wise, so thatThese translations offer some simplifications of the elimination and backsubstitution algorithms. In fact, the row-wise form (6.41) is more useful forelimination where the index of an element is simply incremented to proceedacross a row of the coefficient matrix. For back-substitution, we need to formmatrix-vector products which oblige us to access array elements by marchingsimultaneously across rows and down columns. Implicit pivoting is also possiblewith a one-dimensional storage scheme.
This adds just one more item to thosefrom which a method must be selected.It is probably clear to my readers that I have already decided that simplest isbest and intend to stick with algorithms 5 and 6. My reasons are as follows.(i) Despite the elegance of implicit pivoting, the extra index vector and theprogram code needed to make it work are counter to the spirit of a compactalgorithm.(ii) The implicit interchange only gains in efficiency relative to the direct methodif an interchange is needed; this is without counting the overhead which arrayaccess via q implies.
But in many instances very few interchanges are required andthe whole discussion then boils down to an argument over the likely number ofinterchanges in the problem set to be solved.(iii) In coding Gauss elimination with back-substitution and the Gauss-Jordanreduction with various of the above choices, S G Nash and I (unpublishedwork) found that the implicit pivoting methods were surprisingly prone to ‘bugs’which were difficult to discover.
This applied particularly to the one-dimensionalstorage forms. Most of these errors were simple typographical errors in entry ofthe code. Since it is hoped the algorithms in this book will prove straightforwardto implement, only a direct method has been included.6.4. COMPLEX SYSTEMS OF EQUATIONSConsider the system of equations (where i = (–1)½)( Y+iZ) (u +i v) = g + i h.(6.43)Linear equations—a direct approach83Separating these into real and imaginary components gives the real equations(6.44)Yu – Zv = g(6.45)Yv + Zu = hwhich is a set of linear equations (2.22) with(6.46)(6.47)and(6.48)This is how complex systems of linear equations can be solved using realarithmetic only. Unfortunately the repetition of the matrices Y and Z in (6.46)means that for a set of equations of order n, 2n2 storage locations are usedunnecessarily.
However, the alternative is to recode algorithms 5 and 6 to takeaccount of the complex arithmetic in (6.43). Bowdler et al (1966) give ALGOLprocedures to perform the Crout variant of the elimination for such systems ofequations, unfortunately again requiring double-length accumulation.6.5. METHODS FOR SPECIAL MATRICESThe literature contains a number of methods for solving special systems ofequations. For instance, several contributions in Wilkinson and Reinsch (1971)deal with band matrices, that is, those for whichAij = 0if | i– j | > k(6.49)for some k.
Thus if k = 1, the matrix is tridiagonal. While these methods areundoubtedly useful and save memory, I have not included them in this monograph because I feel any serious user with enough special problems to warrant amethod tailored to the task is likely to find and implement one. Others may onlyfind too many special cases tedious or bewildering. Thus no discussion of bandedor other special forms is given, though the user should be alert to triangular formssince it is very wasteful of effort to apply Gauss elimination to a lower-triangularmatrix when simple forward-substitution will suffice. Likewise, no treatment isincluded of the various iteration methods for the systems of equations arisingfrom partial differential equations (see Varga 1962). It should be pointed out,however, that the Givens’ reduction can often be organised to take advantage ofpatterns of zeros in matrices.
Even as it stands, algorithm 3 is quite efficient forsuch problems, since very little work is done when zero elements are encounteredand no pivot interchanges are made.The only special form which will be considered is a symmetric positive definitematrix. Chapter 7 deals with a decomposition of such a matrix useful for solvingspecial sets of linear equations. Chapter 8 discusses a very compact algorithm forinverting such a matrix in situ, that is, on top of itself.Previous Home NextChapter 7THE CHOLESKI DECOMPOSITION7.1. THE CHOLESKI DECOMPOSITIONWhen the matrix A is symmetric and positive definite (see discussion below fordefinition), it is possible to perform (without pivoting) the symmetric decompositionA = LLT = RT R(7.1)whereL = RT(7.2)is a lower-triangular matrix. In fact, it is also possible to perform Gauss elimination in a symmetric fashion for symmetric positive definite matrices withoutpivoting for stability (see Dahlquist and Björck 1974, pp 162-4).The Choleski algorithm (Wilkinson 1965) is derived directly from (7.1), that is(7.3)Note that the summation runs only from 1 to the minimum of i and j due to thetriangular nature of L.
Thus we have(7.4)so thatL11 = (A11)½.(7.5)Ai1 = LilL11(7.6)L i1 = A i1 /L 11 .(7.7)Furthermoreso that we obtainConsider now the mth column of L which is defined for i > m by(7.8)with the diagonal element determined first by setting i = m. It is straightforward tosee that every element in the right-hand side of equation (7.8) comes fromcolumns 1, 2, . . . , (m – 1) of L or from column m of A.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
