Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 90
Текст из файла (страница 90)
He found that after one step of refinement the Cray Y-MP solutionwas almost as accurate as that from the Cray 2 without refinement.25.4. CompilersSome of the pitfalls a compiler writer should avoid when implementing floating point arithmetic are discussed by Farnum [364, 1988]. The theme of hispaper is that programmers must be able to predict the floating point operations that will be performed when their codes are compiled; this may require,for example, that the compiler writer forgoes certain “optimizations”.
In particular, compilers should not change groupings specified by parentheses. Forexample, the two expressions(1.0E-30 + 1.0E+30) - 1.0E+301.0E-30 + (1.0E+30 - 1.0E+30)will produce different answers on many machines. Farnum explains thatCompiler texts and classes rarely address the peculiar problemsof floating-point computation, and research literature on the topicis generally confined to journals read by numerical analysts, notcompiler writers. Many production-quality compilers that are excellent in other respects make basic mistakes in their compilationof floating-point, resulting in programs that produce patently absurd results, or worse, reasonable but inaccurate results.25.5. Determining Properties of Floating Point ArithmeticClever algorithms have been devised that attempt to reveal the propertiesand parameters of a floating point arithmetic.
The first algorithms were published by Malcolm [724, 1972] (see Problem 25.3). These have been improvedand extended by other authors. Kahan’s paranoia code carries out detailedinvestigations of a computer’s floating point arithmetic; there are Basic, C,Modula, Pascal, and Fortran versions, all available from netlib. In addition tocomputing the arithmetic parameters, paranoia tries to determine how wellthe arithmetic has been implemented (so it can be regarded as a program totest a floating point arithmetic—see the next section). Karpinski [645, 1985]gives an introduction to paranoia for the layman, but the best documentationfor the code is the output it produces.Cody has a Fortran routine macher for determining 13 parameters associated with a floating point arithmetic system [222, 1988] (an earlier versionwas published in the book by Cody and Waite [228, 1980]).
Cody notes some498S OFTWARE I SSUESINF LOATING P OINT A RITHMETICstrange behaviour of compilers and says that he was unable to make his coderun correctly on the Cyber 205. A routine based on machar is given in Numerical Recipes [842, 1992, §20.1].LAPACK contains a routine xLAMCH for determining machine parameters.Because of the difficulty of handling the existing variety of machine arithmeticsit is over 850 lines long (including comments and the subprograms it calls).Programs such as machar, paranoia, and xLAMCH are difficult to write; forexample, xLAMCH tries to determine the overflow threshold without invokingoverflow. The Fortran version of paranoia handles overflow by requiring theuser to restart the program, after which checkpointing information previouslywritten to a file is read to determine how to continue.Fortran 90 contains environmental inquiry functions, which for REAL arguments return the precision (PRECISION23), exponent range (RANGE), machineepsilon (EPSILON), largest number (HUGE), and so on, corresponding to thatargument [749, 1990].
The values of these parameters are chosen by the implementor to best fit a model of the machine arithmetic due to Brown [150,1981] (see §25.7.4). Fortran 90 also contains functions for manipulating floating point numbers: for example, to set or return the exponent or fractionalpart (EXPONENT, FRACTION, SET_EXPONENT) and to return the spacing of thenumbers having the same exponent as the argument (SPACING).25.6. Testing a Floating Point ArithmeticHow can we test whether a particular implementation of floating point arithmetic is correct? It is impractical to test a floating point operation with allpossible arguments because there are so many of them—about 1019 in IEEEdouble precision arithmetic, for example! Special techniques are thereforeneeded that test with a limited number of arguments.A package FPV [784, 1986] from NAG Ltd.
attempts to verify experimentally that a floating point arithmetic has been correctly implemented according to its specification. FPV must be supplied by the user with the arithmeticparameters β, t, emin, emax, and with the rounding rule; it then attempts toverify that the arithmetic conforms to these parameters by probing for errors.The tests cover the basic arithmetic operations (including square root) butnot the elementary functions. FPV, which is available in both Fortran 77 andPascal versions, adopts an approach used in an earlier program FPTST [906,1981], [907, 1984] by Schryer of AT&T Bell Laboratories: it tests arithmeticoperations with a limited set of operands that are regarded as being mostlikely to reveal errors.
This approach is based on the premise that errors aremost likely to occur as “edge effects”, induced by operands close to a discontinuity or boundary of the floating point number system (such as a power of23Intrinsic function names are shown in parentheses.25.7 P ORTABILITY499the base β).FPV and FPTST have both revealed implementation errors in floatingpoint arithmetics on commercial computers. Errors detected include multiplication and negation producing unnormalized results, x * y differing from(–x) * (–y), and the product of two numbers greater than 1 being less than 1!Wichmann [1078, 1992] suggests that it was probably revelations such as thesethat led the UK Ministry of Defence to issue an interim standard prohibitingthe use of floating point arithmetic in safety-critical systems.A Fortran package ELEFUNT by Cody and Waite contains programs totest the elementary functions [221, 1982], [228, 1980]; the package is availablefrom netlib.
It checks identities such as cos(x ) = cos(x/3)(4 cos2(x/3) – 1),taking care to choose arguments x for which the errors in evaluating theidentities are negligible. A package CELEFUNT serves the same purposefor the complex elementary functions [224, 1993]. Tang [988, 1990] developstable-driven techniques for testing the functions exp and log.25.7. PortabilitySoftware is portable if it can be made run on different systems with justa few straightforward changes (ideally, we would like to have to make nochanges, but this level of portability is often impossible to achieve). Sometimesthe word “transportable” is used to describe software that requires certainwell-defined changes to run on different machines. For example, Demmel,Dongarra, and Kahan [290, 1992] describe LAPACK as “a transportable wayto achieve high efficiency on diverse modern machines”, noting that to achievehigh efficiency the BLAS need to be optimized for each machine.
A goodexample of a portable collection of Fortran codes is LINPACK. It containsno machine-dependent constants and uses the PFORT subset of Fortran 77;it uses the level-1 BLAS, so, ideally, optimized BLAS would be used on eachmachine.25.7.1. Arithmetic ParametersDifferences between floating point arithmetics cause portability problems.First, what is meant by REAL and DOUBLE PRECISION in Fortran varies greatlybetween machines, as shown by the following table:Second, for a given level of precision, u, the various arithmetic parameterssuch as base, unit roundoff, largest and smallest machine numbers, and the500S OFTWARE I SSUESINFLOATING P OINT A RITHMETICtype of rounding, can all vary. Some possible ways for a code to acquire theseparameters are as follows.(1) The parameters are evaluated and embedded in the program in PARAMETER and DATA statements.
This is conceptually the simplest approach, butit is not portable.(2) Functions are provided that return the machine parameters. Bell Laboratories’ PORT library [405, 1978] has three functions:R1MACH(k) 1 < k < 5D1MACH(k) 1 < k < 5I1MACH(k) 1 < k < 16real parameters,double precision parameters,integer parameters.R1MACH returns the underflow and overflow thresholds, the smallest and largestrelative spacing (β –t, β 1-t respectively), and log10 β, where β is the baseand t the number of digits in the mantissa.
I1MACH returns standard input,output and error units and further floating point information, such as β, t,and the minimum and maximum exponents emin and emax. These functionsdo not carry out any computation; they contain DATA statements with theparameters for most common machines in comment lines, and the relevantstatements have to be uncommented for a particular machine. This approachis more sensible for a library than the previous one, because only these threeroutines have to be edited, instead of every routine in the library.The NAG Library takes a similar approach to PORT.
Each of the 18routines in its X02 chapter returns a different machine constant. For example,XO2AJF returns the unit roundoff and X02ALF returns the largest positivefloating point number. These values are determined once and for all when theNAG library is implemented on a particular platform using a routine similarto paranoia and machar, and they are hard coded into the Chapter X02routines.(3) The information is computed at run-time, using algorithms such asthose described in §25.5.25.7.2. 2 x 2 Problems in LAPACKLAPACK contains 10 or so auxiliary routines for solving linear equation andeigenvalue problems of dimension 2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
