Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 103
Текст из файла (страница 103)
With n = 3 and almost any starting data, the backward error can easily bemade of order 1, showing that the method is unstable. However, the backward erroris found to be roughly of orderso the method may have a backward errorbound proportional to(this is an open question).25.1. (b) An optimizing compiler might convert the test xp1 > 1 to x+1 > 1 andthen to x > 0. (For a way to stop this conversion in Fortran, see the solution toProblem 25.3.) Then the code would compute a number of order 2 e min instead of anumber of order 2–t.576S OLUTIONSTOP ROBLEMS25.3.
The algorithm is based on the fact that the positive integers that can beexactly represented are 1,2, ..., β t andIn the interval [β t, β t+l] the floating point numbers are spaced β apart. This intervalkmust contain a power of 2, a = 2 . The first while loop finds such an a (or, rather,the floating point representation of such an a) by testing successive powers 2i to seeif both 2i and 2i + 1 are representable.
The next while loop adds successive powersof 2 until the next floating point number is found; on subtracting a the base β isproduced. Finally t is determined as the smallest power of β for which the distanceto the next floating point number exceeds 1.The routine can fail for at least two reasons. First, an optimizing compiler mightsimplify the test while (a+l) - a == 1 to while 1 == 1, thus altering the meaningof the program.
Second, in the same test the result (a+l) might be held in anextra length register and the subtraction done in extra precision. The computationwould then reflect this higher precision and not the intended one. We could tryto overcome this problem by saving the result a+l in a variable, but an optimizingcompiler might undo this effort by dispensing with the variable and storing theresult in a register.
In Fortran, the compiler’s unwanted efforts can be nullified bya test of the form while foo (a+l) - a == 1, where foo is a function that simplyreturns its argument. The problems caused by the use of extra length registers werediscussed by Gentleman and Marovich [437, 1974]; see also Miller [759, 1984, §2.2].25.4. (This description is based on Schreiber [903, 1989].) The random numbergenerator in matgen repeats after 16384 numbers [668, 1981, p. 19]. The dimensionn = 512 divides the period of the generator (16384 = 512 x 32), with the effect thatthe first 32 columns of the matrix are repeated 16 times (512 = 32 x 16), so thematrix has this structure:B = rand(512,32);A = [B, B, B, B, B, B, B, B, B, B, B, B, B, B, B, B];Now consider the first 32 steps of Gaussian elimination.
We apply 32 transformationsto A that have the effect, in exact arithmetic, of making B upper trapezoidal. Infloating point arithmetic, they leave a residue of small numbers (about u10-7in size) in rows 33 onward. Because of the structure of the matrix, identical smallnumbers occur in each of the 15 blocks of A to the right of the first. Thus, theremaining (512 – 32) x (512 – 32) submatrix has the same block structure as A (butwith 15 block columns). Hence this process repeats every 32 steps:after 32 steps the elements drop to10 -7 ;after 64 steps the elements drop to10 –14 ;after 96 steps the elements drop to10 -21 ;after 128 steps the elements drop to10 –28 ;after 160 steps the elements drop to10 –35 ;after 192 steps the elements would drop to 10 –42 ,S OLUTIONSTOP ROBLEMS577but that is less than the under flow threshold.
The actual pivots do not exactly matchthe analysis, which is probably due to rank deficiency of one of the submatricesgenerated. Also, underflow occurs earlier than predicted, apparently because twosmall numbers (both O( 10–21 ) ) are multiplied in a saxpy operation.25.5.
Let si and tt denote the values of s and t at the start of the ith stage of thealgorithm. Then(A.15)If |xi | > ti then the algorithm uses the relationso with si+l = (ti /xi )2si + 1 and ti+l = |xi|, (A.15) continues to hold for the nextvalue of i. The same is true trivially in the case |xi| < ti .This is a one-pass algorithm using n divisions that avoids overflow except, poswhich can overflow only if ||x||2 does.sibly, on the final stage in forming25.7. (a) The matrix A := I + Y2 = I – YTY is symmetric positive semidefinitesince ||Y||2 < 1, hence it has a (unique) symmetric positive semidefinite square rootX satisfying X2 = A = I + Y2.
The square root X is a polynomial in A (see,for example, Higham [532, 1987]) and hence a polynomial in Y; therefore X and Ycommute, which implies that (X+ Y)T(X + Y) = (X – Y)(X + Y) = X2 – Y2 = I,as required.25.8. |x|/3 could underflow to zero, or become unnormalized, whennot.doesPrevious HomeAppendix BSingular Value Decomposition,M-MatricesThe singular value decomposition is a matrix factorization which can produceapproximations to large arrays.Cryptanalysts is the task of breaking coded messages.In this paper, we present an unusual merger of the two in which thesingular value decomposition may aid the cryptanalyst in discoveringvowels and consonants in messages coded in certainvariations of simple substitution ciphers.— CLEVE B.
MOLER and DONALD MORRISON,Singular Value Analysis of Cryptograms (1983)then each of the following 50 conditions is equivalent to the statement:“A is a nonsingular M-matrix”.— ABRAHAM BERMAN and ROBERT J. PLEMMONS,Nonnegative Matrices in the Mathematical Sciences (1994)579Next580SINGULAR VALUE DECOMPOSITION, M-MATRICESIn this appendix we define the singular value decomposition and M-matrices.B.1. Singular Value DecompositionAny matrixhas a singular value decomposition (SVD)are both unitary.Theare the singular values of A and the columns of U and V are the leftand right singular vectors of A, respectively.For more details on the SVD see Golub and Van Loan [470, 1989, pp. 7173], Stewart and Sun [954, 1990, pp.
30-34], and Horn and Johnson [580, 1985,§7.3], [581, 1991, §3.1]. The history of the SVD is described by Stewart [950,1993] and Horn and Johnson [581, 1991, $3.0].B . 2 . M-Matricesand all theA matrixis an M-matrix if aij < 0 for alleigenvalues of A have nonnegative real part. This is one of many equivalentdefinitions [94, 1994, Chap. 6].An M-matrix may be singular. A particularly useful characterization of anonsingular M-matrix is a nonsingular matrixfor which aij < 0for alland A-1 has nonnegative elements (written as A–1 > 0).For more information on ill-matrices see Berman and Plemmons [94, 1994]and Horn and Johnson [581, 1991, §2.5].Previous HomeAppendix CAcquiring SoftwareCaveat receptor .
. .Anything free comes with no guarantee!— JACK DONGARRA and ERIC GROSSE, Netlib mail header581NextA CQUIRING S OFTWARE582In this appendix we provide information on how to acquire software mentionedin the book. First, we describe some basic aspects of the Internet.C.
1. InternetA huge variety of information and software is available over the Internet, theworldwide combination of interconnected computer networks. The location ofa particular object is specified by a URL, which stands for “(Uniform ResourceLocator”. Examples of URLS arehttp://www.netlib.erg/index.htmlftp://ftp.netlib. orgThe first example specifies a World Wide Web server (http = hypertexttransfer protocol) together with a file in hypertext format (html = hypertext markup language), while the second specifies an anonymous ftp site.
Inany URL, the site address may, optionally, be followed by a filename thatspecifies a particular file. For more details about the Internet see on-line information, or one of the many books on the subject, such as Krol [674, 1994].C.2. NetlibNetlib is a repository of freely available mathematical software, documents,and databases of interest to the scientific computing community [316, 1987],[151, 1994]. It includeslresearch codes,lgolden oldies (classic programs that are not available in standard libraries),lthe collected algorithms of the ACM,lprogram libraries such as EISPACK, LINPACK, LAPACK, and MINPACKlback issues of NA-Digest, a weekly digest for the numerical analysiscommunity,ldatabases of conferences and performance data for a wide variety ofmachines.Netlib also enables the user to download technical reports from certain institutions, to download software and errata for textbooks, and to search “theSIAM membership list and a “white pages” database.Netlib can be accessed in several ways.C.3 M ATLAB583lOver the World Wide Web.
The URL ishttp://www.netlib.erg/index.htmllVia Xnetlib, an X Windows application that provides interactive’ accessto netlib.lBy anonymous ftp. The URL isftp://ftp.netlib.orglBy electronic mail. For an introduction and master index, send a oneline email message as follows:mail netlib@ornl. govsend indexNetlib is mirrored at various sites throughout the world.C.3. M ATLABM ATLAB is a commercial program sold by The MathWorks, Inc. It runson a variety of platforms.
The MathWorks maintains a collection of usercontributed M-files, which is accessible over the Internet.For information contactThe MathWorks, Inc.24 Prime Park WayNatick, MA 01760-1500USATel: 5086531415Fax: 5086532997email: info@mathworks.comURL: ftp://ftp.mathworks.comURL: http://www.mathworks.comC.4. NAG Library and FTN90 CompilerThe Numerical Algorithms Group (NAG) produces a variety of software products.
Relevant to this book are the FTN90 compiler and the NAG Library, alarge numerical program library available in several programming languages.For information contactNAG Ltd.Wilkinson HouseA CQUIRING S OFTWARE584Jordan Hill RoadOxford, OX2 8DRUKTel: +44 1865 511245Fax: +44 1865 310139email: infodesk@nag.co.ukURL: http://www.nag.co.uk:70/NAG has subsidiaries and distributors, whose addresses can be obtainedfrom the above sources.Previous Home NextAppendix DProgram LibrariesSince the programming is likely to be themain bottleneck in the use of an electronic computerwe have given a good deal of thought to thepreparation of standard routines of considerable generality for themore important processes involved in computation.By this means we hope to reduce the time takento code up large-scale computing problems,by building them up, as it were,from prefabricated units.— J. H.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
