Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 89
Текст из файла (страница 89)
We give four examples to show the benefits of IEEE arithmetic on software.Suppose we wish to evaluate the dual of the vector p- n o r mthat is, the g-norm, where p –l + q –1 = 1. In MATLAB notation we simplywrite norm (x, 1/(1-l/p)), and the extreme cases p = 1 andcorrectlyyield q = 1/(1 – 1/p) = and 1 in IEEE arithmetic. (Note that the formulaq = p/ (p – 1) would not work, because inf/inf evaluates to a NaN.)The following code finds the biggest element in an array a (1: n), and elegantly solves the problem of choosing a “sufficiently negative number” withwhich to initialize the variable max:max = –inffor j = 1:nif aj > maxmax = ajendendAny unknown or missing elements of the array a could be assigned NaNs(a(j) = NaN) and the code would (presumably) still work, since a NaN compares as unordered with everything.Consider the following continued fraction [630, 1981]:3r(x) = 7 –,1x - 2 –10x - 7 +x - 22x - 3which is plotted over the range [0,4] in Figure 25.1.
Another way to write the25.1 E XPLOITING IEEE A RITHMETIC493Figure 25.1. Rational function r.rational function r(x) is in the formin which the polynomials in the numerator and denominator are written inthe form in which they would be evaluated by Homer’s rule. Examiningthese two representations of r, we see that the continued fraction requiresless arithmetic operations to evaluate (assuming it is evaluated in the obvious“bottom up” fashion) but it incurs division by zero at the points x = 1:4,even though r is well behaved at these points.
However, in IEEE arithmetic revaluates correctly at these points because of the rules of infinity arithmetic.For x = 1077, the rational form suffers overflow, while the continued fractionevaluates correctly to 7.0; indeed, in IEEE arithmetic the continued fraction isimmune to overflow. Figure 25.2 shows the relative errors made in evaluatingr in double precision on an equally spaced grid of points on the range [0,4](many of the errors for the continued fraction are zero, hence the gaps inthe error curve); clearly, the continued fraction produces the more accuratefunction values.
The conclusion is that in IEEE arithmetic the continuedfraction representation is the preferred one for evaluating r.The exception handling provided by IEEE arithmetic can be used to simplify and speed up certain types of software, provided that the programmingenvironment properly supports these facilities.
Recall that the exceptional494S OFTWARE I SSUESINFLOATING P OINT A RITHMETICFigure 25.2. Error in evaluating rational function r. Solid line: continued fraction,dotted line: usual rational form.operations underflow, overflow, divide by zero, invalid operation, and inexactdeliver a result, set a flag, and continue. The flags (one for each type of exception) are “sticky”, remaining set until explicitly cleared, and implementationsof the IEEE standard are required to provide a means to read and write theflags. Unfortunately, compiler writers have been slow to make exception handling, and some of IEEE arithmetic’s other unusual features, available to theprogrammer.
Kahan has commented that “the fault lies partly with the IEEEstandard, which neglected to spell out standard names for the exceptions, theirflags, or evenAs an example of how exception handling can be exploited, consider theuse of the LAPACK norm estimator (Algorithm 14.4) to estimate ||A–1 ||1 .The algorithm requires the solution of linear systems of the form Ax = b andATy = c, which is typically done with the aid of an LU factorization PA =LU. Solution of the resulting triangular systems is prone to both overflow anddivision by zero. In LAPACK, triangular systems are solved in this context notwith the level-2 BLAS routine xTRSV but with a routine xLATRS that containselaborate tests and scalings to avoid overflow and division by zero [18, 1991].In an IEEE arithmetic environment a simpler and potentially faster approachis possible: solve the triangular systems with xTRSV and after computing eachsolution check the flags to see whether any exceptions occurred.
Provided thatsome tests are added to take into account overflow due solely to ||A||1 being25.2 S UBTLETIESOFFLOATING P OINT A RITHMETIC495tiny, the occurrence of an exception allows the conclusion that A is singularto working precision; see Demmel and Li [298, 1994] for the details.For the condition estimator that uses exception handling to be uniformlyfaster than the code with scaling it is necessary that arithmetic with NaNs,infinities , and subnormal numbers be performed at the same speed as conventional arithmetic.
While this requirement is often satisfied, there are machinesfor which arithmetic on NaNs, infinities, and subnormal numbers is betweenone and three orders of magnitude slower than conventional arithmetic [298,1994, Table 2]. Demmel and Li compared the LAPACK condition estimation routines xxxCON with modified versions containing exception handlingand found the latter to be up to 6 times faster and up to 13 times slower,depending on the machine and the matrix.Appendices to the IEEE standards 754 and 854 recommend 10 auxiliaryfunctions and predicates to support portability of programs across differentIEEE systems.
These include nextafter(x,y), which returns the next floatingpoint number to x in the direction of y, and scalb(x,n), which returns x x β nwithout explicitly computing β n, where n is an integer and β the base. Notall implementations of IEEE arithmetic support these functions. Portable Cversions of six of the functions are presented by Cody and Coonen [226, 1993].25.2. Subtleties of Floating Point ArithmeticThe difference x – y of two machine numbers can evaluate to zero even whenbecause of underflow. This makes a test such asf = f / (x - y)endunreliable. However, in a system that supports gradual underflow (such asIEEE arithmetic) x – y always evaluates as nonzero when xy, as is easyto verify. On several models of Cray computer (Cray 1, 2, X-MP, Y-MP,and C90) this code could fail for another reason: they compute f/(x – y) asf * (1/(x – y)) and 1/(x – y) could overflow.It is a general principle that one should not test two floating point numbersfor equality, but rather use a test such as if |x – y| < tol (there are exceptions,such as Algorithm 2 in §1.14.1).
Of course, skill may be required to choosean appropriate value for tel. A test of this form would correctly avoid thedivision in the example above when underflow occurs.496S OFTWARE I SSUESINFLOATING P OINT A RITHMETICTable 25.1. Results from Cholesky factorization.Bits1286464646464ule-29le-17le-16le-16le-15le-15ComputerCray 2IBM 3090Convex 220IRISCray 2Cray Y-MPDisplacement.447440341.447440344.44744033 9.44744033 9.447440303.44743610625.3. Cray PeculiaritiesCarter [188, 1991] describes how a Cray computer at the NASA Ames Laboratories produced puzzling results that were eventually traced to propertiesof its floating point arithmetic. Carter used Cholesky factorization on a CrayY-MP to solve a sparse symmetric positive definite system of order 16146resulting from a finite element model of the National Aerospace Plane.
Theresults obtained for several computers are shown in Table 25.1, where “Displacement” denotes the largest component of the solution and incorrect digitsare set in italics and underlined. Since the coefficient matrix has a conditionnumber of order 1011, the errors in the displacements are consistent with theerror bounds for Cholesky factorization.Given that both the Cray 2 and the Cray Y-MP lack guard digits, it is notsurprising that they give a less accurate answer than the other machines witha similar unit roundoff.
What is surprising is that, even though both machinesuse 64-bit binary arithmetic with a 48-bit mantissa, the result from the CrayY-MP is much less accurate than that from the Cray 2. The explanation(diagnosed over the telephone to Carter by Kahn, as he scratched on the backof an envelope) lies in the fact that the Cray 2 implementation of subtractionwithout a guard digit produces more nearly unbiased results (average errorzero), while the Cray Y-MP implementation produces biased results, causingfl(x – y) to be too big if x > y > 0. Since the inner loop of the Choleskyalgorithm contains an operation of the form a i i = a i i –the errors inthe diagonal elements reinforce as the factorization proceeds on the CrayY-MP, producing a Cholesky factor with a diagonal that is systematicallytoo large.
Similarly, since Carter’s matrix has off-diagonal elements that aremostly negative or zero, the Cray Y-MP produces a Cholesky factor withoff-diagonal elements that are systematically too large and negative. For thelarge value of n used by Carter, these two effects in concert are large enoughto cause the loss of a further two digits in the answer.An inexpensive way to improve the stability of the computed solution is to25.4 COMPILERS497use iterative refinement in fixed precision (see Chapter 11). This was tried byCarter.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
