Nash - Scientific Computing with PCs (523165), страница 21
Текст из файла (страница 21)
InBASIC, the DATA statements require space that can be released by such a restructuring.We recommend small experiments to decide which approach is most economical of space. The resultingcode may be less simple or readable, but necessity takes priority over elegance, though some changes tobranching may save code yet render the structure more transparent. Figure 7.4.1 shows a small example.The truncate function INT replaces a loop structure, but we must add a PRINT to ensure subsequentPRINTs start on a new line when N is not a whole multiple of 5.Figure 7.4.1Simple example of restructuring.
Both program segments will print a vector V of Nelements with 5 elements printed on each line.A. 1000 REM PRINT VECTOR 5 ELEMENTS PER101010201030104010501060LINEFOR I=1 TO N STEP 5FOR J=I TO I+4PRINT V(J),NEXT JPRINTNEXT ISize: 99 bytes7.5B. 1000 REM PRINT VECTOR 5 ELEMENTS PER10101020103010401050LINEFOR I=1 TO NPRINT V(I),IF 5*INT(I/5)=I THEN PRINTNEXT IPRINTSize 92 bytesReprogrammingAn extreme case of restructuring is the complete reprogramming of the computation. This is not nearlythe admission of defeat it may be taken to be. First, programs always grow and change with the tasks towhich they are put. There are nearly always parts that can and should be simplified, especially after"patches" are made.
Second, reprogramming allows variable and array names to be rationalized. Third,the experience gained with the program can be incorporated into improvements in the rewritten version.Patches and fixes can be cleanly made a part of the whole. Finally, documentation can be brought up todate.Obviously, reprogramming is not a task to be undertaken without due consideration.
Nevertheless it maybe better than hours of frustrating work trying to restructure a much altered code, possibly withoutup-to-date documentation.7: SIZE DIFFICULTIES7.655ChainingChaining is an "old" technique for handling size problems. We break up the task into pieces that can berun sequentially in a "chain". Some language processors allow for automatic chaining. That is, an executingprogram may command the PC to load and begin executing another program.
While this technique stillexists in some compilers and interpreters, it has largely disappeared as a formal tool. Users can, however,arrange for similar functionality in a variety of ways, in particular BATch command files in MS-DOS orprograms that execute other programs, such as hypertext systems. We use "chaining" to refer to any formof this idea.The primary requirements for successful chaining are:•Suitable breakpoints in the task and•A mechanism for passing data from one program to another.The first of these requirements is dictated by the task at hand. Note that almost all problems have at leastthree common parts: initialization and input, processing, and reporting of results.
Rarely can we not splita task into pieces.The second requirement, that of data transfer, is sometimes provided for in the language processor.Compilers and linkers that permit overlays are performing such a job. We avoid such techniques in favorof simplicity. It is quite easy to write needed data to a file in one program, then read it again with the nextprogram in the chain. The chaining tasks can then be accomplished by BATch commands or similarmethod.Such "start-from-file" facilities in a program can be very helpful if we must carry out very long programruns. Our programs can be set up to periodically output a file with a standard name, e.g.,PROGSTRT.DTA, that is always read when the program starts and which contains the status of theprogram.
If it is necessary to interrupt the program to perform other tasks, or if power failure or otherbreakdown upset the progress of execution, we can simply restart without losing the whole of our work.This is especially useful for programs that use iterative algorithms.7.7Change of MethodThe last way to overcome size difficulties is to change the method by which a solution is sought. This isreally formulation and choice of method, which will be discussed again in Chapters 12 and 13. Normallywe will not change the method unless there are strong motivations.
In our experience, running out ofspace is such a motivation.When changing methods for this reason, the main criterion is to reduce the overall memory space. Wemake our point by using the example of function minimization methods, for which Table 7.7.1 showssome typical working storage requirements for some well-known methods (Nash J C and Nash S G, 1988).Other factors, such as efficiency in finding a minimum and need for gradient information will influenceour choice, but we can see that the Conjugate Gradient (cg) and quasi-Newton (qn) method both usefunction and gradient information, yet cg only requires a few vectors of storage, while qn needs a matrix.If the methods are not too dissimilar in performance, then cg can save us space as n increases.
GenerallyTable 7.7.1 illustrates some signposts to reduced storage requirements we should watch for:•The use of vectors instead of matrices;•Iteration rather than direct solution;•Short code length, even if there are many loops;•Use of data implicitly rather than explicitly, e.g., a way to compute a matrix-vector product where wedo not use the matrix itself;56Copyright © 1984, 1994 J C & M M NashSCIENTIFIC COMPUTING WITH PCsNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 Canada•Copy for:Dr. Dobb’s JournalSegmentation of data, so we use it piecemeal from storage.The divisions between restructuring, reprogramming and change of method are not always clear. Anexample of space-saving that merges the ideas is Algorithm 15 from Nash J C (1990d).
This code isdesigned to solve the generalized matrix eigenproblemA x = e B x(7.5.1)for symmetric A and B, with B positive definite. That is, we want to obtain all the possible sets ofeigenvalues e and eigenvectors x that satisfy Equation 7.5.1. Several techniques are available but most ofthose that solve for all the eigensolutions first decompose the matrix B, then develop a standardsymmetric matrix eigenvalue problem. In the algorithm in question, we decided to use a Jacobi matrixeigenvalue routine to accomplish both the decomposition and the resulting standard eigenproblem.
Thusthe code length was considerably reduced.Table 7.7.1 Working storage requirements for some function minimization methods. Only the leadingterm in n, the order of the problem, is given for storage requirements.MethodStoragerequirementsHooke and JeevesOrthogonal searchNelder MeadNewtonTruncated NewtonQuasi-NewtonConjugate GradientSteepest Descent2nn2n2n2/29n to 13nn2/23n2nFunctionrequirementsffff,f,f,f,f,Work per stepg, HggggO(n)O(n)O(n)2O(n )2O(n )2O(n )O(n)O(n)f = function, g = gradient, H = Hessian7.8Data Storage Mechanisms for MatricesMatrices often have a special structure that can be exploited to save memory space when solvingparticular problems. Matrices are so important in scientific computing that we include here a few tacticsand structures that are useful.
We will use v[ ] and a[ , ] to represent one and two dimensional storagearrays in computer programming languages.An LU decomposition of a square matrix A, (used in elimination methods for solving linear equations)can often overwrite the original array. This is accomplished easily if one of L or U must always have unitdiagonal elements. Suppose(7.8.1)Lii = 1 for i = 1,2,. . .,nThis can be "remembered" by the program and requires no storage. The other elements of thedecomposition are stored as follows:(7.8.2)Uij is stored in a[i,j](7.8.3)Uij = 0 for j < i (upper triangle)(7.8.4)Lij is stored in a[i,j] for j < i (strict lower triangle) sincefor j >= i, since7: SIZE DIFFICULTIES(7.8.5)57Lij = 0 for j >= iAny triangular matrix can be stored in a linear array (vector). There are several ways of doing this, butone simple one (Nash, 1990d, page 82) is to use a row-wise ordering, which for a lower triangular formhas the following equivalence:Figure 7.8.1Array storage ordering[2 D] Array(1,1)(2,1)(2,2)(3,1)(3,2)(4,1)(4,2)(3,3)(4,3)(4,4)[1 D] Array (vector)12 34 5 67 8 9 10The relationship can be written(7.8.6)v[i*(i-1)/2+j] =AijNote that we could also use a column-wise enumeration1234(7.8.7)5678910In which case the correspondence is(7.8.8)v[j*(j-1)/2+i] =AijFor upper triangular matrices we must transpose these rules.
Programmers making use of such techniquesshould be very careful to make sure that the correspondences are always correct. It is also worth notingthat for matrices less than, say, 15 by 15, the extra program code will offset the storage saving of n(n-1)/2matrix elements, where n is the order of the matrix. The break-even value of n can be found by testing.These triangular schemes are appropriate for symmetric matrices, where(7.8.9)Aij = AjiBanded matrices, which have(7.8.10)Aij = 0 for i-j > k-1where k is the band width, can be stored as a sequence of diagonals.
Consider the case where we havea symmetric band matrix with band width k = 5. Thus we have(7.8.11)Aij = 0 for i-j > 4Since the matrix is symmetric, the subdiagonal elements Aij for j < i are equal to the superdiagonals. We can store the data in a rectangular array a, with the following correspondences(7.8.12)a[i,1] =Aii for i = 1,2,. . .,n(7.8.13)a[i,2] = Ai,i+1 = Ai+1,i for[note: a[n,2] = 0](7.8.14)a[i,3] = Ai,i+2 = Ai+2,i for i = 1,2,.[note: a[n,3] = a[n-1,3] = 0]i = 1,2,. .
.,(n-1). .,(n-2)Aji58Copyright © 1984, 1994 J C & M M NashSCIENTIFIC COMPUTING WITH PCsNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 CanadaCopy for:Dr. Dobb’s Journal2In this case, we store only 3n elements instead of n = n*n elements. Clever programmers probably caneven use the three unused elements of a, i.e., a[n,2], a[n,3] and a[n-1,3].Hermitian matrices are complex matrices defined so thatH = H*(7.8.15)where * represents the conjugate transpose. That is, if Re(x) is the real part of the complex variable x andIm(x) is its imaginary part,Hij = Re(Hji) - Im(Hji)(7.8.16)But the equality imposed means that the real part ofanti-symmetric (Nash J C 1974).
Thus,H is symmetric, while the imaginary part isRe(Hij) = Re(Hji)Im(Hij) = - Im(Hji)(7.8.17)(7.8.18)Note that Im(Hii) = - Im(Hii), so thatIm(Hii) =(7.8.19)0We therefore can store the real part of the Hermitian matrix in one triangular matrix, and the imaginarypart in a triangular matrix with no diagonal.
These can be put back together in a square array as follows.(7.8.20)Re(Hij)= a[i,j]for j <== a[j,i]for j(7.8.21)Im(Hij)= a[j,i]for j <= a[i,j]for=0fori>iij>ij=iThere are so many mechanisms for storing sparse matrices that it is futile to try to survey them. We findtwo simple mechanisms useful.•The array is stored as three linear arrays (vectors) that only mention non-zero elements: r[j] gives therow index of an element, c[j] gives the column index of that element and v[j] gives its value, for j =1,2,. . .,E, where E is the total number of elements. Note that the triples of numbers (r[j], c[j], v[j]) canbe written in any order, but that some orderings may be helpful for particular calculations such asmatrix multiplication.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
